Abstract
We develop a new Monte Carlo method that solves hyperbolic transport equations with stiff terms, characterized by a (small) scaling parameter. In particular, we focus on systems which lead to a reduced problem of parabolic type in the limit when the scaling parameter tends to zero. Classical Monte Carlo methods suffer of severe time step limitations in these situations, due to the fact that the characteristic speeds go to infinity in the diffusion limit. This makes the problem a real challenge, since the scaling parameter may differ by several orders of magnitude in the domain. To circumvent these time step limitations, we construct a new, asymptoticpreserving Monte Carlo method that is stable independently of the scaling parameter and degenerates to a standard probabilistic approach for solving the limiting equation in the diffusion limit. The method uses an implicit time discretization to formulate a modified equation in which the characteristic speeds do not grow indefinitely when the scaling factor tends to zero. The resulting modified equation can readily be discretized by a Monte Carlo scheme, in which the particles combine a finite propagation speed with a timestep dependent diffusion term. We show the performance of the method by comparing it with standard (deterministic) approaches in the literature.
AsymptoticPreserving Monte Carlo methods for transport equations in the diffusive limit
G. Dimarco^{1}, L. Pareschi^{1} and G. Samaey^{2}
^{1}Department of Mathematics and Computer Science, University of Ferrara, Italy
^{2}NUMA (Numerical Analysis and Applied Mathematics), Department of Computer Science, KU Leuven, Celestijnenlaan 200A, 3001 Leuven, Belgium
Key words.
Transport equations, diffusion limit, Monte Carlo methods, asymptoticpreser–
ving schemes.
AMS subject classifications. 65C05, 35B25, 65M75, 76R50, 82C70
1 Introduction
In many situations in which hyperbolic transport equations intervene, such as neutron transport [14], radiative transfer [5, 10, 11], plasma physics [6] or semiconductor simulation [42], the multiscale nature of the phenomena involved causes large difficulties for the development of efficient numerical methods. In fact, the scaling parameters that characterize the relevant time scales which determine evolution of such problems may differ by several order of magnitude, making the problem very stiff [14]. This stiffness limits the maximal timestep that can be taken when using an explicit discretization. In addition, in the limit when the scaling parameter tends to zero, the equations that are used to model these phenomena may change character, passing from a hyperbolic to a parabolic structure [3]. In the particular setting of a diffusive scaling [3, 13] that is considered in this work, the characteristic speeds of the hyperbolic system grow to infinity as the scaling parameter tends to zero, causing severe CFL restrictions in standard explicit numerical methods, both in deterministic (gridbased) and stochastic Monte Carlo (particlebased) approaches.
Thus, the first idea to deal with such systems consists in using implicit time integration methods. Unfortunately, these techniques are only possible to implement for gridbased space discretizations. Moreover, the main drawback of such approaches is that they often lead to large nonlinear systems very hard to invert due to the high dimension of the equations involved. For this reason, it is desirable to develop numerical methods that are not fully implicit, yet able to overcome the computational cost caused by the stiffness by avoiding time step limitations related to the scaling parameter. A particularly appealing class of schemes, developed for facing these situations, are the socalled asymptoticpreserving schemes [4, 7, 9, 15, 23, 27, 28, 29, 30, 32, 33, 37, 39, 40], which degenerate to consistent discretization of the limiting problem when the scaling parameter is set to zero. Recently, this approach has been considered in the framework of ImplicitExplicit (IMEX) RungeKutta schemes with the aim of deriving high order numerical methods that are accurate in all regimes, i.e., regardless of the value of the scaling parameter [23, 7, 8, 22]. In particular, these schemes also preserve the desired order of accuracy, even in the limit when the scaling parameter tends to zero. Recently, via projective integration [25], also fully explicit methods for stiff hyperbolic transport equations have been developed [36, 35, 34].
In this work, we concentrate on the development of asymptoticpreserving Monte Carlo methods for solving hyperbolic transport equations in the diffusion limit. The Monte Carlo approach represents a very popular method to deal with transport equations due to its flexibility and low computational cost [12, 5, 41, 46, 38]. However, one of the main limitation of Monte Carlo methods is the difficulty to construct schemes which work uniformly for all values of the scaling parameters [12, 46]. For this reason, various modified Monte Carlo techniques have been recently proposed in the case of the socalled hydrodynamic scaling [20, 21, 17, 44, 46, 47]. To the best of our knowledge, up to now, a Monte Carlo scheme which is able to work uniformly in the diffusive scaling without causing severe time step limitations has not been constructed. The main motivation for not using particle methods in this framework is that the characteristic speeds grow to infinity in the diffusive limit. Thus, the dominant technique to overcome this problem consists nowadays in employing domain decomposition strategies, in which a Monte Carlo discretization of the hyperbolic transport equation is solved in regions in which the scaling parameter is large, and the limiting equation is solved (deterministically) in regions where the scaling parameter is small. These approaches have been largely studied for kinetic equations both for the diffusive [2, 16, 32] and for the hydrodynamic scaling [18, 19]; see also, e.g., [24] for related ideas. However, even if these methods are very efficient, they are affected by some difficulties due to the fact that it is not always a simple task to define the different regions of the domain in which the use of a macroscopic model is fully justified.
In this paper, our main goal is to develop an asymptoticpreserving Monte Carlo method in the diffusive scaling. The method does not rely on domain decomposition strategies, nor on coupling with a deterministic discretization of the limiting equation. Instead, to overcome the parabolic stiffness, the main ingredients are a suitable reformulation of the original system based on an implicit time discretization [7, 8], which leads to a modified equation where the characteristic speeds are bounded in terms of the scaling parameters. An appropriate splitting strategy for this equation then permits the construction of a Monte Carlo scheme that works independently of the value of the scaling parameter and that automatically degenerates to a classical random walk method for limiting the diffusion equation. The method is first constructed by using the Goldstein–Taylor model [26] under the diffusive scaling. Successively, the scheme is extended to the case of the kinetic radiative transport equation [14].
The remainder of this paper is organized as follows. In Section 2, we discuss the development of the new Monte Carlo method in the case of the GoldsteinTaylor model. First we introduce the model and its diffusion limit. We also recall the standard splitting method that leads to classical Monte Carlo schemes for kinetic equations. Next, we introduce the new asymptotic preserving Monte Carlo method. Subsequently, in Section 3, we extend our methodology to the case of the radiative transport equation. After the introduction of the kinetic model in the diffusive scaling we again discuss its diffusive limit and the corresponding classical Monte Carlo schemes. Then, we then extend the asymptotic preserving Monte Carlo scheme to the case of the radiative transport. We present several numerical tests and analyze the performance of the new methods in Section 4. Finally, in Section 5, we draw some conclusions and give an outlook to future research directions.
2 Asymptoticpreserving Monte Carlo methods for the Goldstein Taylor model
In this section we discuss the construction of asymptoticpreserving Monte Carlo methods using as a prototype problem the GoldsteinTaylor (GT) model in the diffusive scaling. To this aim, we first introduce the prototype problem and emphasize the drawbacks of a standard Monte Carlo approach. Next, by means of a suitable problem reformulation we construct our novel class of asymptotic preserving Monte Carlo (APMC) methods.
2.1 The GT model in the diffusive limit
The GT model [26] can be written as the following system of kinetic equations
(1) 
where , and . The model describes the evolution of two sets of particles: particles with (constant) velocity , with density , and particles with (constant) velocity with density . Equation (1) gives a physically intuitive description of the process: the left hand side denotes transport with the characteristic speeds , whereas the right hand side encodes random velocity changes that can be interpreted as collisions: in both populations particles disappear with a rate and reappear with the opposite velocity.
The diffusive scaling corresponds to take and where is the scaling parameter. The scaled model reads
(2) 
where we have introduced the mass density and the scaled momentum defined as
(3) 
Now, using the above macroscopic variables system (2) can be written equivalently in the form
(4) 
The original kinetic variables are recovered through relations
(5) 
Throughout the manuscript, we will systematically switch back and forth between the systems (2) and (4). The motivation is that (2) is more convenient for the development of the Monte Carlo solver whereas (4) permits to compute the diffusion limit simply using the leading order term.
2.2 A standard Monte Carlo method for the GT model
To construct the Monte Carlo scheme, we define an ensemble of particles , in which represents the position and the velocity of particle at time . For each , the velocity can only take two values, namely
We will approximate the functions (respectively ), representing the density of particles with velocity (respectively ), by an empirical distribution
(7)  
(8) 
in which represents a Dirac delta, a Kronecker delta, and the mass of an individual particle is defined as
(9) 
Note that one can introduce a (possibly timedependent) weight to each particle, and consider the corresponding weighted empirical distributions. We will not pursue this path here, as doing so does not affect the modified Monte Carlo schemes that are the focus of this manuscript.
Remark 1 (From empirical distributions to spacediscretized particle densities)
Given the ensemble of particles , a approximation to the particle densities can be obtained as an histogram by introducing a spatial mesh with centers and mesh spacing , and defining the stochastic approximation by simply counting the particles inside the bin :
(10) 
The density can then be obtained as
(11) 
Clearly, many alternative approaches are available in the literature, see, e.g., [43, 38]. In particular, a popular alternative to the abovedescribed histogram approach is kernel density estimation [48].
To simulate the GoldsteinTaylor model, one needs to define a stochastic process for the evolution of the ensemble , such that the population dynamics corresponds to (2). We introduce the discretized time , with the time step and , and write the timediscretized particle states as and respectively [45]. The standard Monte Carlo method for the GoldsteinTaylor model (2) is based on a splitting between the transport and collision terms [43, 45]

Transport:
(12) 
Collision:
(13)
The splitting (12)–(13) provides a convenient strategy to define the evolution of the particles: the transport step can be seen to affect the particle positions, leaving the velocities untouched, whereas the collision step updates the velocities, leaving the positions unaltered.
Transport step
During the transport step (12), each particle advances from time over a time step of size by changing its position according to
(14) 
This can easily be shown by inserting (7)–(8) in (12). Thus, after the transport step, we obtain the intermediate empirical distributions,
(15)  
(16) 
from which the intermediate particle densities can be computed via (10).
Collision step
We consider now the solution of the collision step (13). First observe that, during the collision step, the density is constant. Hence, the density after the full time step has already been obtained during the transport step, and we have . The effect of the collision step is thus to randomly change the velocities of a fraction of the particles from to , and vice versa, without affecting their positions. Thus, the density after the complete time step is given as
(17) 
and its stochastic approximation follows from equation (10).
Now observe that, since equation (13) is linear and local, its exact solution is known. When neglecting (for now) the stochastic particle discretizations (10), and using as initial conditions the values of and after the transport step, i.e. , we have
(18)  
At the Monte Carlo level, the above formula can be interpreted as the convex combination of two probability distributions:

With probability , the speed of a particle does not change, and the particle is left untouched.

With probability the speed of a particle changes, and the new velocity is chosen to be or with equal probability.
Remark 2 (Computational complexity in the diffusion limit)
Even if there is no time step restriction in the collision part of the algorithm (the step always represents a convex combination of two distributions, regardless of the size of and ), the approach outlined above suffers of severe time step restrictions when . In fact, the main problem arises in the transport phase of the algorithm: when , the scaled particle velocities diverge, since and . This means that we are forced to take , making the Monte Carlo solver unusable for small values of in the diffusion limit.
While the standard Monte Carlo method suffers from a prohibitive computational cost, we can observe that in the limit, the kinetic equation becomes equivalent to the heat equation (6) and consequently a Monte Carlo method which solves this problem is easily available. In this case, the Monte Carlo method consists in first assigning the positions to particles which approximates the function at the initial time by the empirical measure
(19) 
where the particle positions are sampled from the probability distribution with density and the constant is defined as
(20) 
Successively, the position of the particles evolves in time according to
(21) 
where is a standard normally distributed random number. Observing that the Monte Carlo method (21) does not have time step limitations, one would like to construct a Monte Carlo scheme for the GoldsteinTaylor model (4) which in the limit degenerates to (21) without time step limitations induced by the characteristic speeds. This property is refereed to as AsymptoticPreserving (AP) property and we discuss such a Monte Carlo scheme in the next paragraph.
2.3 An AP Monte Carlo method for the GT model
In this section, we introduce a new Monte Carlo approach based on a suitable reformulation of the original system. We first discuss a time discretization that leads to a reformulated GT model (section 2.3.1). Subsequently, we introduce our new scheme based on the reformulated system (section 2.3.2).
2.3.1 An implicit timediscrete reformulation
We start considering the following fully implicit discretization for system (4)
(22) 
Now, solving the second equation for , one obtains
(23) 
or, equivalently,
(24) 
Plugging (23) in the first equation we get
Finally, using the first equation of (22) and filling this into (24) we get the equivalent form
(25) 
which, using the change of variables (5), can be also written as
(26) 
Observe now that, by using
we obtain that system (25) is equivalent up to first order in to
(27) 
and in diagonal form
(28) 
Note that, the left part of system (28) is hyperbolic with characteristic speeds
(29) 
When for a fixed , system (27) converges to the original system (4), while the characteristic speeds converge to the usual ones, i.e.,
On the other hand, for a fixed , the characteristic speeds and are bounded for any value of and converge to zero as , while the diffusion coefficient tends to in that limit. Consequently, the system becomes fully parabolic and converges to the solution of the heat equation
(30) 
2.3.2 The APMC method
We now introduce a Monte Carlo scheme that solves the GoldsteinTaylor model (4) for all choices of the time step and , without any dependent time step restriction. The method is based on the following splitting approach between:

Transport–diffusion:
(31) 
Collision:
(32)
Now, as we did with the standard Monte Carlo method in Section 2.2, we approximate the functions and by a finite set of particles , which correspond to the empirical measures (7) and (8). The velocities now can take the two values
which are bounded for any value of and are such that as . Recall that we can then reconstruct a histogram approximation of the distributions on a mesh via (10).
Transport–diffusion step
To solve the transport–diffusion step (31), we observe that the particle velocities now scale with , and that the diffusion corresponds to a Brownian motion with coefficient . Thus, particles move according to
(33) 
in which is a standard normally distributed random variable. The velocity of the particles does not change in this step, and we again have the intermediate empirical distributions (15) and (16), from which the intermediate particle densities can be computed via (10).
Collision step
We consider now the collision step (32). As in the standard Monte Carlo method in Section 2.2, the density is constant. The effect of this step is to randomly change the velocities of a fraction of the particles from to , keeping the positions untouched. By solving (32) with the forward Euler method we get
(34) 
This forward Euler discretization leads to the following convex combination
(35) 
Compared with collision step (18) in the standard Monte Carlo method, the only change is a slightly changed collision rate, i.e., a slightly different probability of a velocity change. At the Monte Carlo level, equation (35) is thus again interpreted as

With probability , the speed of a particle does not change, and the particle is left untouched.

With probability the speed of a particle changes, and the new velocity is chosen to be or with equal probability.
It is easy to see the asymptoticpreserving property of the new method. In fact, the time step of the transport–diffusion step is now independent of . In particular, in the limit we get a standard Brownian motion for the heat equation
(36) 
Remark 3

The GoldsteinTaylor model is equivalent to the telegrapher’s equation, other probabilistic approach can be derived using this latter form [1, 31]. Note however, that in the diffusion limit the time step in the above approaches has to be taken of the size of and therefore the methods are not asymptoticpreserving.
3 Asymptoticpreserving Monte Carlo methods for the radiative transport
In this section we show how to generalize the above approach to the radiative transport equation [10, 11] under the diffusive scaling.
3.1 The radiative transport equation
Let be the probability density distribution for particles at space point , at time traveling in direction , with . Particles undergo two types of interactions: scattering, with scattering coefficient and absorption, with absorption coefficient . Under the diffusive scaling, solves the radiative transfer equation
(37) 
where is the total transport coefficient, is the source term, is proportional to the mean free path and
(38) 
is the position density.
To study the process in the diffusive limit when tends to zero, we use the following expansion in of the distribution function
(39) 
and we introduce it in (37). Then, considering terms of the same order in , we get at the leading order
(40) 
where represents the density of the gas. Then, to the first order in , we get
(41) 
Now, writing the balance equation in terms of of (37), one gets
(42) 
and the integration in velocity space yields
(43) 
Now, assuming positive and strictly bounded away from zero, since one has from (41) that
(44) 
we finally obtain the following equation
(45) 
where is the socalled diffusion coefficient which, for example, takes the value in onedimensional slab geometry and when is a unit circle in two dimensions.
3.2 A standard Monte Carlo scheme for the radiative transport
Let us now discuss a standard Monte Carlo method for solving the radiative transport equation, highlighting the limitations of this approach when close to the diffusive limit. The starting point is, as for the twospeed case, a time splitting scheme for (37) (where for simplicity we set the source term ). It reads

Transport:
(46) 
Collision:
(47) 
Absorption:
(48)
We again approximate the distribution by an empirical distribution, using a finite set of particles with positions and velocities . The particle velocities are given as
Defining the mass of an individual particle as
(49) 
we obtain the empirical distribution
(50) 
in which is again the Dirac delta.
Remark 4 (Empirical particle densities in phase space)
We restrict for simplicity to the onedimensional case in both and . In analogy with the twospeed case, we can introduce a mesh and compute a histogram on the phase space mesh with cell centers and and mesh widths and as
(51) 
An empirical position density is then obtained as
(52) 
We can now describe the Monte Carlo method that correspond to transport and collisions of the particles. In the sequel we will neglect the presence of the source term which can be easily included in the method.
Transport
As in the twospeed case, each particle advances from time over a time interval of length during the transport step (46) by changing its position according to
(53) 
We then have an intermediate empirical distribution:
(54) 
from which the intermediate particle density can be computed using (51).
Collision
Next, we solve the collision process (47) without absorption
(55) 
which corresponds to
(56) 
At the Monte Carlo level, the above formula can be interpreted in the following way:

With probability , the speed of a particle does not change

With probability , the speed of a particle changes to a new value , in which is a random value with uniform probability in the domain .
Absorption
We consider now the solution of the absorption step (48). Unlike the twospeed case, due to absorption the density is not conserved.
In a time step , we solve the absorption process
(58) 
which allows to compute
(59) 
The above process is easily realized, assuming that with probability , the particle gets absorbed and disappears from the simulation.
Remark 5 (Meshbased approach)
We may consider a method based on a mesh in space. We define the density of the particles in the center of the cells, and solve (47) in with the number of mesh points. In order to compute the integral of the distribution function in the cell centers different techniques can be used. The simplest first order space reconstruction in one dimension, the same used for the two speed case, is given by (52).
Concerning the asymptotic behavior for small values of the same remark can be made as in the twospeed case (see Remark 2): in the limit when tends to zero, the main computational bottleneck is due to transport phase, where the transport speeds approach infinity and hence infinitely small time steps would be required.
Also for the radiative transfer equation, a standard Monte Carlo method for the diffusion equation (45) can be derived. In the one dimensional case, if we neglect the source term and the absorption, it consists in initializing the system by creating an ensemble of particles that are sampled according to the local density , and then by advancing in a time step following the equation
(60) 
where is a standard normally distributed random number. Thus, we want to construct am asymptoticpreserving Monte Carlo method for radiative transfer (37) that automatically degenerates to the above described Monte Carlo method without any time step restriction induced by the unbounded increasing particle speed when .
3.3 An AP Monte Carlo method for the radiative transport
In this section, we generalize the reformulation discussed in Section 2.3 for the GoldsteinTaylor model (4) to the case of the radiative transfer model (37). We start by reformulating the radiative equation by using the even and odd formalism and by introducing a suitable time discretization (section 3.3.1). In section 3.3.2, we then make use of this reformulation for constructing our new scheme.
3.3.1 The reformulated radiative transport equation
In order to emphasize the analogies with the Goldstein–Taylor model we consider radiative transport equation (37) without source term , i.e.,
(61) 
with
We first rewrite the radiative transfer equation as
(62) 
where now and . This permits to define the even and odd parities
(63) 
Then, and satisfy the following equivalent system
(64) 
The inverse transformation of (63) is easily seen to be
(65) 
In order to construct an implicit reformulation of the problem we first split the system into three parts as
(66) 
(67) 
and
(68) 
The first step (I) now has the same structure of the Goldstein–Taylor model and we can follow the approach developed in Section 2.3.1. We consider the implicit discretization of (66) as
(69) 
where and denote the solutions of this first step.
Solving for one gets
(70) 
or, equivalently,
(71) 
Equation (70) can be plugged in the first equation of (69) to give
(72) 
Now, using the first equation of (69) into (71) gives
(73) 
The second part of the splitting, equation (67), can be discretized similarly with an implicit method to give
(74) 
where since the density remains unchanged during this step. We observe now that (73)(74) are, up to an error , equivalent to a time splitting of the reformulated system
(75) 
Using the back transformation (65), equation (75) can be written also as
(76) 
Note that, for fixed values of , the above equations revert to the original system (64) or (61) in the limit when the time step tends to zero when the absorption coefficient . Let observe that up to an error of order system (75) or (76) plus the third step of the splitting (68) represent a first order in time approximation of the original radiative transfer equation. On the other hand, for all , system (75) or (76) together with (68) are a approximation with bounded eigenvalues for every choice of of the original system. In particular, for every finite time step, the system tends to the limiting diffusion equation (45) in the limit when tends to zero.
Remark 6 (Micromacro decomposition)
As an alternative to the oddeven splitting above, one could also consider a micromacro splitting, see, e.g., [32, 39]. Let us illustrate the approach in one space dimension and with for simplicity. In that case, we write
(77) 
with defined as before, from which we naturally derive that . Inserting this expansion in (61) and averaging over velocity space, we get the system