Anaïs Crestetto
email: anais.crestetto@univ-nantes.fr

Nicolas Crouseilles
email: nicolas.crouseilles@inria.fr

Giacomo Dimarco
email: giacomo.dimarco@unife.it

Mohammed Lemou
email: mohammed.lemou@univ-rennes1.fr

Asymptotically complexity diminishing schemes (ACDS) for kinetic equations in the diffusive scaling

Keywords: Mathematics Subject Classification: 65M06, 35B25, 82C80, 82D10, 41A60 Kinetic equations, Diffusion scaling, Asymptotic preserving schemes, Asymptotically complexity diminishing schemes, Micro-macro decomposition, Monte Carlo methods

published or not. The documents may come

Asymptotically complexity diminishing schemes (ACDS) for kinetic equations in the diffusive scaling

Anaïs Crestetto, Nicolas Crouseilles, Giacomo Dimarco,

Introduction

Kinetic equations are commonly used to describe several phenomena arising in physics such as rarefied gas [START_REF] Bird | Molecular gas dynamics and direct simulation of gas flows[END_REF][START_REF] Dimarco | Numerical methods for kinetic equations[END_REF], neutron transport [START_REF] Case | Linear Transport Theory[END_REF], radiative transfer [START_REF] Chandrasekhar | Radiative Transfer[END_REF], plasmas [START_REF] Birsdall | Plasma Physics Via Computer Simulation[END_REF] or electron flow in semiconductors [START_REF] Markowich | The Stationary Semiconductor Device Equations[END_REF]. In the description of these problems, one often encounters different spatiotemporal scales which are due to the coexistence of equilibrium and non equilibrium phenomena. This makes the problem hard to solve from the numerical point of view: the full kinetic model has to be considered everywhere to capture the correct physics. This means that the numerical method should deal with the high dimensionality issue and with different scales at the same time. Then, whenever a direct discretization is considered, strong constraints on the time step are required which involves huge computational cost. In addition, passing from the microscopic to the macroscopic description may change the nature of the equations and the numerical method should be able to follow this transition. This is the specific case considered in this work, in which equations pass from a hyperbolic to a parabolic structure [START_REF] Bardos | Fluid dynamic limit of kinetic equations II: Convergence proofs for the Boltzmann equations[END_REF]. Consequently, the characteristic speeds of the hyperbolic system grow to infinity as the scaling parameter tends to zero [START_REF] Bardos | Fluid dynamic limit of kinetic equations II: Convergence proofs for the Boltzmann equations[END_REF][START_REF] Cercignani | The Boltzmann equation and its applications[END_REF] causing severe CFL restrictions in standard explicit numerical methods. Unfortunately, the alternative consisting in using macroscopic models can not be an option since macroscopic models are not able to capture non equilibrium effects. Thus, a natural idea consists in performing a domain decomposition approach and use the suitable model in the appropriate region [START_REF] Bal | Coupling of transport and diffusion models in linear transport theory[END_REF][START_REF] Bal | Coupling of transport and diffusion models in linear transport theory[END_REF][START_REF] Buet | Asymptotic preserving scheme and numerical methods for radiative hydrodynamic models[END_REF][START_REF] Degond | A smooth transition model between kinetic and diffusion equations[END_REF][START_REF] Degond | A multiscale kinetic-fluid solver with dynamic localization of kinetic effects[END_REF][START_REF] Degond | Fluid simulations with localized Boltzmann upscaling by Direct Simulation Monte-Carlo[END_REF][START_REF] Golse | A domain decomposition analysis for a two-scale linear transport problem[END_REF][START_REF] Jin | Uniformly accurate diffusive relaxation schemes for multiscale transport equations[END_REF][START_REF] Jin | Numerical study of a domain decomposition method for a two-scale linear transport equation[END_REF][START_REF] Klar | Asymptotic induced domain decomposition methods for kinetic and drift diffusion semiconductor equations[END_REF][START_REF] Klar | A numerical method for kinetic semiconductor equations in the drift diffusion limit[END_REF][START_REF] Klar | Transition from kinetic theory to macroscopic fluid equations: a problem for domain decomposition and a source for new algorithms[END_REF][START_REF] Klar | Boundary layers and domain decomposition for radiative heat transfer and diffusion equations: applications to glass manufacturing processes[END_REF][START_REF] Tallec | Coupling Boltzmann and Navier-Stokes by half-fluxes[END_REF][START_REF] Tiwari | A particle-particle hybrid method for kinetic and continuous equations[END_REF]. This strategy is somehow optimal in terms of efficiency since the expensive kinetic model is used only when necessary, but it requires to connect the macro and micro models which is not an easy task at both mathematical and numerical levels, even if efficient approaches have been recently developed [START_REF] Filbet | A hierarchy of hybrid numerical methods for multi-scale kinetic equations[END_REF][START_REF] Xiong | A Hierarchical Uniformly High Order DG-IMEX Scheme for the 1D BGK Equation[END_REF]. Another way to handle such multiscale nature of the phenomena is to use asymptotic preserving (AP) methods [START_REF] Boscarino | Implicit-Explicit Runge-Kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit[END_REF][START_REF] Boscarino | A unified IMEX Runge-Kutta approach for hyperbolic systems with multiscale relaxation[END_REF][START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF][START_REF] Dimarco | Asymptotic Preserving Monte Carlo methods for transport equations in the diffusive limit[END_REF][START_REF] Lafitte | Asymptotic-preserving projective integration schemes for kinetic equations in the diffusion limit[END_REF][START_REF] Jin | Efficient Asymptotic-Preserving (AP) schemes for some multiscale kinetic equations[END_REF][START_REF] Ren | An Asymptotic-Preserving Monte Carlo method for the Boltzmann equation[END_REF][START_REF] Klar | An asymptotic-induced scheme for non stationary transport equations in the diffusive limit[END_REF][START_REF] Krycki | Asymptotic preserving numerical schemes for a nonclassical radiation transport model for atmospheric clouds[END_REF][START_REF] Larsen | Asymptotic solution of neutron transport problems for small mean free paths[END_REF][START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF][START_REF] Lemou | Micro-macro schemes for kinetic equations including boundary layers[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF][START_REF] Naldi | Numerical schemes for kinetic equations in diffusive regimes[END_REF] which enable to overcome the numerical stiffness and to use time and space steps which are independent of the stiffness parameters which characterize the fast scales. However, while this approach permits to solve the problem related to the choice of the small time and space steps, it does not overcome the cost related to the solution of the kinetic model even in the equilibrium regions where a much less complex asymptotic model could be used.

In this work, we propose a numerical method which tries to take the best of the above two approaches: domain decomposition and asymptotic preservation. The method is based on a hybrid Monte Carlo strategy [START_REF] Caflisch | Monte Carlo and Quasi-Monte Carlo Methods[END_REF][START_REF] Liu | Monte Carlo strategies in scientific computing[END_REF] for the underlying kinetic equation which follows the ideas introduced in [START_REF] Dimarco | A Fluid Solver Independent Hybrid method for multiscale kinetic equations[END_REF][START_REF] Dimarco | Hybrid multiscale methods II. Kinetic equations[END_REF][START_REF] Degond | The moment guided Monte Carlo method[END_REF] and more recently in [START_REF] Crouseilles | Asymptotic preserving and time diminishing schemes for rarefied gas dynamic[END_REF][START_REF] Crestetto | A particle micro-macro decomposition based numerical scheme for collisional kinetic equations in the diffusion scaling[END_REF][START_REF] Crestetto | Kinetic/Fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equations using particles[END_REF] in the context of gas dynamics. The problem addressed is the solution of a kinetic type equation under a diffusive scaling which in the limit gives rise to a parabolic type equation. More in details, the equation we consider is the following radiative transport equation (RTE)

∂ t f + 1 ε v • ∇ x f = 1 ε 2 (ρM -f ), f (t = 0, x, v) = f 0 (x, v), (1.1)
where f (t, x, v) represents the density of particles in phase space with x ∈ Ω ⊂ R dx , v ∈ V = R dv , d x ∈ N * is the dimension in space, d v ∈ N * in velocity space and ε the scaling parameter. The density ρ(t, x) of the particles is given by

ρ(t, x) = V f (t, x, v)dv, while M (v) = 1 (2π) dv/2 exp - |v| 2 2 (1.2)
is the so-called absolute Maxwellian. Finally f 0 (x, v) is a given initial condition and we consider periodic boundary conditions in space. To achieve our goal in this context, we proceed as follows: we couple a Monte Carlo method for the solution of the kinetic model with a finite volume method for the solution of the macroscopic one [START_REF] Dimarco | A Fluid Solver Independent Hybrid method for multiscale kinetic equations[END_REF][START_REF] Dimarco | Hybrid multiscale methods II. Kinetic equations[END_REF][START_REF] Degond | The moment guided Monte Carlo method[END_REF]. However, instead of solving directly the kinetic model, we reformulate the original problem by using a micro-macro decomposition strategy [START_REF] Liu | Boltzmann equation: micro-macro decompositions and positivity of shock profiles[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF][START_REF] Bennoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics[END_REF][START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF][START_REF] Crestetto | Kinetic/Fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equations using particles[END_REF][START_REF] Crestetto | A particle micro-macro decomposition based numerical scheme for collisional kinetic equations in the diffusion scaling[END_REF]. This permits to derive an equivalent set of equations which expresses the time evolution of the macroscopic part ρ(t, x)M (v) plus the time evolution of the perturbation g(t, x, v). Our strategy will be shown to have a low computational complexity in the fluid regions, which then makes it competitive with the computational complexity of the domain decomposition methods in these regions. In other words, our method has almost a fluid complexity in the fluid regions. Moreover, it has even lower complexity in the kinetic regions compared to other approaches (like AP or domain decomposition methods), since only the deviation from equilibrium has to be discretized with the Monte Carlo method, not the whole distribution f . Successively, we derive our numerical method by using an asymptotic preserving approach for the microscopic part of the solution which overcomes the parabolic stiffness and freezes the characteristic speeds in terms of the scaling parameter.

Since the kinetic equation which describes the evolution of the perturbation is solved by a particle approach, the result is a scheme in which the computational cost depends on the magnitude of the perturbation. The larger is the fraction of the solution far from equilibrium, the larger is the number of particles used. Moreover, since this fraction changes with time and space, the method automatically adapts itself to optimize the solution in terms of computational cost and reduction of the statistical error. In fact, the part of the solution described by particles is the lowest possible at each instant of time. This means that the proposed method realizes a variance reduction method whose effectiveness depends on the regime studied. As and where the equilibrium is approached, the number of particles in the Monte Carlo method is reduced, which decreases the computational cost enough to make it almost equivalent to the computational cost of a classical numerical method for the diffusion equation. With respect 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 of the problem in the diffusion asymptotics.

The rest of the paper is organized as follows. Section 2 is concerned with the presentation of the model, the derivation of the micro-macro equivalent model and its diffusion asymptotic. In Section 3 we present the asymptotically stable discretization of the microscopic model and a Monte Carlo Lagrangian discretization is proposed. In the same section, the discretization of the macroscopic part by a suitable Eulerian scheme that avoids the parabolic constraint on the time step is discussed. Section 4 is devoted to numerical results and analysis of the performances of the proposed method by treating problems from one up to three spatial dimensions. Finally, some remarks conclude the work in Section 5 together with an outlook on future researches.

Micro-Macro decomposition and the limit diffusion equation

In this section, we introduce the micro-macro decomposition of (1.1) for the distribution f [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF][START_REF] Bennoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics[END_REF][START_REF] Liu | Boltzmann equation: micro-macro decompositions and positivity of shock profiles[END_REF] which permits to separate the time evolution of the equilibrium part from the time evolution of the perturbation part. This reads

f (t, x, v) = ρ(t, x)M (v) + g(t, x, v), ρ(t, x) = f (t, x) (2.1)
with M (v) the Maxwellian equilibrium already introduced, i.e. a Gaussian distribution centered in zero with variance one. The non equilibrium g satisfies by construction for all times

g = 0 where f = V f (v)dv (2.2)
with V = R dv . In the following, we use the notation Πφ = φ M to indicate the orthogonal projection of a function φ into the null space N of the relaxation operator (ρ(t, x)M (v)-f (t, x, v)).

With this notation, the so called microscopic-macroscopic (micro-macro) model reads

     ∂ t ρ + 1 ε v • ∇ x g = 0, ∂ t g + 1 ε [v • ∇ x ρM + (I -Π)(v • ∇ x g)] = - 1 ε 2 g. (2.3)
The above model is equivalent to the original one (1.1) and is obtained as follows. We first plug in (1.1) the decomposition (2.1). This gives

∂ t ρM + ∂ t g + 1 ε v • ∇ x ρM + 1 ε v • ∇ x g = - 1 ε 2 g. (2.4)
Then, we apply (I -Π) with I the identity operator to this equation and use Πg = 0, Π(∂ t g) = 0, Π(vM ) = 0, Π(ρM ) = ρM to get the second equation of (2.3).

On the other hand, applying Π to equation (2.4) gives the first equation of (2.3). Now, when ε goes to zero, one gets from the second equation of (2.3)

g = -ε(v • ∇ x ρM ) + O(ε 2 ).
(2.5)

Then, the macroscopic equation becomes

∂ t ρ -v • ∇ x (v • ∇ x ρM ) = O(ε), which gives ∂ t ρ -∆ x ρ = O(ε).
The limit diffusion equation is then obtained by passing to the limit ε → 0

∂ t ρ -∆ x ρ = 0. (2.6)
3 A new asymptotically complexity diminishing scheme (ACDS)

In this section, we present the new numerical method to solve equation (1.1). The main difficulty related to the derivation of a numerical scheme is due to the fact that the scaling parameter ε can range of several order of magnitude both in time and space. Thus, in the situation in which the scaling parameter is small, the problem becomes stiff and in particular this causes the characteristic speeds to grow to infinity. Among the possible solutions, the class of asymptotic preserving method represents certainly a good choice to tackle the method [START_REF] Dimarco | Asymptotic Preserving Monte Carlo methods for transport equations in the diffusive limit[END_REF][START_REF] Boscarino | A unified IMEX Runge-Kutta approach for hyperbolic systems with multiscale relaxation[END_REF][START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF][START_REF] Crestetto | A particle micro-macro decomposition based numerical scheme for collisional kinetic equations in the diffusion scaling[END_REF], they permit to choose the time step independently of the stiffness of the underlying equation remaining consistent and stable.

Here, we propose an alternative which enjoys the same consistency and stability properties but, in addition, it also permits to reduce the numerical complexity (and then the computational cost) of the problem. Our first step is presented in the following paragraph and follows the strategy of [START_REF] Crestetto | A particle micro-macro decomposition based numerical scheme for collisional kinetic equations in the diffusion scaling[END_REF]. Indeed, we propose a reformulation of the micro-macro system (2.3) which permits to surround the stiffness of the transport term. We then propose a particle method to solve the reformulated microscopic equation.

A time discrete reformulation

The development of a particle-based scheme for equation (1.1) or equivalently for equation (2.3) involves a splitting between the transport and the collision parts. However, while the collision term can be taken implicit, this cannot be done for the transport part. Avoiding the costly inversion of the advection operator, since this transport part is stiff, its exact resolution is faced to the issue of unbounded characteristic speeds. In order to bound the characteristic speeds, we proceed as follows. We first rewrite the micro part of (2.3) by multiplying it by e t/ε 2 . This reads

∂ t (e t/ε 2 g) = - e t/ε 2 ε F (ρ, g) , (3.1)
where F (ρ, g) is given by

F (ρ, g) = v • ∇ x ρM + v • ∇ x g -v • ∇ x g M. (3.2)
We denote now by ∆t > 0 a time step and set t n = n∆t with n ∈ N. By integrating (3.1) on [t n , t n+1 ], one gets

g(t n+1 ) = e -∆t/ε 2 g(t n ) - 1 ε t n+1
t n e -(t n+1 -τ )/ε 2 F (ρ(τ ), g(τ )) dτ, and then up to an error O(∆t 2 )

g(t n+1 ) = e -∆t/ε 2 g(t n ) -ε(1 -e -∆t/ε 2 )F (ρ(t n ), g(t n )) + O(∆t 2 ).
Now, one can observe that the above equation can be rewritten as

g(t n+1 ) -g(t n ) ∆t = e -∆t/ε 2 -1 ∆t g(t n ) -ε 1 -e -∆t/ε 2 ∆t F (ρ(t n ), g(t n )) + O(∆t), (3.3)
which can be recast, up to an additional error of order O(∆t), as

∂ t g(t n ) = e -∆t/ε 2 -1 ∆t g(t n ) -ε 1 -e -∆t/ε 2 ∆t F (ρ(t n ), g(t n )) + O(∆t).
We now define the first-order in time reformulation of (2.3) as

∂ t ρ + 1 ε ∇ x • vg = 0, (3.4)
∂ t g = e -∆t/ε 2 -1 ∆t g -ε 1 -e -∆t/ε 2 ∆t F (ρ, g) , (3.5)
with F (ρ, g) given by (3.2). Let us observe that the evolution equation (3.5) for the perturbation g, unlike the second equation of (2.3), does not contain any stiff term. In particular, the characteristic speeds are bounded with respect to ε. Therefore, the approximated model (3.5) appears to be more suitable for a numerical discretization which is constructed following the characteristics equations. We recall that this strategy bears some similarities with the so-called Implicit Monte Carlo (IMC) approach (see [START_REF] Fleck | An implicit Monte Carlo scheme for calculating time and frequency dependent nonlinear radiation transport[END_REF][START_REF] Long | The iterative thermal emission method: A more implicit modification of IMC[END_REF][START_REF] Gentile | Implicit Monte Carlo Diffusion -An Acceleration Method for Monte Carlo Time-Dependent Radiative Transfer Simulations[END_REF][START_REF] Densmore | Urbatsch A hybrid transport-diffusion Monte Carlo method for frequency-dependent radiative-transfer simulations[END_REF]) which was originally introduced to overcome the stiffness due to the collisional part of the kinetic equation. However, while the IMC method is usually applied to the full unknown f and efficiency is generally lost close to the asymptotic limit, here the Monte Carlo approach is used to solve the perturbation g and both efficiency and accuracy increase when the scaling parameter goes to zero. Finally, let us remark that this reformulation is consistent with the original model (2.3) in the following sense: when ∆t → 0, (3.5) converges towards the original microscopic equation of (2.3) for a fixed ε > 0. Let us finally stress that higher order (in time) reformulations of the microscopic model are possible. We do not discuss this possibility here and we refer to [START_REF] Crestetto | A particle micro-macro decomposition based numerical scheme for collisional kinetic equations in the diffusion scaling[END_REF] for the derivation of a second order approximation of the micro equation of (2.3) which is smooth in ε.

The Monte Carlo discretization for the perturbation equation

This subsection is devoted to the derivation of a particle based numerical scheme for equation (3.5). The particle method is based on a time splitting which separates the advection part from the interaction part. It reads 1. Transport:

∂ t g + ε 1 -e -∆t/ε 2 ∆t v • ∇ x g = 0, (3.6)
2. Interaction:

∂ t g = e -∆t/ε 2 -1 ∆t g -ε 1 -e -∆t/ε 2 ∆t (v • ∇ x ρM -v • ∇ x g M ) .
(3.7)

The type of scheme we consider is based on stochastic approach and it first consists in replacing the continuous function g(t, x, v) by a set of N p (t) particles which occupy random positions x k (t) ∈ R dx with random velocities v k (t) ∈ R dv and weights ω k at time t sampled accordingly to g(t, x, v) as detailed later on. Thus, we consider the following particle approximation of g:

g(t, x, v) ≈ µ(t, x, v) = Np(t) k=1 ω k δ(x -x k (t))δ(v -v k (t)),
in which δ is the usual Dirac distribution and µ(t, x, v) is a so-called empirical distribution. Let us observe that the above choice corresponds to constant weights ω k while the number of particles in the domain may change with time. An alternative particle approach may consist in taking timedependent weights ω k (t) while maintaining constant their number [START_REF] Liu | Boltzmann equation: micro-macro decompositions and positivity of shock profiles[END_REF]. Starting from the empirical distribution, one needs to reconstruct the particle densities in phase space during the time evolution of the solution. Restricting ourselves to the one dimensional case for simplicity, these densities are obtained by introducing a mesh in the phase space and computing the histogram on this mesh in the cell centers x i and v ℓ with corresponding mesh sizes ∆x and ∆v:

g(t, x i , v ℓ ) = 1 ∆x∆v x i +∆x/2 x i -∆x/2 v ℓ +∆v/2 v ℓ -∆v/2 1 • dµ(t, x, v). (3.8)
In practice, we use the following formulae

g(t, x i , v ℓ ) = 1 ∆x∆v k=1,...,Np(t) (x k (t),v k (t))∈C i,ℓ ω k , where C i,ℓ = [x i -∆x/2, x i + ∆x/2] × [v ℓ -∆v/2, v ℓ + ∆v/2].
The extension to the multidimensional case of such reconstruction is straightforward.

It remains to define the weights ω k . This is done as follows. We first define the mass m p of an individual particle representing a portion of the distribution f (t = 0, x, v) as

m p = 1 N p (t = 0) R dx R dv f (t = 0, x, v)dvdx. (3.9)
Successively, we assign to a particle a positive or negative weight. This is a consequence of the fact that the perturbation g is not positive everywhere and this contrasts with the usual particle methods which are applied on the full distribution function f (in this case particles have only positive weights since f is positive by definition). Thus, the weights ω k are chosen either positive ω k = m p if k belongs to the subset of the N p (t) particles used to represent the positive part of the microscopic function g, i.e. k ∈ N + ⊆ 1, N p , or negative in the opposite case ω k = -m p if k ∈ N -⊆ 1, N p : the subset of the particles used for representing the negative part of g. For consistency we have N + ∩ N -= ∅, while the accuracy of the method is related to the choice of the mass m p : smaller is the mass m p , larger is the accuracy of the method. At time t = 0, the number of particles is given by

N 0,± = ± 1 m p R dx R dv g ± (t = 0, x, v)dvdx,
where the positive (respectively the negative) part of g is

g ± = (g ± |g|)/2. (3.10)
Since the integrals of g ± have no reason to be integers we associate to any real number a the integer ⌊a⌋ with probability 1a + ⌊a⌋ and the integer ⌊a⌋ + 1 with probability a -⌊a⌋ where ⌊a⌋ is the integer part of a. All the real expressions defining integers have to be interpreted with the above meaning.

In the sequel, we denote by C i,ℓ ⊂ R dx × R dv the cell centered around the point (x i , v ℓ ) and by N n,+ i,ℓ (resp. N n,- i,ℓ ) the number of positive (resp. negative) particles in C i,ℓ at time t n . Initially, we use in the 1-dimensional framework the following definition

N 0,± i,ℓ = ± 1 m p x i +∆x/2 x i -∆x/2 v ℓ +∆v/2 v ℓ -∆v/2 g ± (t = 0, x, v)dvdx,
and approximate it by N 0,± i,ℓ = ± ∆x∆v mp g ± (t = 0, x i , v ℓ ). The extension of this definition to the multidimensional case is straightforward.

For the positive particles in C i,ℓ , the sampling is then defined as follows: for all k = 1, . . . , N 0,+ i,ℓ , the velocities v 0 k are assigned to v ℓ , the positions x 0 k are uniformly distributed in the cell C i defined for

d x = 1 by C i = [x i -∆x/2, x i + ∆x/2]
(this notation is generalized in the multidimensional case), and the weights ω k = m p . The same strategy is used for negative particles, expects the weights that are defined as ω k = -m p .

Moreover, we introduce

N n,± i = ℓ N n,± i,ℓ and N n i = N n,+ i + N n,- i
the total number of particles in C i at time t n , n ≥ 0, and finally N n = i N n i . Thanks to the above described setting and notations, we are now ready to detail our algorithm step by step. Thus, we start with an initial sampling of the N p (t = 0) particles (x 0 k , v 0 k ) which approximates the initial micro unknown g(t = 0, x, v), as described previously. Then at each time step ∆t, we follow the splitting procedure presented at the beginning of this subsection.

1. We solve the transport part by shifting particles

x n+1 k = x n k + ε(1 -e -∆t/ε 2 )v n k . (3.11)
The velocities of the particles do not change in this step. This defines an intermediate empirical distribution from which an intermediate particle density gn (x, v) can be recovered by the formulae

gn (x, v) ≃ μn (x, v) = N n k=1 ω k δ(x -x n+1 k )δ(v -v n k ) (3.12)
and then used for solving the next step. In particular, this gives the values gn (x i , v ℓ ) in the cell C i,ℓ . In the 1-dimensional case, the formulae writes

gn (x i , v ℓ ) = 1 ∆x∆v k=1,...,N n p (x n+1 k ,v n k )∈C i,ℓ ω k (3.13)
and is easily extended in the multidimensional case.

Since particles have moved in this transport part, we denote by Ñ n,+ i (resp. Ñ n,i

) the number of positive (resp. negative) particles in C i after the transport step, and we define

Ñ n i = Ñ n,+ i + Ñ n,- i .
2. We solve the interaction part by reporting ∂ t g ∼ g n+1 -g n ∆t into (3.7):

g n+1 = e -∆t/ε 2 gn + (1 -e -∆t/ε 2 )ε -v • ∇ x ρ n M + ∇ x • vg n M , (3.14)
where the function gn is defined by (3.13). In the following we give details about how this step is solved with stochastic particles.

One can observe that g n+1 in (3.14) is written as a linear combination of gn and a given analytic function in v. At the Monte Carlo level, the above formulae can be interpreted in the following way:

• with probability e -∆t/ε 2 , the distribution g n+1 remains to be gn ,

• with probability (1e -∆t/ε 2 ), the distribution g n+1 is replaced by a new distribution given by

P n (x, v) = ε -v • ∇ x ρ n M + ∇ x • vg n M .
In practice, this means that

• a fraction e -∆t/ε 2 of particles is kept with unchanged velocity and position. More precisely, in each cell C i we keep e -∆t/ε 2 Ñ n i particles. The particles are picked and discarded randomly since there is no dependence of the coefficient e -∆t/ε 2 on the velocity.

• One samples a fraction (1-e -∆t/ε 2 )M n,± i particles from the distribution P n,± (x i , v)

x,v P n,± dxdv where the positive and negative parts are defined like in (3.10).

More in details, to samples the M n,± i particles, we consider a piecewise constant discretization of P n,± (x, v), where gradients ∇ x ρ n and ∇ x • vg n are approximated by central differences. In practice, in each cell C i,ℓ we create M n,± i,ℓ particles with weights ±m p . In the 1-dimensional case, for example, we have that

M n,± i,ℓ = 1 m p x i +∆x/2 x i -∆x/2 v ℓ +∆v/2 v ℓ -∆v/2 ±P n,± (x, v)dvdx, is approximated by M n,± i,ℓ = ± ∆x∆v mp P n,± (x i , v ℓ ).
To these particles, we assign velocities v n k = v ℓ , positions x n k uniformly distributed in the cell C i , and weights ω k = ±m p . Moreover, since only a fraction (1e -∆t/ε 2 ) of the sampled particles is needed for the next time step, we directly sample this smaller amount of particles instead of sampling the full set and then discarding a part. This automatically allows to diminish the number of used particles as ε → 0. Let us finally remark that a simple sampling technique is used here but other more efficient sampling strategies can be used (such as acceptance/rejection sampling for example (see [START_REF] Pareschi | An introduction to Monte Carlo methods for the Boltzmann equation[END_REF] for instance)). Finally, for the next time step the number of particles becomes

N n+1 i = e -∆t/ε 2 Ñ n i + M n i
and the above procedure can be iterated to advance in time.

The method described deserves some remarks. First, it permits to overcome the stiffness related to the characteristic speeds allowing for time steps independent of the scaling parameter ε. Second, it automatically diminishes the computational complexity as ε → 0: the number of particles diminishes exponentially with the scaling parameter. This is opposite to standard Monte Carlo or more general particle approaches for kinetic equations, where the computational complexity remains unchanged with ε. This highlights the asymptotically complexity diminishing (ACD) property of our strategy.

A finite volume method for the macroscopic equation

In this subsection, we focus on the space-time discretization of the macroscopic equation (3.4). One natural choice would be the following time discretization (see [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF])

ρ n+1 -ρ n ∆t + 1 ε ∇ x • vg n+1 = 0, (3.15)
where the value g n+1 would be recovered from the Monte Carlo solution of equation (3.5). This would permit to overcome the stiffness in (3.15) and would lead to a numerical scheme which is stable for all time steps and is consistent with the limit diffusion equation scheme. However, due to the stochastic approximation, g n+1 may suffer from noise which will be passed on the macroscopic unknown ρ n+1 . This statistical noise, inherent to particle methods, will be amplified by the factor 1/ε so that from one side it would be more difficult to stabilize a finite volume method and from the other side the asymptotic behavior could be lost. To overcome this drawback, we then propose to separate in the expression of g n+1 given by (3.14) the macroscopic term containing ρ n which gives rise to the limit diffusion equation (2.6) from the microscopic ones. Thus, starting from (3.14), let us first compute v • ∇ x g n+1 :

v • ∇ x g n+1 = e -∆t/ε 2 v • ∇ x gn -ε(1 -e -∆t/ε 2 ) ( v • ∇ x (v • ∇ x ρ n M ) -v • ∇ x ( v • ∇ x gn M ) ) .
Since v • ∇ x g n does not depend on v and vM = 0, we get

v • ∇ x g n+1 = e -∆t/ε 2 v • ∇ x gn -ε(1 -e -∆t/ε 2 ) v • ∇ x (v • ∇ x ρ n M )
and finally

v • ∇ x g n+1 = e -∆t/ε 2 v • ∇ x gn -ε(1 -e -∆t/ε 2 )∆ x ρ n . (3.16)
Plugging expression (3.16) into the macro equation (3.15), we get

ρ n+1 -ρ n ∆t + 1 ε e -∆t/ε 2 v • ∇ x gn -(1 -e -∆t/ε 2 )∆ x ρ n = 0.
We now observe that the above time discretization permits to remove the stiffness that appeared in (3.15) when ε → 0, but the diffusion turns out to be explicit. Therefore, this requires a parabolic CFL condition of type ∆t ≤ C∆x 2 . To avoid such restrictive condition, an implicit treatment of the asymptotic heat equation may be adopted so that the time discretization of the macro unknown becomes

ρ n+1 -ρ n ∆t + 1 ε e -∆t/ε 2 ∇ x • vg n -(1 -e -∆t/ε 2 )∆ x ρ n+1 = 0. (3.17)
Note that this choice allows larger time steps, and that equation (3.14) may be modified accordingly in this case:

g n+1 = e -∆t/ε 2 gn + (1 -e -∆t/ε 2 )ε -v • ∇ x ρ n+1 M + ∇ x • vg n M . (3.18)
Luckily, equations (3.17) and (3.18) can still be solved separately by first solving implicitly equation (3.17) which needs only the knowledge of gn which is obtained by the advection step of the time splitting algorithm described in (3.2). Successively one can solve (3.18) with the value ρ n+1 found in (3.17). We finally discuss the space discretization. We choose for (3.17) a simple centered finite difference approximation of ∇ x • vg n . Concerning the space approximation of the Laplacian term, we use an ADI (Alternating Direction Implicit) method (see [START_REF] Peaceman | The numerical solution of parabolic and elliptic differential equations[END_REF]) to avoid the inversion of large matrices in the multidimensional cases, still avoiding parabolic CFL condition. We give details in the two-dimensional case d x = 2 in space (x = (x, y)). In this case, the method reduces to the following splitting, starting from ρ n . 1. We solve over a time step ∆t the equation

∂ t ρ + 1 2ε e -∆t/ε 2 v • ∇ x gn -(1 -e -∆t/ε 2 )∂ xx ρ = 0, using a Crank-Nicolson time discretization to get ρ ⋆ .
2. Starting from ρ ⋆ , we solve over a time step ∆t the equation

∂ t ρ + 1 2ε e -∆t/ε 2 v • ∇ x gn -(1 - e -∆t/ε 2 )∂ yy ρ = 0, using a Crank-Nicolson time discretization to get ρ n+1 .
Given x i = x min + (i -1)∆x, y j = y min + (j -1)∆y, a cartesian mesh of size (∆x, ∆y) with N x and N y points, we define ρ n i,j an approximation of ρ(t n , x i , y j ), and vg n i,j an approximation of the first moment vg(t n , x i , y j ) , with v = (v x , v y ). Then, the corresponding time-space discretization is the following.

1. For each fixed j such that 1 ≤ j ≤ N y , solve the one-dimensional system of size N

x - ∆t 2∆x 2 (1 -e -∆t/ε 2 )ρ ⋆ i+1,j + 1 + ∆t ∆x 2 (1 -e -∆t/ε 2 ) ρ ⋆ i,j - ∆t 2∆x 2 (1 -e -∆t/ε 2 )ρ ⋆ i-1,j = ρ n i,j - ∆t 2ε e -∆t/ε 2 v x gn i+1,j -v x gn i-1,j 2∆x + v y gn i,j+1 -v y gn i,j-1 2∆y + ∆t 2 (1 -e -∆t/ε 2 ) ρ n i+1,j -2ρ n i,j + ρ n i-1,j ∆x 2 ,
which furnishes the intermediate values ρ ⋆ 1≤i≤Nx,j .

2. For each fixed i such that 1 ≤ i ≤ N x , solve the one-dimensional system of size N y

- ∆t 2∆y 2 (1 -e -∆t/ε 2 )ρ n+1 i,j+1 + 1 + ∆t ∆y 2 (1 -e -∆t/ε 2 ) ρ n+1 i,j - ∆t 2∆y 2 (1 -e -∆t/ε 2 )ρ n+1 i,j-1 = ρ ⋆ i,j - ∆t 2ε e -∆t/ε 2 v x gn i+1,j -v x gn i-1,j 2∆x + v y gn i,j+1 -v y gn i,j-1 2∆y + ∆t 2 (1 -e -∆t/ε 2 ) ρ ⋆ i,j+1 -2ρ ⋆ ij + ρ ⋆ i,j-1 ∆y 2 ,
which gives the value of the density at time n + 1:

ρ n+1 i,1≤j≤Ny .
It is known that the natural extension to the three dimensional case of the above ADI method is conditionally stable (i.e. ∆t = O ∆x 2 ). Consequently the resulting micro-macro method is also conditionally stable in the full three dimensional case. However, some unconditionally stable variants of the ADI method have been derived, for instance in [START_REF] Sharma | A fast semi-implicit method for anisotropic diffusion[END_REF]. Alternatively, a fully implicit scheme (which involves the resolution of a linear system in ρ) can be used. In both situations, the resulting full micro-macro method can be modified to become unconditionally stable.

We conclude this section by stressing that the proposed method is globally unconditionally stable, independently of ε, thanks to the two following ingredients: (i) the approximation of the microscopic part is unconditionally stable due to the reformulation and the use of a particle method; (ii) the approximation of the macroscopic part is free from the CFL parabolic condition due to the use of a implicit or semi-implicit method for the diffusion term.

Numerical results

We present in this section several numerical results with the scope of understanding the behaviors and the characteristics of the proposed method. The two first subsections aim at validating the Hybrid Micro-Macro Monte Carlo scheme described in the previous section. They concern four dimensional phase space experiments (d x = d v = 2). In these cases, our method is compared to other known approaches: namely a micro-macro grid method (see [START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF]) considered as a reference, a full particle scheme performed on the original distribution function f , or a finite difference method for the asymptotic diffusion equation. Diagnostics on macroscopic variables as well as v-dependent quantities are presented. Numerical cost reduction measures -linked to the ACD property -are discussed, in order to point out the efficiency of our method. After these validating tests, we highlight other advantages of our method, namely the possibility of considering space-dependent scaling parameter ε and the performances in higher dimension. To that aim, in the third subsection, we present a slightly different model, in which a space-dependent collisional frequency is considered. We show how the Hybrid Micro-Macro Monte Carlo method is automatically able to adapt the number of particles to the collision frequency. Finally, we propose a full 3D in space and 3D in velocity test case to show the ability of our approach to deal with problems of physical relevance.

Two-dimensional test cases with constant scaling parameter and equilibrium initial condition

This subsection is devoted to 2D-2D numerical results in the constant ε case and with an equilibrium initial condition given by

f (t = 0, x, v) = ρ(t = 0, x)M (v), x ∈ [0, 4π] 2 , v ∈ R 2 (4.19) with ρ(t = 0, x) = 1+ 1 2 cos x 2 cos y 2 , M (v) = 1 2π exp -|v| 2 2
, so that g(t = 0, x, v) = 0. Periodic boundary conditions are imposed in space. The aim is to compare our Hybrid Micro-Macro Monte Carlo approach, referred to as "MM-MC" in the sequel, with other methods. Namely, we compare the MM-MC with the following methods.

• A micro-macro grid approach (referred to as "MM-G"). Considering the model described by equations (3.4)-(3.5), g is discretized on a grid in space and in velocity. This method has been shown to be uniformly stable independently of the scaling parameter ε and asymptotically consistent with the diffusion equation (see [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF] for details). When the meshes are refined enough, we will take the solution provided by MM-G as a reference solution for comparisons. • A full particle method (referred to as "Full MC"). We use a classical particle method directly on model (1.1). We refer to [START_REF] Birsdall | Plasma Physics Via Computer Simulation[END_REF] for details about this approach. At the limit ε → 0, this method suffers from a severe constraint on the time step: it has to be of order ε 2 . For this reason, we shall only use it in regimes where ε is of order 1.

• A finite difference discretization for the limit equation (referred to as "Limit"). When the scaling parameter ε is small, we use finite differences in space and an ADI Crank-Nicolson scheme in time for discretizing the limit diffusion equation, as explained in Subsection 3.3.

The solution of a refined MM-G will be taken as a reference solution.

Limit and intermediate regimes. First, we check the asymptotic behavior of our MM-MC method in the limit (ε = 0.0001) and intermediate (ε = 0.1) regimes and verify that the right limit is reached. The deviation from equilibrium for the different values of the scaling parameter can be measured by performing a diagnostic on the non equilibrium part of the momentum variable v x g (t, x). To that aim, in Figure 1, we report the momentum at the final time of the simulation v x g (T = 2, x). This quantity is fixed to zero at the beginning, in fact equilibrium initial data are chosen, and it should remain zero in the limit ε → 0. Table 1 reports the numerical parameters chosen for the simulations.

For ε = 0.0001 (respectively ε = 0.1), results given by the MM-MC approach are presented in Figure 1a (respectively 1b) while results given by the MM-G approach are presented in Figure 1c (respectively 1d). We observe from the simulations that the component of the momentum given by the perturbation g diminishes as ε goes to zero, as theoretically expected. So we can claim that both approaches share a good asymptotic behavior, i.e. are asymptotically consistent. Moreover, we can observe that the MM-MC method produces some statistical oscillations which automatically diminish when the scaling parameter decreases, since the number of particles used in the simulation is proportional to the perturbation g (let us remark that when ε = 10 -6 there is almost no particles at T = 2). We emphasize that this noise can also be reduced by imposing a higher number of particles, or equivalently a smaller typical weight m p for each particle.

We finally compare the results by considering the slice y = 0. Figure 2a shows the slice density ρ(T = 2, x, y = 0) while in Figure 2b the slice momentum v x g (T = 2, x, y = 0) for ε = 0.1 and ε = 0.0001 is shown. We can finally claim that for these test cases, the MM-MC method is able to reproduce the correct solution identified by the reference furnished by the MM-G scheme.

Concerning the choice of the time step, the MM-MC method is stable independently of the time step ∆t. On the contrary, the MM-G scheme, as other grid based schemes, suffers from a constraint on ∆t related to the speed of the particles. We also stress that these results can not be obtained by a standard Monte Carlo approach since the stiffness of the equation will cause too severe time step limitations while statistical noise will remain at the contrary to our approach in which this error disappears in the asymptotic limit.

Kinetic regime. In this part, we are interested in the numerical results when ε is of order 1, we refer to this regime as the kinetic regime. The initialization is the same as for the previous case (see (4.19)). Table 2 reports the numerical parameters chosen for the following simulations. Figure 3 (respectively Figure 4) shows the density profile (respectively the error on the density profile) at time T = 2, whereas the momentum v x f is reported in Figure 5. For ε = 0.5 (respectively ε = 1), we report the density obtained when a full Monte Carlo method (referred to as "Full MC") is employed to compute the solution. In Figure 3a (respectively 3d) these density profiles are shown for the two different values of ε. As well known, a full Monte Carlo method suffers from large numerical noise especially when non stationary problems are solved. This fact is made clear by the simulation reported in which 10 8 particles are used in the simulations. This corresponds to a particle weight of order 10 -6 . On the contrary with respect to the full particle scheme, the MM-MC produces solutions with much less oscillations as shown in Figure 3b for ε = 0.5 (respectively 3e for ε = 1) with the same typical weight m p . In both cases, the results of the MM-MC method are very close to the ones of the MM-G reference method shown in Figure 3c for ε = 0.5 (respectively 3f for ε = 1). These observations are made more clear thanks to Figure 4 where the difference between the density obtained by the reference solution and the density obtained by Full MC on the the left (Figures 4a and4c) and the difference between the reference solution and the density obtained by MM-MC on the right (Figures 4b and4d) are displayed for ε = 0.5 and ε = 1. One can see that the error produced by Full MC is one order of magnitude larger than the error produced by MM-MC.

Let us also remark that the reference solution and the Monte Carlo methods employ different time steps as shown in Table 2. This causes different numerical diffusions in the solutions responsible for an increase of the difference between the reference and the Monte Carlo solutions.

Concerning the momentum, we show the quantity v x f given by the Full MC method in Figure 5a for ε = 0.5 and in Figure 5d for ε = 1. The numerical noise strongly affects this quantity. For the MM-MC method, the corresponding quantity v x g is presented in Figure 5b and in Figure 5e for the two different values of the scaling parameters while the results for the MM-G approach are presented in Figure 5c and in Figure 5f. As for the case of the limit and intermediate regimes, we report the density and the momentum profiles for a slice in y = 0. The density ρ(T = 2, x, y = 0) is shown in Figure 6a, while the momentum v x f (T = 2, x, y = 0) is shown in Figure 6b for the case ε = 1. Let us remark that the density slice of MM-MC is slightly more damped than the one of MM-G. This has also been observed in the error plots in Figure 3e. In fact, as already stated, this behavior is due to the time discretization and in particular to the different time steps employed. Indeed, decreasing the time step for MM-MC enables to decrease the numerical diffusion.

Asymptotically complexity diminishing property. In this part, we are interested in measuring the computational cost of the method as well as measuring the so-called ACD property of our MM-MC method. This latter means the capacity of the scheme to diminish its computational complexity as the limit diffusion equation is approached. We highlight these properties through the evolution in time of the total number of particles N p (t). In Figure 7a, we present the results obtained in the 2D-2D case where the initial condition is given by (4.19) and where ε takes the values 1, 0.5, 0.2 and 0.1. As the perturbation g is initially zero, the number of particles at the beginning increases due to the fact that the transport phase induces departure from equilibrium. Then, due to collisions, the system approaches the equilibrium again and the number of particles decreases. The smaller is ε, the smaller is g and the smaller is the number of particles we need to represent this component of the solution. This behavior is confirmed by Figure 7b, where results are given for ε = 0.01, ε = 0.001 and ε = 0.0001. However, it is important to stress that the global cost does not only depend on the number of particles since the sampling procedure has a cost which is not negligible compared to the rest of the scheme. More in details, there are two additional steps in the scheme compared to a standard Monte Carlo approach which increase the cost. These are:

• sampling from the function P n,± (x, v),

• solve the macroscopic equation.

Even if standard numerical techniques are employed in the whole algorithm (LU solver for the ADI step, classical techniques for the particles part), it is never simple to estimate the efficiency of the current implementation. More generally, the computational cost of a numerical method may depend on different sources, and instead of doing a direct measure of the computational cost, we try in the following to give a measure of the complexity of our method in comparison with the complexity of a standard Monte Carlo approach. In order to do that, we start from giving the cost needed to compute the solution of the kinetic equation in the diffusive scaling by using a standard Monte Carlo method. This can be estimated by the simple following formula

C M C (f ) = C f N M
where C f is a suitable constant, N is the number of particles and M is the number of time steps. The time step of standard MC approach is of the order of ε∆x/ max(v) due to the type of scaling. On the other hand, the cost related to the computation of a solution with the MM-MC method is

C M M -M C (f ) = C g N g M g + C ρ N dx x M g + C g N dx x N dv v M g
where C g is a constant very close to C g ≃ C f since the transport parts of the two Monte Carlo methods are equivalent, N g is the number of particles used to describe the perturbation, N x (resp. N v ) the number of grid points in physical (resp. velocity) space and d v and d x respectively the dimensions in velocity and space. Now, since the cost related to the solution of the macroscopic equation C ρ N dx x M g is in general much lower than the cost of any type of solver applied to kinetic equations, i.e. C ρ ≪ C g , then, at least as a first approximation, the computational overhead can be considered negligible for this part of the method. Concerning M g , this is the number of time steps used in our method which is always M g ≤ M . In particular, the ratio between the number of time steps needed to reach the final solution for the two methods is such that

M g = εM
which means that we experiment very large gains when close to the asymptotic limit and smaller ones when close to the kinetic regime. Concerning the ratio between the number of particles of the standard Monte Carlo N and the number of particles employed in the solution of the perturbation g, i.e. N g /N , it is always in favor of the MM-MC method here presented. In particular, it goes to zero as ε → 0 and the numerical tests show that it decreases exponentially in the space homogeneous case and faster than linearly in the space non homogeneous case with ε. Finally, we consider the quantity C g N dx

x N dv v M g . This cost is due to the sampling of new particles from the distribution P n,± (x, v). This sampling should be done at each time step in each spatial cell and it depends on a grid in velocity space which should be adopted to successively perform the sampling. This part of the numerical method is comparable to a standard grid method for solving the kinetic equation. However, it still depends on the scaling through the term (1e -dt/ε 2 ) which means that C g depends on ε and it goes to zero as ε → 0. Thus, to resume, the cost of the MM-MC method, as it has been presented, is comparable to the cost of a grid method when ε is of order one and comparable to a cost of a macroscopic scheme when ε → 0, with a cost which decreases with the scaling parameter. Let us also remark that the computational performances of the MM-MC strategy shown here can be additionally improved by avoiding the construction of a grid in velocity space with the use of acceptance/rejection techniques.

Two-dimensional test cases with non-equilibrium initial condition and constant scaling parameter

In this part, we consider the following non-equilibrium initial data for the distribution function f with ρ(t = 0, x) = 1 + 1 2 cos x 2 cos y 2 . The tests presented use the parameters detailed in Table 3. In this part, we look at quantities depending only on the velocity space: i.e. v x and v y . More precisely, we compute f (t, x, v)dx and g(t, x, v)dx from the MM-MC method and we compare them to the results from the MM-G one. We present in Figures 8 the quantity f (t, x, v)dx at different times for ε = 1. In details, Figures 8a, 8b and 8c (resp. 8d, 8e and 8f) report the results for the MM-MC method (resp. the MM-G method) at times T = 0, T = 0.5 and T = 2. When ε diminishes equilibrium is reached faster, results are not reported. Finally, in Figure 9, we compare the integral of the perturbation g(T, x, v)dx at T = 0.2 for different values of the scaling parameter ε. Results given by the MM-MC method, respectively the MM-G method, are presented in Figure 9a, respectively in Figure 9d, for ε = 1, in Figure 9b, respectively in Figure 9e, for ε = 0.5 and in Figure 9c, respectively in Figure 9f, for ε = 0.1. We summarize by observing from these test cases that the MM-MC method furnishes results which are in a very good agreement with those obtained by the reference MM-G scheme. From a computational point of view, in this situation the initial number of particles is as large as the one used for a standard Full MC method. This is due to the fact that the initial number of particles has to be large enough to represent g(T = 0, x, v) = 0 with reasonable accuracy. Thus, the first iterations in time require the same computational cost as a classical full particle method. Then collisions drive the system to its equilibrium and the number of particles diminishes. We finally report the evolution in time of the number of particles for different values of ε in Figure 10: ε = 0.1, 0.2, 0.5, 1. The numerical parameters are: N x = N y = 128, m p = 5 × 10 -6 and ∆t = 0.1.

f (t = 0, x, v) = 1 4π exp - |v -2| 2 2 + exp - |v + 2| 2 2 ρ(t = 0, x), x ∈ [0, 4π] 2 , v ∈ R 2 , (4.20

Two-dimensional test cases with non constant scaling parameter

In this part, we focus on a more realistic multi-scale kinetic problem where a space dependent function a(x) is introduced to enrich the model. Our aim is to study the behavior of the previously described scheme in such a case. We expect that the method is able to automatically adapt the number of particles in each cell, according to the values of a(x) ≥ 0. In details, we consider the following equation

∂ t f + 1 ε v • ∇ x f = a 2 (x) ε 2 (ρM -f ), (4.21)
where (x, v) ∈ [0, 4π] 2 × IR 2 and t ≥ 0. Periodic boundary conditions are imposed in space. We choose for a(x) a step function whose values are 10 or 0.5 (see Figure 11). We recall that the asymptotic diffusion model in this context becomes

∂ t ρ -∇ x • 1 a 2 (x) ∇ x ρ = 0.
We now briefly present the main steps of our method, which is slightly modified to take into account this new context • micro part: advance the particles

x n+1 k = x n k + ε a 2 (x) (1 -e -a 2 (x)∆t/ε 2 )v n k ,
• micro part: solve the modified interaction part by sampling and discarding particles • macro part: advance the macroscopic part

g n+1 = e -a 2 (x)∆t/ε 2 gn + ε a 2 (x) (1 -e -a 2 (x)∆t/ε 2 ) -v • ∇ x ρ n M + ∇ x • vg n M , 0 2
ρ n+1 -ρ n ∆t + 1 ε ∇ x • ve -a 2 (x)∆t/ε 2 gn -∇ x • (1 -e -a 2 (x)∆t/ε 2 ) a 2 (x) ∇ x ρ n+1 = 0.
Obviously, when a(x) = 1, we recover the scheme presented in Section 3.

Our Micro-Macro Monte Carlo (MM-MC) method is compared, as previously, to a grid based reference solution referred to as "MM-G" where the function a(x) has been included. The results obtained by MM-G use the following numerical parameters: N x = N y = 64 for the spatial grid, N vx = N vy = 32 for the velocity grid and ∆t = 10 -4 . For our MM-MC solver, the numerical parameters are: N x = N y = 64, N p = 5000 (hence m p ≈ 7.7 × 10 -6 ), ∆t = 0.1. The initial condition is chosen as

f 0 (x, v) = 1 4π exp - |v -u| 2 2 + exp - |v + u| 2 2 ρ(0, x),
with u = (2, 2) and ρ(0, x) = 1 + 1 2 cos( x 2 ) cos( y 2 ). In Figure 12, we first plot the space dependent density ρ(T = 0.5, x) at the final time T = 0.5 obtained by the two different solvers (with ε = 0.5). We can see that the two solutions are very close, which confirms the good behavior of the MM-MC method in this context. Keeping ε = 0.5, we then plot the velocity repartition of the microscopic part g at position x = (1, 1), which corresponds to a small value of the a function (a(1, 1) = 0.5) and then the point x = (1, 1) belongs to a microscopic region. Indeed, we can observe in Figure 13 that the distribution function has not reached the equilibrium state whereas the equilibrium state is reached for x = (2π, 2π) (which corresponds to a(2π, 2π) = 10), which belongs to a macroscopic region. Next, in Figure 14, we show the spatial distribution of the particles used to sample the perturbation g, for different values of ε (ε = 0.1, 10 -3 , 10 -4 ). From this figure, it is clear that particles are automatically created in the zones in which the departure from the equilibrium state is large and are discarded otherwise. Typically, particles leave the central macroscopic zone for which the function a(x) is quite large whereas the zone in which a(x) is small still contains particles even when ε = 10 -3 . Obviously, when ε is small enough (ε = 10 -4 ), the non homogeneity of a(x) has essentially no influence on the microscopic part which is then sampled using very few particles. Finally, in Figure 15, the time history (semi-logscale) of the number of particles is displayed. Initially, the number of particles is chosen as 5000 per cell (hence the mass of a particle is close to m p ≈ 7.7×10 -6 ). The corresponding number of particles in the whole domain is around 3.2 × 10 7 . After one iteration, we can see that this number is decreased by a factor of order 1/ε and then a linear decrease is observed (which corresponds to an exponential dependence of the number of particles with respect to time). Let us notice that the number of particles employed is a direct consequence of the magnitude of the perturbation. On the contrary a standard Monte Carlo approach would require a constant number of particles for all the simulation, which means around 15 times more particles than the one used by MM-MC method.

Three-dimensional test cases

We conclude this section with a full three dimensional in space and three dimensional in velocity simulation. We consider the following initial condition

f 0 (x, v) = 1 2 (2π) 3/2 exp - |v -u| 2 2 + exp - |v + u| 2 2 ρ(0, x), (4.22)
with u = (2, 2, 2), ρ(0, x) = 1 + 1 2 cos( x 2 ) cos( y 2 ) cos( z 2 ) and x = (x, y, z) ∈ [0, 4π] 3 , v = (v x , v y , v z ) ∈ R 3 . Periodic boundary conditions are considered in space. In this framework, we only consider our MM-MC solver, since both grid and full particle approaches are too expensive in this case. The numerical parameters of the MM-MC method are the following: N x = N y = N z = 16, a typical weight of order 6 × 10 -5 , ∆t = 0.1. We present the evolution in time of the quantity f (t, x, v)dx, obtained by our MM-MC scheme. Results obtained for ε = 1 are plotted in Figure 16, for T = 0, T = 0.2 and T = 1. In addition to the presented figures, we have verified that the equilibrium is reached in shorter time when ε is reduced. We finally report in Figure 17 the evolution of the number of particles in this full 3D-3D case with equilibrium initial condition f (t = 0, x, v) = ρ(t = 0, x)M (v), x ∈ [0, 4π] 3 , v ∈ R 3 , with ρ(t = 0, x) = 1 + 1 2 cos x 2 cos y 2 cos z 2 , M (v) = 1 (2π) 3/2 exp -|v| 2

2

. The space discretization is: N x = N y = N z = 16 with a typical weight for the particles of order 5 × 10 -4 and ∆t = 0.1.

From the results, we observe the same behaviors as for the 2D-2D case both for high values of the scaling parameter ε (Figure 17a) and in the limit regime (Figure 17b). The number of particles really employed in the construction of the solution is a function of the scaling parameter and it initially increases since the motion creates a departure from the equilibrium and successively it diminishes. This is a sign that a global equilibrium is about to be reached. In the last Figure 18, we show the evolution of the number of particles for a non-equilibrium initial condition (4.22). In this last case, the numerical parameters are N x = N y = N z = 16, a typical weight for the particles of order 5 × 10 -4 and ∆t = 0.1. The same behaviors as for the 2D-2D case are observed: the number of particles used to compute the solution decreases as a function of time and as a function of the scaling parameter ε.

Conclusion

In this work, we have presented a new numerical method for kinetic equations in the diffusion scaling. The main lines of the proposed method are the following. (i) It is based on a suitable merging between Monte Carlo approach and a finite volume method. (ii) The Monte Carlo part of the solution is constructed to solve the perturbation part of the model while the finite volume part is designed to solve the macroscopic equilibrium part of the solution. The method enjoys some nice properties which we summarize as follows. (a) Compared to a standard Monte Carlo approach for solving the kinetic equation its statistical noise is smaller and it diminishes when the scaling parameter ε decreases. (b) It is uniformly stable with respect to the scaling parameter as well as with respect to the space mesh size. (c) Its computational cost as well as its variance diminish as the equilibrium is approached. (d) It does not need artificial transitions to pass from the microscopic description to the macroscopic one. The numerical results illustrate the efficiency of the proposed method compared to the standard Monte Carlo approaches as well as compared to standard deterministic methods based on a grid in the phase space.

In the future, we intend to work more extensively on the full three dimensional case and to consider problems with more physical relevence. In this respect, the extension of the present approach to the case of the Boltzmann operator is also under study.

1 Figure 1 :

 11 Figure1: The component of the momentum given by the perturbation g, v x g , is shown for the MM-MC and the MM-G methods and for different values of the scaling parameter ε in the limit and intermediate regimes.

Figure 2 :

 2 Figure2: Slices of the density ρ(T = 2, x, y = 0) and of the momentum v x g (T = 2, x, y = 0) in the limit and intermediate regimes. Comparisons between the MM-MC, the MM-G and the limit diffusion equation.

Figure

  Figure Scheme N x , N y Typical weight N vx , N vy ∆t ε T 3a, 5a Full MC 128 10 -6 0.1 0.5 2 3d, 5d Full MC 128 10 -6 0.1 1 2 3b, 5b MM-MC 128 10 -6 0.1 0.5 2 3e, 5e MM-MC 128 10 -6 0.1 1 2 3c, 5c MM-G 128 10 0.001 0.5 2 3f, 5f MM-G 128 10 0.005 1 2

  MM-G method, T = 2, ε = 1

Figure 3 :

 3 Figure 3: The density profile ρ is shown for the MM-MC, the MM-G and the Full MC methods in the kinetic regime.

1 Figure 4 :

 14 Figure 4: Plots of the difference between the density between Full MC and MM-G (figures (a) and (c)) and between MM-MC and MM-G (figures (b) and (d)) for ε = 0.5 and ε = 1.

1 Figure 5 :Figure 6 :

 156 Figure 5: The component of the momentum v x g (resp. v x f ), given by the perturbation is shown for the MM-MC, the MM-G and the Full MC methods, for different values of the scaling parameter ε in the kinetic regime.

)

Figure 7 :

 7 Figure 7: Total number of particles as a function of time for a 2D-2D case with equilibrium initial condition.

Figure Scheme N x , N y 10

 10 Figure Scheme N x , N y Typical weight N vx , N vy ∆t ε T 8a, 8b,8c, 9a MM-MC 128 5 × 10 -6 0.1 1 various 9b MM-MC 128 5 × 10 -6 0.1 0.5 0.2 9c MM-MC 128 5 × 10 -6 0.1 0.1 0.2 8d, 8e, 8f, 9d MM-G 128 20 0.01 1 various 9e MM-G 128 20 0.001 0.5 0.2 9f MM-G 128 20 0.0005 0.1 0.2

1 - 6 1 - 6 1 - 6 1 - 6 1 - 6

 1616161616 MM-MC, T = 0, ε = MM-MC, T = 0.5, ε = MM-MC, T = 2, ε = dy (d) MM-G, T = 0, ε = dy (e) MM-G, T = 0.5, ε = dy (f) MM-G, T = 2, ε = 1

Figure 8 :

 8 Figure 8: Integral of the distribution function in space f (T, x, v)dx at times T = 0, T = 0.5 and T = 2 for ε = 1.

1 - 6 5 - 6 1 - 6 1 - 6 5 - 6 1 Figure 9 :

 165616165619 Figure 9: Integral of the perturbation in space g(T, x, v)dx at time T = 0.2 for ε = 1, ε = 0.5 and ε = 0.1.

Figure 10 :

 10 Figure 10: Total number of particles as a function of time for a non-equilibrium initial condition and different values of the scaling parameter.

Figure 11 :

 11 Figure 11: Profile of the function a(x).

Figure 13 :

 13 Figure 13: Velocity dependence of the distribution function at position (x, y) = (1, 1). Left: MM-MC method; right: MM-G method.

Figure 14 :Figure 15 :

 1415 Figure 14: Spatial repartition of the particles at time T = 0.5. Left: ε = 0.1; middle: ε = 10 -3 ; right: ε = 10 -4 .

1 Figure 16 :

 116 Figure 16: The integral of the distribution function in space x f (T, x, v)dx for ε = 1 and different times. Full three dimensional case.

Figure 17 :

 17 Figure 17: Total number of particles as a function of time for the full 3D-3D case with equilibrium initial condition.

Figure 18 :

 18 Figure 18: Total number of particles as a function of time for the full 3D-3D case with a nonequilibrium initial condition.

Table 1 :

 1 Numerical parameters for the MM-MC and the MM-G methods.

			N x , N y m p	N vx , N vy ∆t	ε	T
	1a	MM-MC 128	10 -6		0.1	0.0001 2
	1b	MM-MC 128	10 -6		0.1	0.1	2
	1c	MM-G	128		10	0.01	0.0001 2
	1d	MM-G	128		10	0.0005 0.1	2

Table 2 :

 2 Scheme N x , N y Typical weight N vx , N vy ∆t Numerical parameters for the Full MC, the MM-MC and the MM-G methods for the kinetic regime.

	ε	T

Table 3 :

 3 Numerical parameters for the MM-MC and the MM-G methods for non equilibrium initial data.

Acknowledgements

N. Crouseilles and M. Lemou are supported by the French ANR project MOONRISE ANR-14-CE23-0007-01 and by the Enabling Research EUROFusion project CfP-WP14-ER-01/IPP-03. A. Crestetto is supported by the French ANR projects MoHyCon ANR-17-CE40-0027-01 and ACHYLLES ANR-14-CE25-0001.