# Addressing the gas kinetics Boltzmann equation with branching paths
statistics

Guillaume Terrée  Mouna El Hafi  RAPSODEE, UMR CNRS 5302, Mines Albi, Campus
Jarlard; 81013 Albi CT cedex 09, France  Stéphane Blanco  Richard Fournier
LAPLACE, Université de Toulouse, CNRS, INPT, UPS, France  Jérémi Dauchet
Université Clermont Auvergne, CNRS, SIGMA Clermont, Institut Pascal, F-63000
Clermont–Ferrand, France  Jacques Gautrais  Centre de Recherches sur la
Cognition Animale (CRCA), Centre de Biologie Intégrative (CBI), Université de
Toulouse, CNRS, UPS; 118 Route de Narbonne, Toulouse F-31062 cedex 4, France

###### Abstract

Several works have shown that it is possible to conceive non-linear Monte-
Carlo algorithms, but these works have remained overlooked so far. It is
commonly admitted that in non-linear physics it is necessary to start with a
linearization step before then calling on the linear Monte-Carlo method. We
show here how it is possible to avoid this linearization step in the
particular case of the Boltzmann equation. Our method retains the exact
character of the Monte-Carlo method (unbiased estimation with a reliable
evaluation of its uncertainty). We also extend to non-linear transport a very
strong property of Monte-Carlo algorithms in linear transport: the ability to
perform probe calculus (evaluate a quantity at a given position and time
without computing all the field at previous times). Another striking feature
of this algorithm is its ability to focus on rare events, however scarce they
may be, including very high energy events that are important from the
applicative point of view.

Monte Carlo method; branching random process; integral formulation; unbiased
estimator; non-linear Boltzmann equation; gas kinetics

##  I Introduction

### In linear transport physics,

and especially in radiative transfer physics Howell (1998); Metropolis and
Ulam (1949), the Monte-Carlo method (MCM) consists in simulating the history
of numerous independent particles, from which mean observables can be deduced.
Since it’s initial development, this numerical method has become extremely
successful, because of its many qualities: a natural management of multi-
dimensional phase space, a null systematic error compared to the mathematical
and physical model, the confidence intervals given with the results, an
ability to take into account numerous physical phenomena simultaneously, an
efficient management of complex geometries, the abilty to calculate
sensitivities simultaneously , and easy parallelization Dunn and Shultis
(2012); Dimov (2008). This success extends to all linear physics, not only to
transport Karatzas and Shreve (1991); Feynman  _et al._ (2010); Umrigar  _et
al._ (2007).

### In non-linear physics,

the Monte-Carlo method is commonly believed to be unusable:

> “So far as the author is aware, the extension of Monte Carlo methods to
> nonlinear processes has not yet been accomplished and may be impossible.”
> Curtiss (1954)
>
> “Monte Carlo methods are not generally effective for nonlinear problems
> mainly because expectations are linear in character.” Kalos and Whitlock
> (2008)
>
> “A nonlinear problem must usually be linearized in order to use Monte Carlo
> technique.” Kalos and Whitlock (2008)

However several examples lead us to moderate these assertions.

A first group of such examples is the lineage of the tests of Ermakov _et al_
on the MCM used for solving non-linear kinetic equations Ermakov  _et al._
(1976). At least twenty-five years ago, Gurov Gurov (1992); Dimov and Gurov
(2000) presented the calculation in the MCM of Fredholm integrals containing
polynomial non-linearities, by branching estimation processes. Rasulov _et al_
Rasulov  _et al._ (2004, 2010); Raimova (2016) showed how to use the MCM to
resolve parabolic equations like the equation of heat, with a source term
depending non-linearly on the density. Although barely cited or reused Henry-
Labordère  _et al._ (2014), these techniques were shown to work by their
authors, who do not hide their intention to extend the MCM towards a broader
class of non-linear problems.

A second group of examples Dauchet  _et al._ (2016) comes from the recent PhD
works of J. Dauchet Dauchet (2012) and O. Farges Farges (2014), dedicated
respectively to radiative transfer in photobioreactors and to the optimal
design of concentrated solar power plants. The authors were confronted with
numerical problems in which, on the one hand the geometrical complexity was so
high that only statistical methods were looking affordable De la Torre  _et
al._ (2014), and on the other hand non-linearities were present – especially
in the coupling law between the local productivity inside the set-up and the
radiative transfer. Dauchet has also tackled the problem of evaluating the
radiative properties of microalgae by means of electromagnetics, accounting
for their size and shape distribution Dauchet  _et al._ (2015): here the non-
linearity came from the fact that radiative properties are quadratic in the
electric field amplitude, which he was basically able to calculate Charon  _et
al._ (2015). Confronted with these difficulties, J. Dauchet developed a
technique to handle multi-folded integrals with one single stage of non-
linearity in between Dauchet  _et al._ (2016), drawing his inspiration from
Null Collision Algorithms (NCAs).

NCAs are the latest example we know of handling non-linearities in the MCM,
although the example is drawn from linear transport – _eg_ radiative transfer
physics Galtier  _et al._ (2013) and image synthesis Novák  _et al._ (2014);
Szirmay-Kalos  _et al._ (2017), semiconductor physics Rees (1968), neutronics
Woodcock  _et al._ (1965); Spanier (1966); Coleman (1968), or even plasma
physics Skullerud (1968). These algorithms serve to account for the extinction
properties of a medium, when they are highly variable. Their principle is to
add virtual and ineffective colliders in order to make constant the total
collision frequency. The NCAs are justified by the rigorous equivalence
between the original and the modified transport problems; but they can also be
viewed as a way to handle the non-linearity of the exponential Beer extinction
law Longo (2002).

### A trick

shared in all these works is to convert the non-linearity into an increase in
the dimension of the sampling space Dauchet  _et al._ (2016).

As a theoretical illustration, let us write that a searched quantity $q$ is
expressed as an expectation (either by mimicry of the physical process under
study, or through a more mathematical reasoning) using the random variable
(RV) ${\vec{X}}$, under the form:

$\displaystyle q$
$\displaystyle=\operatorname{E}\Bigl{(}a({\vec{X}})h({\vec{X}})+b({\vec{X}})\Bigr{)}$
(1a) where $\operatorname{E}$ denotes the statistical expectation, and $a$ and
$b$ are two deterministic functions. Suppose also that $h({\vec{X}})$ is only
known as the expectation of an estimator $\tilde{h}({\vec{X}};{\vec{Y}})$
using the sampling of a RV ${\vec{Y}}$: $\displaystyle q$
$\displaystyle=\operatorname{E}\Bigl{(}a({\vec{X}})\operatorname{E}_{|{\vec{X}}}\bigl{(}\tilde{h}({\vec{X}};{\vec{Y}})\bigr{)}+b({\vec{X}})\Bigr{)}$
(1b) the subscript $\bullet_{|{\vec{X}}}$ meaning that the value of
${\vec{X}}$ is known. As the statistical expectation is linear and projective:
$\displaystyle q$
$\displaystyle=\operatorname{E}\Bigl{(}a({\vec{X}})\tilde{h}({\vec{X}};{\vec{Y}})+b({\vec{X}})\Bigr{)}$
(1c)

It is thus useless to evaluate $h({\vec{X}})$ exactly (_ie_ without
statistical noise) at each value of ${\vec{X}}$ in order to estimate $q$
correctly (_ie_ without bias, in MCM). Conversely, an algorithm estimating
$h({\vec{x}})$ for any ${\vec{x}}$ can be turned into an algorithm estimating
$q$, simply by adding a first sampling of ${\vec{X}}$. This is one of the
principal advantages of the Monte-Carlo method: its complexity increases
linearly (and not exponentially) with the total dimension of the configuration
space – here configurations are described by the couple of RVs
$({\vec{X}};{\vec{Y}})$. But this becomes false if the coupling law $h\to q$
is non-linear. Imagine now that it is a square, then in a general manner

$\operatorname{E}\Bigl{(}a({\vec{X}})h({\vec{X}})^{2}\Bigr{)}\neq\operatorname{E}\Bigl{(}a({\vec{X}})\tilde{h}({\vec{X}};{\vec{Y}})^{2}\Bigr{)}$
(2)

the same way as $(u+v)^{2}\neq u^{2}+v^{2}$.

However, in this quadratic case, we could replace ${\vec{Y}}$ with two
independent and identically distributed RVs, ${\vec{Y}}\\!_{1}$ and
${\vec{Y}}\\!_{2}$, and then consider that:

$\operatorname{E}\Bigl{(}a({\vec{X}})h({\vec{X}})^{2}\Bigr{)}=\operatorname{E}\Bigl{(}a({\vec{X}})\tilde{h}({\vec{X}};{\vec{Y}}\\!_{1})\tilde{h}({\vec{X}};{\vec{Y}}\\!_{2})\Bigr{)}$
(3)

The same technique is usable with monomial and polynomial functions, as:

$\forall
n\in\mathbb{N},\operatorname{E}\Bigl{(}a({\vec{X}})h({\vec{X}})^{n}\Bigr{)}=\operatorname{E}\Bigl{(}a({\vec{X}}){\textstyle\prod_{i=1}^{n}\tilde{h}({\vec{X}};{\vec{Y}}\\!_{i})}\Bigr{)}$
(4)

with the ${\vec{Y}}\\!_{i}$ all independent and identically distributed. This
also extends to analytic functions, because in a Monte-Carlo calculation an
infinite sum is no more than an integral.

### Despite the non-linearity of the gas kinetics Boltzmann equation,

the Monte-Carlo method has been commonly used to solve it. But the trick
revealed in the previous paragraph has never been used, and the non-linearity
has always been circumvented by using linearization techniques based on
discretization of time and/or space. The numerical schemes obtained share the
following operating principles: they build the story of a (big) swarm of
statistically representative molecules, and because these cannot cross
exactly, collisions are introduced using a proximity criterion Bird (1994,
1998); Oran  _et al._ (1998); Murrone and Villedieu (2011). There are a lot of
variants being used today, depending on the collision quadrature used in time
and space Khisamutdinov (2004); Stefanov (2011), or on whether an Eulerian
fluid model is used in tandem with the statistical model Homolle and
Hadjiconstantinou (2007); Degond  _et al._ (2010)…

### The author’s proposal is

to solve the Boltzmann equation with Monte-Carlo algorithms like those used in
linear transport, without linearization or discretization, using the trick
revealed in the previous paragraph along with the null collision technique.
The trick alone will deal with the source term in the Boltzmann collision
operator, and the null collisions are useful to handle the extinction term.

### What we have obtained

is a statistical numerical method in gas kinetics, the statistical
formulations of which are given below. The method retains the advantages of
the MCM from linear transport, its most striking properties being that:

  * •

No mesh or time discretization is necessary in the method.

  * •

The method enables probe calculus, _ie_ it is possible to calculate the
distribution function at a point of the phase space, without computing the
rest of the field.

  * •

Rarefaction of the gas at the probe does not compromise the relative precision
of the calculations.

  * •

Similarly, the frequency of rare events (in the space of speeds) can be
estimated however scarce they may be.

  * •

The method is limited in (spatial or temporal) Knudsen number; because it uses
a branching estimation process, the mean complexity of which grows
exponentially with the mean number of collisions.

  * •

The branchings in the estimation process also create a variance problem, which
is very sensitive to algorithmic choices. Today, a theoretical framework for
understanding this is lacking.

  * •

As a result of the two previous points, we are today essentially unable to
perform a calculus on a steady state.

### This paper is organized as follows.

In section II, we show how to handle the non-linearity of the source term in
the Boltzmann collision operator, in an academic case called the BKW mode
Hauge and Præstgaard (1981) which has a symbolic expression. In section III,
we shall see how to deal with extinction in the Boltzmann collision operator,
in a very simplified toy model with a symbolic solution. In sections IV and V
we deal with the full Boltzmann equation, firstly in a particular case Guéry-
Odelin  _et al._ (2014) with a symbolic solution (and local equilibrium) and
secondly in a more general case. We conclude in section VI.

##  II Example: the BKW mode

### The Bobylev-Krook-Wu mode

was discovered simultaneously by these authors in 1976 Bobylev (1976); Krook
and Wu (1976). It is an explicit solution of the Boltzmann equation,
associated with a particular set of physical conditions.

We place ourselves in ordinary 3-dimensional space 111The BKW mode can be
generalized to any dimensionality Ernst (1979).. The collision model consists
of Maxwell molecules with isotropic scattering: the differential cross section
${\sigma_{F}}$ is inversely proportional to the relative speed $g$, with no
angular dependence. The gas is uniform, with a unit density and a peculiar
velocity of zero. Thus $f\equiv f(\vec{c};t)$, and with the adequate non-
dimensionalization of time the Boltzmann equation reads:

$\partial_{t}f(\vec{c};t)=\rule{0.0pt}{4.30554pt}-f(\vec{c};t)+\int_{E_{c}}{\mathrm{d}}{\vec{c}_{*}}\int_{E_{u}}\dfrac{{\mathrm{d}}{\vec{u}^{\mspace{1.5mu}\prime}}}{4\pi}\,f({\vec{c}^{\mspace{3.5mu}\prime}};t)f({\vec{c}^{\mspace{3.5mu}\prime}_{*}};t)$
(5)

where ${E_{c}}$ and ${E_{u}}$ denote respectively the spaces of velocities and
directions (so ${E_{u}}$ is the unit sphere of ${E_{c}}$), and
${\vec{c}^{\mspace{3.5mu}\prime}}$ and ${\vec{c}^{\mspace{3.5mu}\prime}_{*}}$
state

$\left\\{\begin{array}[]{@{\
}l}{\vec{c}^{\mspace{3.5mu}\prime}}=\tfrac{1}{2}(\vec{c}+{\vec{c}_{*}}+g{\vec{u}^{\mspace{1.5mu}\prime}})\\\
\\\\[-4.30554pt]
{\vec{c}^{\mspace{3.5mu}\prime}_{*}}=\tfrac{1}{2}(\vec{c}+{\vec{c}_{*}}-g{\vec{u}^{\mspace{1.5mu}\prime}})\\\
\\\\[-4.30554pt] g=\|\vec{c}-{\vec{c}_{*}}\|\end{array}\right.$ (6)

In this case, the Boltzmann equation admits the particular solution known as
BKW mode

$f(\vec{c};t)=\dfrac{\exp\bigl{(}-\vec{c}^{\>2}/(2K{c_{q}}^{2})\bigr{)}}{2\bigl{(}\sqrt{2\pi
K}{c_{q}}\bigr{)}\rule{0.0pt}{4.30554pt}^{3}}\biggl{(}\dfrac{5K-3}{K}+\dfrac{1-K}{K^{2}}\times\dfrac{\vec{c}^{\>2}}{{c_{q}}^{2}}\biggr{)}$
(7a) with $K=1-\tfrac{2}{5}\exp(-t/6)$ (7b)

and ${c_{q}}$ the Root Mean Square velocity on each axis (RMS velocity). With
the $\frac{2}{5}$ factor in the expression of $K$, $f$ is definite and
positive for $t\geqslant 0$. Fig. 1 gives a graphical overview of the BKW
mode.

(a) Linear scale on $y$-axis

(b) Logarithmic scale on $y$-axis

Figure 1: The distribution function $f$ along the $(Ox)$ axis of the velocity
space, according to the BKW mode described in Eq. (7), at several times $t$.
$\vec{u\mspace{-1.0mu}_{x}}$ denotes the $x$ unit vector.

### In order to build a Monte-Carlo calculation

in the physical situation of the BKW mode, let us rewrite the evolution Eq.
(5).

$\partial_{t}f(\vec{c};t)=\rule{0.0pt}{4.30554pt}-f(\vec{c};t)+{s_{\mathit{is}}}(\vec{c};t)$
(8a) with
${s_{\mathit{is}}}(\vec{c};t)=\int_{E_{c}}{\mathrm{d}}{\vec{c}_{*}}\int_{E_{u}}\dfrac{{\mathrm{d}}{\vec{u}^{\mspace{1.5mu}\prime}}}{4\pi}\,f({\vec{c}^{\mspace{3.5mu}\prime}};t)f({\vec{c}^{\mspace{3.5mu}\prime}_{*}};t)$
(8b)

If it is assumed that the initial condition is “$f$ is known at $t=0$”, the
time-differential Eq. (8a) can be formally solved, giving:

$\forall t\geqslant
0,f(\vec{c};t)=\int_{-\infty}^{t}{\mathrm{d}}t^{\prime}\,\exp(-(t-t^{\prime}))\times\rule{0.0pt}{4.30554pt}\\\
\\\\[-8.61108pt]
\Bigl{(}\operatorname{H}(t^{\prime}){s_{\mathit{is}}}(\vec{c};t^{\prime})+\operatorname{H}(-t^{\prime})f(\vec{c};0)\Bigr{)}$
(9)

where $\operatorname{H}$ is the Heaviside step function:
$\operatorname{H}(t)=1$ if $t>0$, $\operatorname{H}(t)=0$ if $t<0$.

Together, Eqs. (9) and (8b) give a second kind Fredholm writing of $f$. By
choosing probability densities in order to sample $t^{\prime}$,
${\vec{c}_{*}}$, and ${\vec{u}^{\mspace{1.5mu}\prime}}$, this can be converted
into the Monte-Carlo algorithm 1, intended to estimate $f(\vec{c};t)$ with
$(\vec{c};t)\in{E_{c}}\times\mathbb{R}^{+}$. Because algorithm 1 is recursive,
its ability to complete will be known when setting how $t^{\prime}$,
${\vec{c}_{*}}$, and ${\vec{u}^{\mspace{1.5mu}\prime}}$ are sampled.

Input: A point $(\vec{c};t)$ in ${E_{c}}\times\mathbb{R}^{+}$

Output: A point estimate of $f(\vec{c};t)$

Sample $T^{\prime}$ in $(-\infty;t]$: $t^{\prime}$ is obtained;

if  _$t^{\prime}\leqslant 0$_ then

return
_$\dfrac{\exp(-(t-t^{\prime}))\,f(\vec{c};0)}{{p_{T^{\prime}}}(t^{\prime})}$_
;

// $f(\vec{c};0)$ is known

else

Sample ${\vec{C}_{*}}$: ${\vec{c}_{*}}$ is obtained;

Sample ${\vec{U}^{\prime}}$: ${\vec{u}^{\mspace{1.5mu}\prime}}$ is obtained;

${\vec{c}^{\mspace{3.5mu}\prime}}\leftarrow\tfrac{1}{2}\bigl{(}\vec{c}+{\vec{c}_{*}}+\|\vec{c}-{\vec{c}_{*}}\|\,{\vec{u}^{\mspace{1.5mu}\prime}}\bigr{)}$;

${\vec{c}^{\mspace{3.5mu}\prime}_{*}}\leftarrow\tfrac{1}{2}\bigl{(}\vec{c}+{\vec{c}_{*}}-\|\vec{c}-{\vec{c}_{*}}\|\,{\vec{u}^{\mspace{1.5mu}\prime}}\bigr{)}$;

Estimate $f({\vec{c}^{\mspace{3.5mu}\prime}};t^{\prime})$ using this
algorithm: $\widetilde{f}_{1}({\vec{c}^{\mspace{3.5mu}\prime}};t^{\prime})$ is
obtained;

Estimate $f({\vec{c}^{\mspace{3.5mu}\prime}_{*}};t^{\prime})$ using this
algorithm:
$\widetilde{f}_{2}({\vec{c}^{\mspace{3.5mu}\prime}_{*}};t^{\prime})$ is
obtained;

// independently to
$\widetilde{f}_{1}({\vec{c}^{\mspace{3.5mu}\prime}};t^{\prime})$

return
_$\dfrac{\exp(-(t-t^{\prime}))\,\widetilde{f}_{1}({\vec{c}^{\mspace{3.5mu}\prime}};t^{\prime})\,\widetilde{f}_{2}({\vec{c}^{\mspace{3.5mu}\prime}_{*}};t^{\prime})}{4\pi\,{p_{T^{\prime}}}(t^{\prime})\,{p_{{\vec{C}_{*}}}}({\vec{c}_{*}})\,{p_{{\vec{U}^{\prime}}}}({\vec{u}^{\mspace{1.5mu}\prime}})}$_
;

Algorithm 1 Algorithm for estimating $f(\vec{c};t)$, valid in the conditions
of the BKW mode described in section II (uniform and unit density, Maxwell
molecules…)

### The recursive character of algorithm 1

is a feature common to all Monte-Carlo algorithms based on second kind
Fredholm formulations. It is found in all multiple scattering problems of
transport physics. But contrary to what happens in linear transport, each
recursive stage of an estimation of $f$ by algorithm 1 is likely to call
itself _more than once_. So, to be exact, the estimation process no longer
consists in following a (virtual) particle path, but rather in following a
particle _tree_.

We began the formal development of our proposal with a recursive integral
formulation, rather than with a path integral. We will continue in this way,
and invoke the notion of path (or tree) as little as possible. However, for
the sake of understanding, it seems useful to illustrate here what the
branching estimation paths will look like.

At the same time, we propose a way to locate in such estimation trees. We take
our inspiration from a previous proposal of Gurov Gurov (1992); Dimov and
Gurov (2000); Moradi  _et al._ (2014). Any intermediate estimate of $f$ in a
estimation tree is indexed with an integer sequence, according to the
following rules:

  * •

Only the sequences $\mathbb{N}^{*}\to\mathbb{N}$, stationary at 0, in which
all terms are null since the first null term, are used.

  * •

The degree of the sequence (the index of its furthest non-zero term) indicates
the depth of the estimate inside the estimation tree; _ie_ the number of
$(T^{\prime};{\vec{C}_{*}};{\vec{U}^{\prime}})$ samplings which lead to the
point where the considered intermediary estimate is needed.

Thus, the final estimate given by the first call of algorithm 1 is indexed by
the null sequence.

  * •

Among the non-zero terms, the $i$o term equals 1 if at the $i$o depth of
recursion the branch of estimating ${\vec{c}^{\mspace{3.5mu}\prime}}$ has been
followed, while it equals 2 if the branch of estimating
${\vec{c}^{\mspace{3.5mu}\prime}_{*}}$ has been followed.

An example of an estimation tree, with this location indexing, is given in
Fig. 2.

(a) Display of an estimation tree, in one velocity coordinate $c_{y}$ _versus_
the time $t$.

(b) Display of an estimation tree, projected on the $(xOy)$ plane of the
ordinary space.

Figure 2: Illustration of an estimation tree as generated by algorithm 1, with
locations indexed as proposed in paragraph II. The same tree is displayed, in
(a) with one velocity coordinate _versus_ time, in (b) with 2 ordinary space
coordinates. For the sake of readability, the speeds have been indexed
identically to the estimations that are made there. The collision locations
are denoted in (b) by circles and in (a) by rounded cap vertical lines.
Squares denote where the tree meets the initial condition.

### We propose

the following choice for sampling $t^{\prime}$, ${\vec{c}_{*}}$, and
${\vec{u}^{\mspace{1.5mu}\prime}}$:

$\left\\{\begin{array}[]{@{\
}c}{p_{T^{\prime}}}:\left\\{\begin{array}[]{@{}l}(-\infty;t]\to\mathbb{R}^{+}\\\
t^{\prime}\mapsto\exp(-(t-t^{\prime}))\end{array}\right.\\\ \\\\[-4.30554pt]
{p_{{\vec{C}_{*}}}}:\left\\{\begin{array}[]{@{}l}{E_{c}}\to\mathbb{R}^{+}\\\
{\vec{c}_{*}}\mapsto(\sqrt{2\pi}{c_{q}})^{-3}\exp\bigl{(}-{\vec{c}_{*}}\rule{0.0pt}{4.30554pt}^{2}/(2\,{c_{q}}^{2})\bigr{)}\end{array}\right.\\\
\\\\[-4.30554pt]
{p_{{\vec{U}^{\prime}}}}:\left\\{\begin{array}[]{@{}l}{E_{u}}\to\mathbb{R}^{+}\\\
{\vec{u}^{\mspace{1.5mu}\prime}}\mapsto(4\pi)^{-1}\end{array}\right.\end{array}\right.$
(10)

$T^{\prime}$ follows a unit exponential law to the left of $t$,
${\vec{C}_{*}}$ has a Maxwellian distribution with null peculiar velocity and
${c_{q}}$ RMS speed, and ${\vec{U}^{\prime}}$ is isotropically distributed.
Algorithm 1 with these probability densities will be named hereafter algorithm
2: it is displayed in a separate figure.

Input: A point $(\vec{c};t)$ in ${E_{c}}\times\mathbb{R}^{+}$

Output: A point estimate of $f(\vec{c};t)$

Sample $T^{\prime}$ following a unit exponential law to the left of $t$:
$t^{\prime}$ is obtained;

if  _$t^{\prime}\leqslant 0$_ then return _$f(\vec{c};0)$_ ;

// $f(\vec{c};0)$ is known else

Sample ${\vec{C}_{*}}$ following a Maxwellian distribution, with null peculiar
velocity and ${c_{q}}$ RMS speed: ${\vec{c}_{*}}$ is obtained;

Sample ${\vec{U}^{\prime}}$ isotropically: ${\vec{u}^{\mspace{1.5mu}\prime}}$
is obtained;

${\vec{c}^{\mspace{3.5mu}\prime}}\leftarrow\tfrac{1}{2}\bigl{(}\vec{c}+{\vec{c}_{*}}+\|\vec{c}-{\vec{c}_{*}}\|\,{\vec{u}^{\mspace{1.5mu}\prime}}\bigr{)}$;

${\vec{c}^{\mspace{3.5mu}\prime}_{*}}\leftarrow\tfrac{1}{2}\bigl{(}\vec{c}+{\vec{c}_{*}}-\|\vec{c}-{\vec{c}_{*}}\|\,{\vec{u}^{\mspace{1.5mu}\prime}}\bigr{)}$;

Estimate $f({\vec{c}^{\mspace{3.5mu}\prime}};t^{\prime})$ using this
algorithm: $\widetilde{f}_{1}({\vec{c}^{\mspace{3.5mu}\prime}};t^{\prime})$ is
obtained;

Estimate $f({\vec{c}^{\mspace{3.5mu}\prime}_{*}};t^{\prime})$ using this
algorithm:
$\widetilde{f}_{2}({\vec{c}^{\mspace{3.5mu}\prime}_{*}};t^{\prime})$ is
obtained;

// independently of
$\widetilde{f}_{1}({\vec{c}^{\mspace{3.5mu}\prime}};t^{\prime})$

return
_$(\sqrt{2\pi}{c_{q}})^{3}\exp\bigl{(}{\vec{c}_{*}}\rule{0.0pt}{4.30554pt}^{2}/(2\,{c_{q}}^{2})\bigr{)}\,\widetilde{f}_{1}({\vec{c}^{\mspace{3.5mu}\prime}};t^{\prime})\,\widetilde{f}_{2}({\vec{c}^{\mspace{3.5mu}\prime}_{*}};t^{\prime})$_
;

Algorithm 2 Algorithm for estimating $f(\vec{c};t)$, valid in the conditions
of the BKW mode described in section II (uniform and unit density, Maxwell
molecules…). The sampling choices listed in Eq. (10) are used.

This choice, although guided by a desire for simplicity, brings two
interesting qualities to algorithm 1 (now 2):

  1. 1.

The expectation of the recursion, defined as the number of $T^{\prime}$
samplings per final estimation (_ie_ the number of lines in an estimation
tree), is finite. This implies that algorithm 2 surely completes.

  2. 2.

If at $t=0$ the distribution $f$ is at equilibrium:

$f(\vec{c};t=0)=\dfrac{\exp\bigl{(}-\vec{c}^{\>2}/(2\,{c_{q}}^{2})\bigr{)}}{\bigl{(}\sqrt{2\pi}{c_{q}}\bigr{)}\rule{0.0pt}{4.30554pt}^{3}}$
(11)

then algorithm 2 has a null variance.

(a) Results obtained at the centre of the velocity space

(b) Results obtained at $3\times$ the RMS speed away from

the centre of the velocity space

Figure 3: Results obtained by algorithm 2, applied to the BKW mode described
in Eq. (7). Two probe points in the velocity space are considered, where we
follow the distribution function $f$ along the time $t$. The graphs at the top
give a graphical comparison between the results (given with confidence
intervals of 1 standard deviation) and the expected values. The graphs at the
bottom show, with the same results, the relative standard deviation (circles)
(given by the calculations) beside the error actually made (triangles). Each
displayed point was obtained through running e4 realizations of algorithm 2.

(a) Results obtained after 1 mean collision time

(b) Results obtained after 3 mean collision times

Figure 4: Results obtained by algorithm 2, applied to the BKW mode described
in Eq. (7). The probe points are spread on the $(Ox)$ axis of the velocity
space, and only two dates are considered. The graphs at the top give a
graphical comparison between the results (given with confidence intervals of 1
standard deviation) and the expected values. The graphs at the bottom show,
with the same results, the relative standard deviation (circles) (given by the
calculations) beside the error actually made (triangles). Each displayed point
was obtained through running e4 realizations of algorithm 2.

###### Principle of proof (1st point).

Let us consider the expectation of the recursion of algorithm 2 starting at
$(\vec{c};t)$, defined as the expectation of the number of times $T^{\prime}$
random variables will be sampled. This defines a function
${\mathit{rec}}:{E_{c}}\times\mathbb{R}^{+}\to[1;+\infty]$.

Reading algorithm 1, an integral equation on ${\mathit{rec}}$ can be obtained
(as for the calculated quantity $f$). It reads:

${\mathit{rec}}(\vec{c};t)=\int_{-\infty}^{t}{p_{T^{\prime}}}(t^{\prime}){\mathrm{d}}t^{\prime}\biggl{(}1+\operatorname{H}(-t^{\prime})\times
0+\rule{0.0pt}{4.30554pt}\\\
\operatorname{H}(t^{\prime})\int_{E_{c}}{p_{{\vec{C}_{*}}}}({\vec{c}_{*}}){\mathrm{d}}{\vec{c}_{*}}\int_{E_{u}}{p_{{\vec{U}^{\prime}}}}({\vec{u}^{\mspace{1.5mu}\prime}}){\mathrm{d}}{\vec{u}^{\mspace{1.5mu}\prime}}\times\rule{0.0pt}{4.30554pt}\\\
\Bigl{(}{\mathit{rec}}({\vec{c}^{\mspace{3.5mu}\prime}};t^{\prime})+{\mathit{rec}}({\vec{c}^{\mspace{3.5mu}\prime}_{*}};t^{\prime})\Bigr{)}\biggr{)}$
(12)

Then the sampling choices we have made in Eqs. (10) are entered, which gives:

${\mathit{rec}}(\vec{c};t)=\int_{-\infty}^{t}\exp(-(t-t^{\prime})){\mathrm{d}}t^{\prime}\biggl{(}1+\rule{0.0pt}{4.30554pt}\\\
\operatorname{H}(t^{\prime})\int_{E_{c}}\dfrac{\exp\bigl{(}-{\vec{c}_{*}}\rule{0.0pt}{4.30554pt}^{2}/(2\,{c_{q}}^{2})\bigr{)}{\mathrm{d}}{\vec{c}_{*}}}{\bigl{(}\sqrt{2\pi}{c_{q}}\bigr{)}^{3}}\int_{E_{u}}\dfrac{{\mathrm{d}}{\vec{u}^{\mspace{1.5mu}\prime}}}{4\pi}\times\rule{0.0pt}{4.30554pt}\\\
\Bigl{(}{\mathit{rec}}({\vec{c}^{\mspace{3.5mu}\prime}};t^{\prime})+{\mathit{rec}}({\vec{c}^{\mspace{3.5mu}\prime}_{*}};t^{\prime})\Bigr{)}\biggr{)}$
(13)

In order to solve this second kind Fredholm expression manually, it can be
derived with respect to $t$, which gives a differential system on
${\mathit{rec}}$:

$\left\\{\begin{array}[]{@{\ }l@{}}\forall t\geqslant
0,\left\\{\begin{multlined}\partial_{t}{\mathit{rec}}(\vec{c};t)=\rule{0.0pt}{4.30554pt}-{\mathit{rec}}(\vec{c};t)+1+\rule{0.0pt}{4.30554pt}\\\
\displaystyle\int_{E_{c}}\dfrac{\exp\bigl{(}-{\vec{c}_{*}}\rule{0.0pt}{4.30554pt}^{2}/(2\,{c_{q}}^{2})\bigr{)}{\mathrm{d}}{\vec{c}_{*}}}{\bigl{(}\sqrt{2\pi}{c_{q}}\bigr{)}^{3}}\int_{E_{u}}\dfrac{{\mathrm{d}}{\vec{u}^{\mspace{1.5mu}\prime}}}{4\pi}\times\rule{0.0pt}{4.30554pt}\\\
\Bigl{(}{\mathit{rec}}({\vec{c}^{\mspace{3.5mu}\prime}};t)+{\mathit{rec}}({\vec{c}^{\mspace{3.5mu}\prime}_{*}};t)\Bigr{)}\end{multlined}\partial_{t}{\mathit{rec}}(\vec{c};t)=\rule{0.0pt}{4.30554pt}-{\mathit{rec}}(\vec{c};t)+1+\rule{0.0pt}{4.30554pt}\\\
\displaystyle\int_{E_{c}}\dfrac{\exp\bigl{(}-{\vec{c}_{*}}\rule{0.0pt}{4.30554pt}^{2}/(2\,{c_{q}}^{2})\bigr{)}{\mathrm{d}}{\vec{c}_{*}}}{\bigl{(}\sqrt{2\pi}{c_{q}}\bigr{)}^{3}}\int_{E_{u}}\dfrac{{\mathrm{d}}{\vec{u}^{\mspace{1.5mu}\prime}}}{4\pi}\times\rule{0.0pt}{4.30554pt}\\\
\Bigl{(}{\mathit{rec}}({\vec{c}^{\mspace{3.5mu}\prime}};t)+{\mathit{rec}}({\vec{c}^{\mspace{3.5mu}\prime}_{*}};t)\Bigr{)}\right.\\\
\\\\[0.0pt] {\mathit{rec}}(\vec{c};t=0)=1\end{array}\right.$ (14)

${\mathit{rec}}$ is uniform at $t=0$, and its derivative preserves this
uniformity: so here ${\mathit{rec}}(\vec{c};t)\equiv{\mathit{rec}}(t)$.
Finally the system is easily solved in:

${\mathit{rec}}(\vec{c};t)=2\exp(t)-1$ (15)

which is finite for all $(\vec{c};t)\in{E_{c}}\times\mathbb{R}^{+}$. ∎

###### Principle of proof (2nd point).

We call $\widetilde{F}(\vec{c};t)$ an estimator of $f(\vec{c};t)$, obtained
through algorithm 2. We assume that the initial distribution $f(\vec{c};t=0)$
is the equilibrium distribution described in Eq. (11).

Now consider the following induction hypothesis, to be applied anywhere in an
estimation tree:
$\widetilde{F}(\vec{c};t)=(\sqrt{2\pi}{c_{q}})^{-3}\exp\bigl{(}-\vec{c}^{\>2}/(2\,{c_{q}}^{2})\bigr{)}$,
_ie_ $f$ is always estimated as the equilibrium distribution with no variance.

Starting algorithm 2 there are two possibilities: $T^{\prime}\leqslant 0$ and
we are at a leaf of the tree, or $T^{\prime}>0$ and we are at a non-leaf node
of the tree.

If $T^{\prime}\leqslant 0$ then:

$\displaystyle\widetilde{F}(\vec{c};t)$ $\displaystyle=f(\vec{c};0)$
$\displaystyle=\dfrac{\exp\bigl{(}-\vec{c}^{\>2}/(2\,{c_{q}}^{2})\bigr{)}}{\bigl{(}\sqrt{2\pi}{c_{q}}\bigr{)}\rule{0.0pt}{4.30554pt}^{3}}$
(16)

the induction hypothesis is valid.

If $T^{\prime}>0$ then:

$\widetilde{F}(\vec{c};t)=\dfrac{\bigl{(}\sqrt{2\pi}{c_{q}}\bigr{)}\rule{0.0pt}{4.30554pt}^{3}\,\widetilde{F}_{1}({\vec{C}^{\prime}};T^{\prime})\,\widetilde{F}_{2}({\vec{C}^{\prime}_{*}};T^{\prime})}{\exp\bigl{(}-{\vec{C}_{*}}\rule{0.0pt}{4.30554pt}^{2}/(2\,{c_{q}}^{2})\bigr{)}}$

where ${\vec{C}^{\prime}}$ and ${\vec{C}^{\prime}_{*}}$ are the random
variables (RVs) corresponding to the values ${\vec{c}^{\mspace{3.5mu}\prime}}$
and ${\vec{c}^{\mspace{3.5mu}\prime}_{*}}$ in algorithm 2. If the induction
hypothesis holds for $\widetilde{F}_{1}({\vec{C}^{\prime}};T^{\prime})$ and
$\widetilde{F}_{2}({\vec{C}^{\prime}_{*}};T^{\prime})$, then:

$\displaystyle\widetilde{F}(\vec{c};t)$
$\displaystyle=\dfrac{\exp\bigl{(}-{\vec{C}^{\prime}}\rule{0.0pt}{4.30554pt}^{2}/(2\,{c_{q}}^{2})\bigr{)}\,\exp\bigl{(}-{\vec{C}^{\prime}_{*}}\rule{0.0pt}{4.30554pt}^{2}/(2\,{c_{q}}^{2})\bigr{)}}{\bigl{(}\sqrt{2\pi}{c_{q}}\bigr{)}^{3}\,\exp\bigl{(}-{\vec{C}_{*}}\rule{0.0pt}{4.30554pt}^{2}/(2\,{c_{q}}^{2})\bigr{)}}$
$\displaystyle=\dfrac{\exp\bigl{(}-\vec{c}^{\>2}/(2\,{c_{q}}^{2})\bigr{)}}{\bigl{(}\sqrt{2\pi}{c_{q}}\bigr{)}\rule{0.0pt}{4.30554pt}^{3}}$
(17)

because
${\vec{C}^{\prime}}\rule{0.0pt}{4.30554pt}^{2}+{\vec{C}^{\prime}_{*}}\rule{0.0pt}{4.30554pt}^{2}=\vec{c}^{\>2}+{\vec{C}_{*}}\rule{0.0pt}{4.30554pt}^{2}$.
The induction hypothesis applies in this case also.

Because the estimation tree is finite (see the previous point), the induction
hypothesis applies from the leaves of the tree to any node of the tree, root
included. ∎

### Numerical experiments

were carried out in which algorithm 2 was operated to calculate $f$ at several
points of the phase spacetime. The physical situation considered is the BKW
situation expressed in Eq. (7), _ie_ the initial distribution is at maximal
disequilibrium:

$f(\vec{c};0)=\dfrac{5\,\vec{c}^{\>2}}{9\,{c_{q}}^{2}}\,\biggl{(}\sqrt{\dfrac{6\pi}{5}}\,{c_{q}}\biggr{)}^{\\!\\!-3}\exp\biggl{(}-\vec{c}^{\>2}\bigg{/}\biggl{(}\dfrac{6}{5}\,{c_{q}}^{2}\biggr{)}\biggr{)}$
(18)

For each point where $f$ was calculated, e4 realizations were used. The
results are displayed in Figs. 3 and 4.

A first feature of the performance of algorithm 2 is that the variance of its
result increases with the time elapsed since the initial condition. This is
logical. When the elapsed time increases, on average the estimation trees are
padded. Because the relation between the result and the equilibrium value is
exactly the product of these relations at every leaf of the estimation tree
(this can be understood through the 2nd proof just above), more leaves means
more variance in the final result. Otherwise, the modelling process we applied
to the mean recursion ${\mathit{rec}}(\vec{c};t)$ through Eqs. (12) to (14)
can be applied to the estimator variance – this will be the subject of a
future paper – , which does indeed show that this variance increases with the
time $t$.

(a) The threshold is set as $2\times$ the RMS speed.

(b) The threshold is set as $6\times$ the RMS speed.

Figure 5: Results obtained by algorithm 3 combined with the sampling laws
listed in Eqs. (10) and (19), and applied to the BKW mode described in Eq.
(7). Two constant speeds are considered: along the time $t$, we follow how
many particles have a speed above these thresholds. The graphs at the top give
a graphical comparison between the results (given with confidence intervals of
1 standard deviation) and the expected values. The graphs at the bottom show,
with the same results, the relative standard deviation (circles) (given by the
calculations) beside the error actually made (triangles). Each displayed point
was obtained through running e4 realizations of algorithm 3.

(a) Results obtained after 1 mean collision time

(b) Results obtained after 3 mean collision times

Figure 6: Results obtained by algorithm 3 combined with the sampling laws
listed in Eqs. (10) and (19), and applied to the BKW mode described in Eq.
(7). We calculate the distribution of particles along the speed $\|\vec{c}\|$,
_ie_ how many particles have a speed above any given threshold $c_{0}$; only
two dates are considered. The graphs at the top give a graphical comparison
between the results (given with confidence intervals of 1 standard deviation)
and the expected values. The graphs at the bottom show, with the same results,
the relative standard deviation (circles) (given by the calculations) beside
the error actually made (triangles). Each displayed point was obtained through
running e4 realizations of algorithm 3.

Input: A speed $c_{0}$ and a time $t$

Output: A point estimate of
${\mathit{Fra\mspace{-2.5mu}c}}\,(\|\vec{c}\|\geqslant c_{0};t)$, the fraction
of particles of which speed exceeds $c_{0}$ at time $t$

Sample $C\\!_{f}$: $c_{f}$ is obtained;

// $c_{f}\geqslant c_{0}$

Sample $\vec{U}\\!_{f}$: $\vec{u}_{f}$ is obtained;

Estimate $f(c_{f}\vec{u}_{f};t)$ using algorithm 1:
$\widetilde{f}(c_{f}\vec{u}_{f};t)$ is obtained;

return
_$\dfrac{c_{f}\rule{0.0pt}{4.30554pt}^{2}\,\widetilde{f}(c_{f}\vec{u}_{f};t)}{{p_{C\\!_{f}}}(c_{f})\,{p_{\vec{U}\\!_{f}}}(\vec{u}_{f})}$_
;

Algorithm 3 Algorithm for estimating fractions of particles with high kinetic
energy, valid in the conditions of the BKW mode described in section II
(uniform and unit density, Maxwell molecules…)

A more interesting behaviour of algorithm 2 is that the variance of its result
depends very little on the position in the phase space. In particular, it is
nearly independent of the rarefaction: in Fig. 4, we see that when $f$ is
divided by e6 while placing in the velocity space, the relative error
increases only by a factor of 5.

In order to verify this last feature in more depth, we decided to compute
high-energy fractions of the gas. The quantity we calculated was
${\mathit{Fra\mspace{-2.5mu}c}}\,(\|\vec{c}\|\geqslant c_{0};t)$, the fraction
of molecules of which the total speed exceeds a given threshold $c_{0}$ at
time $t$. It is a direct integral of $f$ on ${E_{c}}$: we computed it using
algorithm 2, with the additional initial sampling of a RV
$\vec{C}\\!_{f}=C\\!_{f}\vec{U}\\!_{f}$. Algorithm 3 explains this; the added
sampling law is:

$\left\\{\begin{array}[]{@{\
}l@{\quad}l@{}}{p_{C\\!_{f}}}(c_{f})=\dfrac{c_{0}\Bigl{(}\bigl{(}\tfrac{c_{0}}{{c_{q}}}\bigr{)}^{\\!2}-2\Bigr{)}}{\Bigl{(}\Bigl{(}\bigl{(}\tfrac{c_{0}}{{c_{q}}}\bigr{)}^{\\!2}-2\Bigr{)}(c_{f}-c_{0})+c_{0}\Bigr{)}\big{.}^{\\!2}}&{\text{if
}}c_{0}>2{c_{q}}\\\ \\\\[-4.30554pt]
{p_{C\\!_{f}}}(c_{f})=\dfrac{2{c_{q}}}{(2{c_{q}}+c_{f}-c_{0})^{2}}&{\text{if
}}c_{0}\leqslant 2{c_{q}}\\\ \\\\[-4.30554pt]
{p_{\vec{U}\\!_{f}}}(\vec{u}_{f})=(4\pi)^{-1}&\end{array}\right.$ (19)

The results of algorithms 3 and 2, with sampling laws listed in Eqs. (19) and
(10), are displayed in Figs. 5 and 6. There seems to be no problem in
quantifying tiny parts of the gas, especially those located at high kinetic
energy, even fractions as small as one billionth of the total.

##  III Toy model, with null collisions

### In the previous section,

we used the BKW mode as an example of gas kinetics. This theoretical case is
built on a very particular combination: the uniform density and the Maxwell
collision model together make the collision frequency constant, in the whole
phase space. Thus the extinction in the Boltzmann dynamics becomes
artificially linear. We have written $f(\vec{c};t)$ as a sum of contributions
of the source term due to inward collisions, weighted with the Beer extinction
law, in Eq. (9). This is in fact always feasible, but there the exponential of
the Beer law uses the time only multiplied with a constant coefficient (which
can be set to unity thanks to proper non-dimensionalizing, as was done in this
case), which is very convenient for sampling. In a more general case, the
collision frequency depends on $f$, which introduces through the exponential
term in Eq. (9) an additional non-linearity to the one contained in the source
term.

We show hereafter how to deal with the non-linearity of the extinction, in the
Boltzmann equation. This is done on a very simplified case, built as the
extinction of the Boltzmann dynamics isolated without source term or phase
space. It is only the quadratic ordinary differential system, here written in
$f(t)$:

$\left\\{\begin{array}[]{@{\ }l}f^{\prime}(t)=-\alpha\,f(t)^{2}\\\
f(0)=f_{0}\end{array}\right.$ (20)

where $\alpha$ and $f_{0}$ are arbitrary positive constants. It admits on
$\mathbb{R}^{+}$ the single solution:

$f(t)=\dfrac{f_{0}}{\alpha\,f_{0}\,t+1}$ (21)

### A way to solve

Eq. (20) numerically, with the Monte-Carlo method, would be to use a linear
solution of Eq. (20), as we did in the previous section in Eq. (9). Indeed
here:

$f^{\prime}(t)=-\bigl{(}\alpha\,f(t)\bigr{)}\times f(t)$ (22)

and so

$\forall t\geqslant
0,f(t)=\int_{-\infty}^{t}{\mathrm{d}}t^{\prime}\,\alpha\check{f}(t^{\prime})\exp\biggl{(}-\alpha\int_{t^{\prime}}^{t}{\mathrm{d}}t^{\prime\prime}\,\check{f}(t^{\prime\prime})\biggr{)}\times\rule{0.0pt}{4.30554pt}\\\
\Bigl{(}\operatorname{H}(t^{\prime})\times
0+\operatorname{H}(-t^{\prime})\,f_{0}\Bigr{)}$ (23)

where

$\check{f}(t)=\operatorname{H}(t)\,f(t)+\operatorname{H}(-t)\,f_{0}$ (24)

Combined with the use of independent and identically distributed RVs applied
to a Taylor series of the exponential function, as proposed in section I and
detailed in Dauchet (2012), Eq. (23) leads to a recursive integral writing of
$f(t)$, which can be turned into a Monte-Carlo algorithm. However, this is not
what we are going to use.

### The technique we rely on

is called, among other names, a Null Collision Algorithm (NCA) Boeuf and
Marode (1982). On the examples in this paper, it will remove the need for a
Taylor series. In fact, NCAs are formally equivalent to the use of a Taylor
series of the exponential extinction law, which we mentioned just above; this
has been explained for example in Longo (2002). But the way the NCAs are built
relates to a physical interpretation, and their convergence is very simple to
obtain, which is why we have chosen this approach.

The principle of a NCA is to introduce an arbitrary collision frequency, here
noted $\hat{\nu}$. We thus transform the Eq. (20) in:

$f^{\prime}(t)=-\hat{\nu}f(t)+\Bigl{(}\hat{\nu}f(t)-\alpha f(t)^{2}\Bigr{)}$
(25)

Considering $-\hat{\nu}f(t)$ as an extinction and
$\bigl{(}\hat{\nu}f(t)-\alpha f(t)^{2}\bigr{)}$ as a source term, we now
express:

$\forall t\geqslant
0,f(t)=\int_{-\infty}^{t}{\mathrm{d}}t^{\prime}\,\hat{\nu}\exp\bigl{(}-\hat{\nu}(t-t^{\prime})\bigr{)}\times\rule{0.0pt}{4.30554pt}\\\
\biggl{(}\operatorname{H}(t^{\prime})\biggl{(}f(t)-\dfrac{\alpha\,f(t)^{2}}{\hat{\nu}}\biggr{)}+\operatorname{H}(-t^{\prime})\,f_{0}\biggr{)}$
(26)

This is where the name “Null Collision Algorithm” comes from. The extinction
frequency is set to an arbitrary value $\hat{\nu}$. This is compensated by
adding a source term of particles having encountered null collisions –
collisions without effect. Thus the transport problem remains the same, only
the algorithmic is modified.

In the following we call $\hat{\nu}$ the “raised collision frequency”. This
name arises from the fact that in early implementations of NCAs Skullerud
(1968); Andreucci (1985); Rees (1968), $\hat{\nu}$ was actually having to
exceed the real extinction frequency everywhere for the algorithms to work. It
has been demonstrated recently that this constraint is generally not mandatory
Galtier  _et al._ (2013) (but $\hat{\nu}$ has, however, to be positive);
nevertheless, complying with it improves the convergence in most cases.

(a) Results obtained with $\hat{\nu}=1\,\alpha\,f_{0}$

(b) Results obtained with $\hat{\nu}=\num{2.5}\,\alpha\,f_{0}$

Figure 7: Results obtained by algorithm 4, solving the Eq. (20). $f(t)$ was
calculated at different $t$, using two different values of $\hat{\nu}$ in
columns (a) and (b). The four bottom graphs compare the predictions in Eqs.
(28) and (29) with what algorithm 4 has done. Each point was obtained through
running e4 realizations of algorithm 4.

### By setting a probability density

for the integration variable $t^{\prime}$ in Eq. (26) at every $t\geqslant 0$,
we obtain a recursive Monte-Carlo algorithm for evaluating $f(t)$. We choose
this density arbitrary as the exponential term in Eq. (26):

${p_{T^{\prime}}}:\left\\{\begin{array}[]{@{\
}l}(-\infty;t]\to\mathbb{R}^{+*}\\\
t^{\prime}\mapsto\hat{\nu}\exp\bigl{(}-\hat{\nu}(t-t^{\prime})\bigr{)}\end{array}\right.$
(27)

for the sake of simplicity, and this choice will appear reasonable through the
section. The result is algorithm 4.

Input: A time $t\geqslant 0$

Output: A point estimate of $f(t)$

Sample $T^{\prime}$ following an exponential law to the left of $t$, with
constant $\hat{\nu}$: $t^{\prime}$ is obtained;

if  _$t^{\prime}\leqslant 0$_ then return _$f_{0}$_ ;

else

Estimate $f(t^{\prime})$ using this algorithm: $\widetilde{f}_{1}(t^{\prime})$
is obtained;

Estimate at new $f(t^{\prime})$ using this algorithm:
$\widetilde{f}_{2}(t^{\prime})$ is obtained;

// independently of $\widetilde{f}_{1}(t^{\prime})$

return
_$\widetilde{f}_{1}(t^{\prime})\biggl{(}1-\dfrac{\alpha\,\widetilde{f}_{2}(t^{\prime})}{\hat{\nu}}\biggr{)}$_
;

Algorithm 4 Algorithm for estimating $f(t)$, as described in Eq. (20)

In the same way as shown in the previous section, we can prove that the
expected recursion of algorithm 4, starting at any arbitrary time, is finite.
Defining the recursion ${\mathit{rec}}$ as the total number of $T^{\prime}$
sampling – as in the previous section – we obtain a very similar result:

${\mathit{rec}}(t)=2\exp(\hat{\nu}t)-1$ (28)

The same modelling approach can be applied to the second moment of
$\widetilde{F}(t)$ the estimator of $f(t)$, in order to obtain its variance.
The example given in this section is simple enough to derive a symbolic
expression of the variance. It states that:

$\operatorname{Var}\bigl{(}\widetilde{F}(t)\bigr{)}=\dfrac{\alpha\,f_{0}\,t}{\dfrac{\hat{\nu}}{\alpha\,f_{0}}+\biggl{(}\dfrac{\hat{\nu}}{\alpha\,f_{0}}-1\biggr{)}\,\alpha\,f_{0}\,t}\,f(t)^{2}$
(29a) if
$\hat{\nu}\geqslant\alpha\,f_{0}\quad{\text{or}}\quad\biggl{(}\frac{\alpha\,f_{0}}{\hat{\nu}}-1\biggr{)}\,\alpha\,f_{0}\,t<1$
(29b) , otherwise $\operatorname{Var}\bigl{(}\widetilde{F}(t)\bigr{)}$ is
infinite.

Interestingly, the behaviour of $\widetilde{F}$ varies qualitatively with the
quality of $\hat{\nu}$ as a global bound of the extinction frequency (of which
the maximum is $\alpha\,f_{0}$):

  * •

If $\hat{\nu}>\alpha\,f_{0}$, the relative variance of $\widetilde{F}(t)$,
defined as
$\operatorname{Var}\bigl{(}\widetilde{F}(t)\bigr{)}\big{/}f(t)^{2}$, tends to
a finite value when $t\to+\infty$. So the absolute variance of
$\widetilde{F}(t)$ tends toward 0 while $t\to+\infty$. $\widetilde{F}(t)$ is a
usable estimator of $f(t)$ for every $t$.

  * •

If $\hat{\nu}=\alpha\,f_{0}$, the relative variance of $\widetilde{F}(t)$
grows linearly to infinity with $t$; the absolute variance of
$\widetilde{F}(t)$ still tends towards 0 while $t\to+\infty$.
$\widetilde{F}(t)$ is thus still a reasonable estimator of $f(t)$ regardless
of $t$.

  * •

If $\hat{\nu}<\alpha\,f_{0}$, the variance of $\widetilde{F}(t)$ (absolute or
relative) reaches infinity for a finite value of the time $t$, and for all
subsequent times. In these conditions, obtaining the convergence of a Monte-
Carlo procedure estimating $f(t)$ can be very difficult.

### Algorithm 4 was actually run

with one aim being to confirm what is affirmed in Eqs. (28) and (29). Some
results are given in Fig. 7.

The indications stated above about the recursion and the variance of algorithm
4 are confirmed. Noticeably, we do not have that sharp increase in the
estimator variance with simulated time, which is observed in algorithm 2 aimed
at solving the BKW mode.

Furthermore, in Fig. 7, as in Eqs. (28) and (29), it can be understood that a
compromise needs to be found in the parameter $\hat{\nu}$: between a low
variance associated with a high recursion, and a low recursion associated with
a high variance. This kind of compromise is in fact present when designing
recursive MC algorithms solving any problem presented in this paper. For
example, with oversampling strategies we have been able to increase
simulatable physical times in the BKW mode of the previous section by several
mean free flight times. This will be the subject of a future work.

##  IV Example: harmonic trap

### In this section we solve

the Boltzmann equation with the MCM, in a second academic case where a
symbolic solution is available. In this case, the gas is no longer uniform,
and forms a cloud around the origin which swells and contracts periodically.
Compared to the BKW mode, we have now to account for the ballistic transport
and for a variable collision frequency.

The physical case under consideration here is nevertheless very particular,
because it belongs to the kernel of the collision operator – _ie_ collisions
have no influence –, although it is not the barometric equilibrium. This
possibility has been known since the works of Boltzmann himself Boltzmann
(1909), and has been revisited recently in the field of cold atom gas
manipulation Guéry-Odelin  _et al._ (2014).

We chose our precise test case in the collision operator kernel, built as
simply as possible:

  * •

The molecules are subjected to a constant and purely harmonic force calling
them back to the origin. This force subjects the particles to an acceleration:

$\vec{a}\equiv\vec{a}(\vec{r})=-\omega^{2}\,\vec{r}$ (30)

where $\vec{r}$ is the position.

  * •

The gas has no global kinetic moment around the origin.

Knowing this, the equilibrium distributions constitute a set parametrized by
two quantities, like the total amount of matter $n$ and the thermal RMS speed
on each axis ${c_{q}}$. In this case of elastic confinement, the collision
operator kernel extends to breathing modes oscillating at twice the trap
frequency, with two additional parameters, such as $\Delta{c_{q}}\\!^{2}$ the
amplitude of the thermal energy oscillations and ${\phi_{0}}$ a phase at the
origin of time. In these breathing modes, the distribution of molecules is
always and everywhere Maxwellian, with density, peculiar speed, and thermal
energy:

$\displaystyle\eta_{\text{cok}}(\vec{r};t)=n\,{p_{\mathcal{N}\left(\vec{0};\tfrac{(1+\epsilon\sin\phi(t))\,{c_{q,\,{\text{eq}}}}^{2}}{\omega^{2}}\,\bar{\bar{I}}\right)}}(\vec{r})$
(31a)
$\displaystyle\vec{v}_{\text{cok}}(\vec{r};t)=\dfrac{\epsilon\cos\phi(t)}{1+\epsilon\sin\phi(t)}\,\omega\vec{r}$
(31b)
$\displaystyle{c_{q,\,{\text{cok}}}}(t)^{2}=\dfrac{(1-\epsilon^{2})\,{c_{q,\,{\text{eq}}}}^{2}}{1+\epsilon\sin\phi(t)}$
(31c)

where $\phi(t)$ stands for $2\omega t+{\phi_{0}}$, ${c_{q,\,{\text{eq}}}}$ is
the thermal speed in the equilibrium distribution with same total mass and
same total mechanical energy, $\epsilon$ equals
$\frac{\Delta{c_{q}}\\!^{2}}{{c_{q,\,{\text{eq}}}}^{2}}$, and
${p_{\mathcal{N}(\vec{0};\overline{\overline{V}})}}$ denotes the density of
the centred multidimensional normal probability law with covariance matrix
$\overline{\overline{V}}$.The resulting physical example is presented in Fig.
8.

(a) Density of matter along the $(Ox)$ axis of the ordinary space, as a
function of the time. The physical situation has spherical symmetry. The
dashed lines are density level sets, spaced with a factor 2 or a factor 100.

(b) Distibution of the total mechanical energy of the gas, as a function of
the time.

Figure 8: Presentation of the physical situation of the harmonic trap
described in Eq. (31). The parameters $\epsilon=\frac{9}{41}$ and
${\phi_{0}}=0$ are chosen, as in section V.

Because we do not want to have only ballistic dynamics, we need to choose a
collision model. For the sake of simplicity, we chose the isotropic Maxwell
collision model: the differential collision cross-section
${\sigma_{F}}=\kappa/(4\pi g)$ where $g$ is the relative velocity of particles
and $\kappa$ is a constant.

### Introducing null-collisions,

the Boltzmann equation can be written in our case:

$\partial_{t}f(\vec{r};\vec{c};t)+\vec{c}\cdot\vec{\nabla}_{\\!\vec{r}}\,f(\vec{r};\vec{c};t)-\omega^{2}\,\vec{r}\cdot\vec{\nabla}_{\\!\vec{c}}\,f(\vec{r};\vec{c};t)=\\\\[8.61108pt]
-\hat{\nu}(\vec{r};\vec{c};t)f(\vec{r};\vec{c};t)\\\\[4.30554pt]
\rule{0.0pt}{4.30554pt}+\biggl{(}\hat{\nu}(\vec{r};\vec{c};t)-\kappa\int_{E_{c}}{\mathrm{d}}{\vec{c}_{*}}\,f(\vec{r};{\vec{c}_{*}};t)\biggr{)}f(\vec{r};\vec{c};t)\\\
\rule{0.0pt}{4.30554pt}+\kappa\int_{E_{c}}{\mathrm{d}}{\vec{c}_{*}}\int_{E_{u}}\dfrac{{\mathrm{d}}\vec{u}}{4\pi}\,f(\vec{r};{\vec{c}^{\mspace{3.5mu}\prime}};t)f(\vec{r};{\vec{c}^{\mspace{3.5mu}\prime}_{*}};t)$
(32)

where ${\vec{c}^{\mspace{3.5mu}\prime}}$ and
${\vec{c}^{\mspace{3.5mu}\prime}_{*}}$ are given in Eq. (6), and $\hat{\nu}$
is an arbitrary raised collision frequency. Considering that the distribution
is known at $t=0$ (as in section II), and using Liouville’s theorem, Eq. (32)
can be turned into a Fredholm counterpart:

$f(\vec{r};\vec{c};t)=\int_{-\infty}^{t}{\mathrm{d}}t^{\prime}\,\hat{\nu}\bigl{(}{\vec{r}_{b}}(t^{\prime});{\vec{c}_{b}}(t^{\prime});t^{\prime}\bigr{)}\times\rule{0.0pt}{4.30554pt}\\\
\exp\biggl{(}-\int_{t^{\prime}}^{t}{\mathrm{d}}t^{\prime\prime}\,\hat{\nu}\bigl{(}{\vec{r}_{b}}(t^{\prime\prime});{\vec{c}_{b}}(t^{\prime\prime});t^{\prime\prime}\bigr{)}\biggr{)}\times\rule{0.0pt}{4.30554pt}\\\
\Biggl{(}\operatorname{H}(-t^{\prime})f({\vec{r}_{b}}(0);{\vec{c}_{b}}(0);0)+\operatorname{H}(t^{\prime})\times\rule{0.0pt}{4.30554pt}\\\
\biggl{(}\biggl{(}1-\dfrac{{\nu_{t}}({\vec{r}_{b}}(t^{\prime});{\vec{c}_{b}}(t^{\prime});t^{\prime})}{\hat{\nu}({\vec{r}_{b}}(t^{\prime});{\vec{c}_{b}}(t^{\prime});t^{\prime})}\biggr{)}f\bigl{(}{\vec{r}_{b}}(t^{\prime});{\vec{c}_{b}}(t^{\prime});t^{\prime}\bigr{)}+\rule{0.0pt}{4.30554pt}\\\
\dfrac{{s_{\mathit{is}}}({\vec{r}_{b}}(t^{\prime});{\vec{c}_{b}}(t^{\prime});t^{\prime})}{\hat{\nu}({\vec{r}_{b}}(t^{\prime});{\vec{c}_{b}}(t^{\prime});t^{\prime})}\biggr{)}\Biggr{)}$
(33a) with ${\vec{r}_{b}}$ and ${\vec{c}_{b}}$ describing the ballistic
trajectory flowing past $(\vec{r};\vec{c})$ at $t$: $\left\\{\begin{aligned}
\partial_{t^{\prime}}{\vec{r}_{b}}(t^{\prime})&={\vec{c}_{b}}(t^{\prime})\\\
\partial_{t^{\prime}}{\vec{c}_{b}}(t^{\prime})&=-\omega^{2}{\vec{r}_{b}}(t^{\prime})\end{aligned}\right.\quad{\text{with
the IC}}\ \left\\{\begin{aligned} {\vec{r}_{b}}(t)&=\vec{r}\\\
{\vec{c}_{b}}(t)&=\vec{c}\end{aligned}\right.$ (33b) and ${\nu_{t}}$ and
${s_{\mathit{is}}}$ are the real extinction and source term:
$\displaystyle{\nu_{t}}(\vec{r};\vec{c};t)=\kappa\int_{E_{c}}{\mathrm{d}}{\vec{c}_{*}}\,f(\vec{r};{\vec{c}_{*}};t)$
(33c)
$\displaystyle{s_{\mathit{is}}}(\vec{r};\vec{c};t)=\kappa\int_{E_{c}}{\mathrm{d}}{\vec{c}_{*}}\int_{E_{u}}\dfrac{{\mathrm{d}}{\vec{u}^{\mspace{1.5mu}\prime}}}{4\pi}\,f(\vec{r};{\vec{c}^{\mspace{3.5mu}\prime}};t)f(\vec{r};{\vec{c}^{\mspace{3.5mu}\prime}_{*}};t)$
(33d)

The raised extinction frequency $\hat{\nu}$ can be used as a sampling guide
for $t^{\prime}$. Then, by choosing probability densities in order to sample
${\vec{c}_{*}}$ and ${\vec{u}^{\mspace{1.5mu}\prime}}$, Eq. (33) can be
converted into Monte-Carlo algorithm 5, intended to estimate
$f(\vec{r};\vec{c};t)$ with
$(\vec{r};\vec{c};t)\in{E_{r}}\times{E_{c}}\times\mathbb{R}^{+}$ – where
${E_{r}}$ denotes the space of positions.

Input: A point $(\vec{r};\vec{c};t)$ in
${E_{r}}\times{E_{c}}\times\mathbb{R}^{+}$

Output: A point estimate of $f(\vec{r};\vec{c};t)$

Sample $T^{\prime}$ in $(-\infty;t]$, with the density
${p_{T^{\prime}}}:t^{\prime}\mapsto\hat{\nu}\bigl{(}{\vec{r}_{b}}(t^{\prime});{\vec{c}_{b}}(t^{\prime});t^{\prime}\bigr{)}\exp\bigl{(}-\int_{t^{\prime}}^{t}{\mathrm{d}}t^{\prime\prime}\,\hat{\nu}\bigl{(}{\vec{r}_{b}}(t^{\prime\prime});{\vec{c}_{b}}(t^{\prime\prime});t^{\prime\prime}\bigr{)}\bigr{)}$:
$t^{\prime}$ is obtained;

if  _$t^{\prime}\leqslant 0$_ then

return _$f({\vec{r}_{b}}(0);{\vec{c}_{b}}(0);0)$_ ;

else

Sample ${\vec{C}_{*}}$: ${\vec{c}_{*}}$ is obtained;

Sample ${\vec{U}^{\prime}}$: ${\vec{u}^{\mspace{1.5mu}\prime}}$ is obtained;

${\vec{c}^{\mspace{3.5mu}\prime}}\leftarrow\tfrac{1}{2}\bigl{(}{\vec{c}_{b}}(t^{\prime})+{\vec{c}_{*}}+\|{\vec{c}_{b}}(t^{\prime})-{\vec{c}_{*}}\|\,{\vec{u}^{\mspace{1.5mu}\prime}}\bigr{)}$;

${\vec{c}^{\mspace{3.5mu}\prime}_{*}}\leftarrow\tfrac{1}{2}\bigl{(}{\vec{c}_{b}}(t^{\prime})+{\vec{c}_{*}}-\|{\vec{c}_{b}}(t^{\prime})-{\vec{c}_{*}}\|\,{\vec{u}^{\mspace{1.5mu}\prime}}\bigr{)}$;

Estimate $f({\vec{r}_{b}}(t^{\prime});{\vec{c}_{b}}(t^{\prime});t^{\prime})$
using this algorithm:
$\widetilde{f}_{1}({\vec{r}_{b}}(t^{\prime});{\vec{c}_{b}}(t^{\prime});t^{\prime})$
is obtained;

Estimate $f({\vec{r}_{b}}(t^{\prime});{\vec{c}_{*}};t^{\prime})$ using this
algorithm:
$\widetilde{f}_{2}({\vec{r}_{b}}(t^{\prime});{\vec{c}_{*}};t^{\prime})$ is
obtained;

// independently of
$\widetilde{f}_{1}({\vec{r}_{b}}(t^{\prime});{\vec{c}_{b}}(t^{\prime});t^{\prime})$

Estimate
$f({\vec{r}_{b}}(t^{\prime});{\vec{c}^{\mspace{3.5mu}\prime}};t^{\prime})$
using this algorithm:
$\widetilde{f}_{3}({\vec{r}_{b}}(t^{\prime});{\vec{c}^{\mspace{3.5mu}\prime}};t^{\prime})$
is obtained;

Estimate
$f({\vec{r}_{b}}(t^{\prime});{\vec{c}^{\mspace{3.5mu}\prime}_{*}};t^{\prime})$
using this algorithm:
$\widetilde{f}_{4}({\vec{r}_{b}}(t^{\prime});{\vec{c}^{\mspace{3.5mu}\prime}_{*}};t^{\prime})$
is obtained;

// independently of
$\widetilde{f}_{3}({\vec{r}_{b}}(t^{\prime});{\vec{c}^{\mspace{3.5mu}\prime}};t^{\prime})$

return
_$\widetilde{f}_{1}({\vec{r}_{b}}(t^{\prime});{\vec{c}_{b}}(t^{\prime});t^{\prime})+\kappa\,\dfrac{\widetilde{f}_{3}({\vec{r}_{b}}(t^{\prime});{\vec{c}^{\mspace{3.5mu}\prime}};t^{\prime})\,\widetilde{f}_{4}({\vec{r}_{b}}(t^{\prime});{\vec{c}^{\mspace{3.5mu}\prime}_{*}};t^{\prime})-\widetilde{f}_{1}({\vec{r}_{b}}(t^{\prime});{\vec{c}_{b}}(t^{\prime});t^{\prime})\,\widetilde{f}_{2}({\vec{r}_{b}}(t^{\prime});{\vec{c}_{*}};t^{\prime})}{4\pi\,{p_{{\vec{C}_{*}}}}({\vec{c}_{*}})\,{p_{{\vec{U}^{\prime}}}}({\vec{u}^{\mspace{1.5mu}\prime}})\,\hat{\nu}\bigl{(}{\vec{r}_{b}}(t^{\prime});{\vec{c}_{b}}(t^{\prime});t^{\prime}\bigr{)}}$_
;

Algorithm 5 Algorithm for estimating $f(\vec{r};\vec{c};t)$, valid in the
conditions of the harmonic trap described in sections IV and V

Algorithm 5 has an amazing feature: if the initial distribution is in the
kernel of the collision operator – as described in Eq. (31) –, then whatever
the choices retained for $\hat{\nu}$, ${\vec{C}_{*}}$, and
${\vec{U}^{\prime}}$, the result of algorithm 5 has a null variance. To our
knowledge, this is the first Monte-Carlo algorithm which exhibits a zero
variance result independently of the sampling choices222One could object that
the sampling of $t^{\prime}$ is imposed, because it is linked to the raised
collision frequency $\hat{\nu}$; but $\hat{\nu}$ is itself arbitrary..

###### Principle of proof (Null variance of algorithm 5).

We call $\widetilde{F}(\vec{r};\vec{c};t)$ an estimator of
$f(\vec{r};\vec{c};t)$, obtained through algorithm 5. We assume that the
initial distribution $f(\vec{r};\vec{c};t=0)$ is everywhere Maxwellian, with
the first moments described in Eq. (31).

An essential lemma for this proof is the following: for any arbitrary choice
of
$(\vec{r};\vec{c};t;t^{\prime})\in{E_{r}}\times{E_{c}}\times\mathbb{R}\times\mathbb{R}$,
a density value at $(\vec{r};\vec{c};t)$ is compatible with a Maxwellian
distribution of which the first moments are described in Eq. (31), if and only
if the same density value at
$({\vec{r}_{b}}(t^{\prime});{\vec{c}_{b}}(t^{\prime});t^{\prime})$ is
compatible with a Maxwellian distribution of which the first moments are
described in Eq. (31) with substitution of $t$ by $t^{\prime}$. This means, by
virtue of Liouville’s theorem, that an oscillating Maxwellian distribution, as
described in Eq. (31), is compatible with the ballistic transport in the given
force field. The proof of this, and its extension to more general expressions
of the collision operator kernel, can be found in Guéry-Odelin  _et al._
(2014) and will not be repeated here.

Now consider the following induction hypothesis, to apply anywhere in an
estimation tree: $\widetilde{F}(\vec{r};\vec{c};t)$ equals with no variance
$f(\vec{r};\vec{c};t)$, which describes a Maxwellian distribution, the first
moments of which are given in Eq. (31).

Starting algorithm 5, there are two possibilities: $T^{\prime}\leqslant 0$ and
we are at a leaf of the estimation tree, or $T^{\prime}>0$ and we are at a
non-leaf node of the tree.

If $T^{\prime}\leqslant 0$ then:

$\widetilde{F}(\vec{r};\vec{c};t)=f({\vec{r}_{b}}(0);{\vec{c}_{b}}(0);0)$ (34)

Because of the initial condition and of the lemma, the induction hypothesis is
valid.

If $T^{\prime}>0$ then:

$\widetilde{F}(\vec{r};\vec{c};t)=\widetilde{F}_{1}({\vec{r}_{b}}(T^{\prime});{\vec{c}_{b}}(T^{\prime});T^{\prime})+\rule{0.0pt}{4.30554pt}\\\\[4.30554pt]
\dfrac{\kappa}{4\pi\,{p_{{\vec{C}_{*}}}}({\vec{C}_{*}})\,{p_{{\vec{U}^{\prime}}}}({\vec{U}^{\prime}})\,\hat{\nu}\bigl{(}{\vec{r}_{b}}(T^{\prime});{\vec{c}_{b}}(T^{\prime});T^{\prime}\bigr{)}}\times\rule{0.0pt}{4.30554pt}\\\\[4.30554pt]
\Bigl{(}\widetilde{F}_{3}({\vec{r}_{b}}(T^{\prime});{\vec{C}^{\prime}};T^{\prime})\,\widetilde{F}_{4}({\vec{r}_{b}}(T^{\prime});{\vec{C}^{\prime}_{*}};T^{\prime})-\rule{0.0pt}{4.30554pt}\\\\[4.30554pt]
\widetilde{F}_{1}({\vec{r}_{b}}(T^{\prime});{\vec{c}_{b}}(T^{\prime});T^{\prime})\,\widetilde{F}_{2}({\vec{r}_{b}}(T^{\prime});{\vec{C}_{*}};T^{\prime})\Bigr{)}$

where ${\vec{C}^{\prime}}$ and ${\vec{C}^{\prime}_{*}}$ are the RVs
corresponding to the values ${\vec{c}^{\mspace{3.5mu}\prime}}$ and
${\vec{c}^{\mspace{3.5mu}\prime}_{*}}$ in algorithm 5. If the induction
hypothesis holds for $\widetilde{F}_{1}$, $\widetilde{F}_{2}$,
$\widetilde{F}_{3}$ and $\widetilde{F}_{4}$, then:

$\widetilde{F}(\vec{r};\vec{c};t)=f({\vec{r}_{b}}(T^{\prime});{\vec{c}_{b}}(T^{\prime});T^{\prime})+\rule{0.0pt}{4.30554pt}\\\\[4.30554pt]
\dfrac{\kappa}{4\pi\,{p_{{\vec{C}_{*}}}}({\vec{C}_{*}})\,{p_{{\vec{U}^{\prime}}}}({\vec{U}^{\prime}})\,\hat{\nu}\bigl{(}{\vec{r}_{b}}(T^{\prime});{\vec{c}_{b}}(T^{\prime});T^{\prime}\bigr{)}}\times\rule{0.0pt}{4.30554pt}\\\\[4.30554pt]
\Bigl{(}f({\vec{r}_{b}}(T^{\prime});{\vec{C}^{\prime}};T^{\prime})\,f({\vec{r}_{b}}(T^{\prime});{\vec{C}^{\prime}_{*}};T^{\prime})-\rule{0.0pt}{4.30554pt}\\\\[4.30554pt]
f({\vec{r}_{b}}(T^{\prime});{\vec{c}_{b}}(T^{\prime});T^{\prime})\,f({\vec{r}_{b}}(T^{\prime});{\vec{C}_{*}};T^{\prime})\Bigr{)}$

Considering now that $f$ at $({\vec{r}_{b}}(T^{\prime});T^{\prime})$ follows a
Maxwellian distribution, which respects the detailed balance, we conclude
that:

$\widetilde{F}(\vec{r};\vec{c};t)=f({\vec{r}_{b}}(T^{\prime});{\vec{c}_{b}}(T^{\prime});T^{\prime})$
(35)

Following the lemma, the induction hypothesis is valid.

As long as the estimation tree is finite (which brings conditions on
$\hat{\nu}$), the induction hypothesis applies from the leaves of the tree to
any node of the tree, root included. ∎

Some numerical experiments (not displayed here) have confirmed this result.

##  V Example: harmonic trap without local equilibrium

### There is not a lot to illustrate

in the previous section: in the physical case described, the proposed
algorithm has zero variance. There are of course a lot of possible
modifications of algorithm 5 which break with this state of affairs. We have
however chosen to modify the physical case under study, as it seemed pointless
to present a new algorithm, with the only noticeable characteristic being less
efficient than the previous one.

The physical case of this section is a mix between the case of the previous
section and the BKW mode. We suppose the same external force field, and the
same first moments at the initial time, as in section IV. But now, at the
initial condition, at each point of space the distribution of speeds is the
BKW mode distribution at maximal disequilibrium. This means that at every
point of space the initial distribution of speeds is the one described in Eq.
(18), translated and scaled to accord with the density, peculiar speed, and
temperature prescribed at initial time $t=0$ by Eq. (31).

After the initial condition, the system will evolve in a different way to that
explained in section IV. The differences appear even in macroscopic
quantities. As far as the authors know, there is today no symbolic description
available of the evolution of our system. Nevertheless, two limiting
situations are easily described:

  * •

If the cross-sections are null (_ie_ there are no collisions), the evolution
of the system is made only by ballistic transport. This is symbolically
calculable, in a harmonic force field.

  * •

If the cross-sections are not null, the final state of the gas is exactly the
oscillating state described in the previous section IV, in Eq. (31). This can
be demonstrated using the collisional invariants existing in our system.
Considering these, the only state belonging to the collision operator kernel
that the system can reach is the one described in Eq. (31).

### In order to build simulations,

the last remaining parameters must be set.

Input: A point $(\vec{r};\vec{c};t)$ in
${E_{r}}\times{E_{c}}\times\mathbb{R}^{+}$

Output: A point estimate of $f(\vec{r};\vec{c};t)$

$t^{\prime}\leftarrow t$;

repeat

Sample ${T^{\diamond}}$ in $(-\infty;t^{\prime}]$ following an exponential law
to the left of $t^{\prime}$, with constant
$\hat{\vphantom{\rule{1.0pt}{6.14584pt}}\smash{\hat{\nu}}}$: ${t^{\diamond}}$
is obtained;

if  _${t^{\diamond}}\leqslant 0$_ then

return _$f({\vec{r}_{b}}(0);{\vec{c}_{b}}(0);0)$_ ;

Calculate $\eta_{\text{cok}}({\vec{r}_{b}}({t^{\diamond}});{t^{\diamond}})$;
${\eta^{\diamond}}$ is obtained;

$\hat{\nu}^{\prime}\leftarrow\frac{1}{e}\bigl{(}\frac{5}{3}\bigr{)}\\!^{\frac{5}{2}}\,\kappa\,{\eta^{\diamond}}$;

$t^{\prime}\leftarrow{t^{\diamond}}$;

Sample $R$ following a uniform standard law; $r$ is obtained;

until  _$r\hat{\vphantom{\rule{1.0pt}{6.14584pt}}\smash{\hat{\nu}}}
<\hat{\nu}^{\prime}$_;

Calculate ${c_{q,\,{\text{cok}}}}(t^{\prime})$; ${c_{q}}^{\prime}$ is
obtained;

Calculate $\vec{v}_{\text{cok}}({\vec{r}_{b}}(t^{\prime});t^{\prime})$;
${\vec{v}\mspace{2.0mu}^{\prime}}$ is obtained;

Sample ${\vec{C}_{*}}$ following a Maxwellian distribution, with peculiar
speed ${\vec{v}\mspace{2.0mu}^{\prime}}$ and RMS speed ${c_{q}}^{\prime}$:
${\vec{c}_{*}}$ is obtained;

Sample ${\vec{U}^{\prime}}$ isotropically: ${\vec{u}^{\mspace{1.5mu}\prime}}$
is obtained;

${\vec{c}^{\mspace{3.5mu}\prime}}\leftarrow\tfrac{1}{2}\bigl{(}{\vec{c}_{b}}(t^{\prime})+{\vec{c}_{*}}+\|{\vec{c}_{b}}(t^{\prime})-{\vec{c}_{*}}\|\,{\vec{u}^{\mspace{1.5mu}\prime}}\bigr{)}$;

${\vec{c}^{\mspace{3.5mu}\prime}_{*}}\leftarrow\tfrac{1}{2}\bigl{(}{\vec{c}_{b}}(t^{\prime})+{\vec{c}_{*}}-\|{\vec{c}_{b}}(t^{\prime})-{\vec{c}_{*}}\|\,{\vec{u}^{\mspace{1.5mu}\prime}}\bigr{)}$;

Estimate $f({\vec{r}_{b}}(t^{\prime});{\vec{c}_{b}}(t^{\prime});t^{\prime})$
using this algorithm:
$\widetilde{f}_{1}({\vec{r}_{b}}(t^{\prime});{\vec{c}_{b}}(t^{\prime});t^{\prime})$
is obtained;

Estimate $f({\vec{r}_{b}}(t^{\prime});{\vec{c}_{*}};t^{\prime})$ using this
algorithm:
$\widetilde{f}_{2}({\vec{r}_{b}}(t^{\prime});{\vec{c}_{*}};t^{\prime})$ is
obtained;

// independently of
$\widetilde{f}_{1}({\vec{r}_{b}}(t^{\prime});{\vec{c}_{b}}(t^{\prime});t^{\prime})$

Estimate
$f({\vec{r}_{b}}(t^{\prime});{\vec{c}^{\mspace{3.5mu}\prime}};t^{\prime})$
using this algorithm:
$\widetilde{f}_{3}({\vec{r}_{b}}(t^{\prime});{\vec{c}^{\mspace{3.5mu}\prime}};t^{\prime})$
is obtained;

Estimate
$f({\vec{r}_{b}}(t^{\prime});{\vec{c}^{\mspace{3.5mu}\prime}_{*}};t^{\prime})$
using this algorithm:
$\widetilde{f}_{4}({\vec{r}_{b}}(t^{\prime});{\vec{c}^{\mspace{3.5mu}\prime}_{*}};t^{\prime})$
is obtained;

// independently of
$\widetilde{f}_{3}({\vec{r}_{b}}(t^{\prime});{\vec{c}^{\mspace{3.5mu}\prime}};t^{\prime})$

return
_$\widetilde{f}_{1}({\vec{r}_{b}}(t^{\prime});{\vec{c}_{b}}(t^{\prime});t^{\prime})+\dfrac{\kappa\,(\sqrt{2\pi}{c_{q}}^{\prime})^{3}}{\hat{\nu}^{\prime}}\exp\biggl{(}\dfrac{({\vec{c}_{*}}-{\vec{v}\mspace{2.0mu}^{\prime}})^{2}}{2\,{c_{q}}^{\prime\,2}}\biggr{)}\,\Bigl{(}\widetilde{f}_{3}({\vec{r}_{b}}(t^{\prime});{\vec{c}^{\mspace{3.5mu}\prime}};t^{\prime})\,\widetilde{f}_{4}({\vec{r}_{b}}(t^{\prime});{\vec{c}^{\mspace{3.5mu}\prime}_{*}};t^{\prime})-\widetilde{f}_{1}({\vec{r}_{b}}(t^{\prime});{\vec{c}_{b}}(t^{\prime});t^{\prime})\,\widetilde{f}_{2}({\vec{r}_{b}}(t^{\prime});{\vec{c}_{*}};t^{\prime})\Bigr{)}$_
;

Algorithm 6 Algorithm for estimating $f(\vec{r};\vec{c};t)$, valid in the
conditions of the harmonic trap described in sections IV and V. The sampling
choices detailed in section V are used.

We chose $\epsilon=\frac{9}{41}$ and ${\phi_{0}}=0$, and this set the initial
condition. The other parameters $\omega$, $n$, and ${c_{q,\,{\text{eq}}}}$ are
irrelevant since they are squeezed during non-dimensionalization. This is the
choice we made while expounding the physical case of section IV in Fig. 8. We
must also specify the collision model: we chose $\kappa$ such as
$\frac{n\,\kappa\,\omega^{2}}{{c_{q,\,{\text{eq}}}}^{3}}=3$.

Finally, we need to specify $\hat{\nu}$, ${\vec{C}_{*}}$, and
${\vec{U}^{\prime}}$. ${\vec{U}^{\prime}}$ will follow an isotropic law, as in
section II. ${\vec{C}_{*}}$ is set to a Maxwellian sampling; the Maxwellian
bell is adjusted according to the final distribution, _ie_ it is centred on
$\vec{v}_{\text{cok}}(\vec{r};t^{\prime})$ and reduced to
${c_{q,\,{\text{cok}}}}(t^{\prime})$ standard deviation given in Eq. (31),
where $t^{\prime}$ is the collision date. As we wanted $\hat{\nu}$ as a
credible bound of the collision frequency, we made the following choice:

$\hat{\nu}\equiv\hat{\nu}(\vec{r};t)=\dfrac{1}{e}\biggl{(}\dfrac{5}{3}\biggr{)}^{\\!\raisebox{-4.52083pt}{$\frac{5}{2}$}}\,\kappa\,\eta_{\text{cok}}(\vec{r};t)$
(36)

The pre-factor $\frac{1}{e}\bigl{(}\frac{5}{3}\bigr{)}\\!^{\frac{5}{2}}$ is
the maximal ratio between the BKW mode initial distribution and the
equilibrium distribution with the same first moments. Indeed:

$\dfrac{1}{e}\biggl{(}\dfrac{5}{3}\biggr{)}^{\\!\raisebox{-4.52083pt}{$\frac{5}{2}$}}=\max_{\vec{c}\in{E_{c}}}\,\dfrac{f_{\text{BKW}}(\vec{c};0)}{(\sqrt{2\pi}{c_{q}})^{-3}\exp\bigl{(}-\vec{c}^{\>2}/(2\,{c_{q}}^{2})\bigr{)}}$
(37)

where $f_{\text{BKW}}(\vec{c};0)$ is given by Eq. (18).

Sampling a collision date $t^{\prime}$ using the raised collision frequency
$\hat{\nu}$ given in Eq. (36) can be tricky: we need to apply the Beer
extinction law to the gas density described in Eq. (31), along an elliptic
ballistic trajectory…To achieve this, we use $\hat{\nu}$ as the true
extinction in an internal Null Collision Algorithm. This NCA uses as a raised
collision frequency the global maximum of $\hat{\nu}$ defined in Eq. (36),
reached in the centre of the gas cloud at its maximum contraction. We denote
it $\hat{\vphantom{\rule{1.0pt}{6.14584pt}}\smash{\hat{\nu}}}$, and it equals:

$\displaystyle\hat{\vphantom{\rule{1.0pt}{6.14584pt}}\smash{\hat{\nu}}}$
$\displaystyle=\dfrac{1}{e}\biggl{(}\dfrac{5}{3}\biggr{)}^{\\!\raisebox{-4.52083pt}{$\frac{5}{2}$}}\,\kappa\,\eta_{\text{cok}}(\vec{0};t)\big{|}_{\sin\phi(t)=-1}$
(38)
$\displaystyle=\dfrac{1}{e}\biggl{(}\dfrac{5}{3}\biggr{)}^{\\!\raisebox{-4.52083pt}{$\frac{5}{2}$}}\,\kappa\,n\,\biggl{(}\dfrac{\omega}{\sqrt{2\pi}\sqrt{1-\epsilon}\,{c_{q,\,{\text{eq}}}}}\biggr{)}^{\\!3}$
(39)

With all these choices, algorithm 5 becomes algorithm 6.

(a) Results obtained at the centre of the phase space

(b) Results obtained at $2\times$ the RMS radius away the centre of the space
of positions

Figure 9: Results obtained by algorithm 6, with the choices presented in
subsection V, and applied to the harmonic trap problem described in section IV
and V. Two probe points in the phase space are considered, where we follow the
distribution function $f$ along the time $t$. These probe points
$(\vec{r};\vec{c})$ comply with the constraint
$\vec{c}=\vec{v}_{\text{cok}}(\vec{r};0)$, such that there $f=0$ at the
initial instant. The graphs on the top give a graphical comparison between the
results (given with confidence intervals of 1 standard deviation), the values
predicted in the final oscillating state (dashed lines), and the values
predicted by ballistic transport of the initial condition (dotted lines). The
graphs of the middle line show, with the same results, the relative standard
deviation given by the calculation. The graphs of the bottom display the mean
recursivity of algorithm 6 computing the points on the graphs above. Each
displayed points has been obtained through running e4 realizations of
algorithm 6.

(a) Results obtained at $\omega t=2\pi$

(b) Results obtained at $\omega t=6\pi$

Figure 10: Results obtained by algorithm 7 combined with the sampling laws
listed in Eq. (40) and section V, and applied to the harmonic trap problem
described in sections IV and V. We calculate the distibution of particles
along the radius $\|\vec{r}\|$, _ie_ how many particles are more distant from
the origin than any given threshold $r_{0}$; only two dates are considered.
The graphs at the top give a graphical comparison between the results (given
with confidence intervals of 1 standard deviation), and the values predicted
in the final oscillating state (dashed lines) (ballistic transport of the
initial condition leads to the same results, in these precise cases). The
graphs in the middle line show, with the same results, the relative standard
deviation given by the calculation. The graphs at the bottom display the mean
recursivity of algorithm 6 called by algorithm 7, when computing the points on
the graphs above. Each displayed point was obtained through running e4
realizations of algorithm 7.

### We then ran algorithm 6

with all the physical and numerical settings stated previously in this
section. The results are displayed in Figs. 9 and 10. What is noted on these
figures as ${\mathit{rec}}$ is the mean recursion of algorithm 6, which means
the mean number of internal calls of algorithm 6 as it is written, when
calling it to estimate $f$ at a particular point $(\vec{r};\vec{c};t)$ of the
phase space-time.

In Fig. 9 the distribution function $f$ is followed at two constant probe
points of the phase space. As in the previous sections, all calculated points
are independent, which is why we can focus on any particular region of the
phase space-time, as we did around $\omega t=3$.

Fig. 10 shows the distribution of the mass along the potential energy (or
equivalently, along the distance about the origin). The fraction of the mass
more distant from the origin than a fixed threshold $r_{0}$ is computed, as a
function of this threshold. This is an integral of the distribution function
$f$ over the phase space. It is calculated using algorithm 6, with an
additional initial sampling of a radius $\vec{r}\\!_{f}$ and a velocity
$\vec{c}_{f}$. The result is algorithm 7. The sampling law of
$(\vec{R}_{f};\vec{C}\\!_{f})$ is similar to the one given in Eq. (19),
indeed:

$\left\\{\begin{array}[]{@{\,}l@{}}\vec{R}_{f}=R_{f}\vec{U}\\!_{f}\quad{\text{with}}\\\\[4.30554pt]
\left\\{\begin{array}[]{@{\,}l@{}}{p_{R_{f}}}(r\\!_{f})=\rule{0.0pt}{4.30554pt}\\\
\quad\left\\{\begin{array}[]{@{\,}l@{\
}l@{}}\dfrac{r_{0}\Bigl{(}\bigl{(}\tfrac{\omega
r_{0}}{{c_{q,\,{\text{eq}}}}}\bigr{)}^{\\!2}-2\Bigr{)}}{\Bigl{(}\Bigl{(}\bigl{(}\tfrac{\omega
r_{0}}{{c_{q,\,{\text{eq}}}}}\bigr{)}^{\\!2}-2\Bigr{)}(r\\!_{f}-r_{0})+r_{0}\Bigr{)}\big{.}^{\\!2}}&{\text{if
}}\omega r_{0}>2{c_{q,\,{\text{eq}}}}\\\ \\\\[-4.30554pt]
\dfrac{2{c_{q,\,{\text{eq}}}}/\omega}{(2{c_{q,\,{\text{eq}}}}/\omega+r\\!_{f}-r_{0})^{2}}&{\text{if
}}\omega r_{0}\leqslant 2{c_{q,\,{\text{eq}}}}\\\ \end{array}\right.\\\
\\\\[-4.30554pt] {p_{\vec{U}\\!_{f}}}(\vec{u}_{f})=(4\pi)^{-1}\\\
\end{array}\right.\\\ \\\\[-4.30554pt] \vec{C}\\!_{f}{\text{ of law
}}\mathcal{N}\Bigl{(}\vec{v}_{\text{cok}}\bigl{(}\vec{R}_{f};t\bigr{)};{c_{q,\,{\text{cok}}}}(t)^{2}\bar{\bar{I}}\Bigr{)}\end{array}\right.$
(40)

Input: A radius $r_{0}$ and a time $t$

Output: A point estimate of
${\mathit{Fra\mspace{-2.5mu}c}}\,(\|\vec{r}\|\geqslant r_{0};t)$, the fraction
of particles more distant than $r_{0}$ from the origin at time $t$

Sample $R_{f}$: $r\\!_{f}$ is obtained;

// $r\\!_{f}\geqslant r_{0}$

Sample $\vec{U}\\!_{f}$: $\vec{u}_{f}$ is obtained;

Sample $\vec{C}\\!_{f}$: $\vec{c}_{f}$ is obtained;

Estimate $f(r\\!_{f}\vec{u}_{f};\vec{c}_{f};t)$ using algorithm 6:
$\widetilde{f}(r\\!_{f}\vec{u}_{f};\vec{c}_{f};t)$ is obtained;

return
_$\dfrac{r_{f}\rule{0.0pt}{4.30554pt}^{2}\,\widetilde{f}(r\\!_{f}\vec{u}_{f};\vec{c}_{f};t)}{{p_{R_{f}}}(r\\!_{f})\,{p_{\vec{U}\\!_{f}}}(\vec{u}_{f})\,{p_{\vec{C}\\!_{f}}}(\vec{c}_{f})}$_
;

Algorithm 7 Algorithm for estimating fractions of particles with high
potential energy, valid in the conditions of the harmonic trap described in
sections IV and V

We observe the same behaviours in algorithm 6 as in algorithm 2. The variance
of the obtained estimator of $f$ increases with the simulated physical time.
It is however practically insensitive to rarefaction; we can evaluate without
difficulty fractions of the mass as small as one millionth.

There is a noticeable decrease in the recursivity of algorithm 6 when the
probe points move away from the origin. This was an expected result, as the
raised collision frequency $\hat{\nu}$ is higher near the origin. Also, we can
appreciate now the usefulness of our somewhat complex raised extinction
frequency field $\hat{\nu}$ set in Eq. (36). If we had instead used a constant
raised collision frequency such as
$\hat{\vphantom{\rule{1.0pt}{6.14584pt}}\smash{\hat{\nu}}}$ defined in Eq.
(38), the mean recursion would have been higher than shown in Figs. 9 and 10.
However, it is already very high there, especially considering that at the
maximum physical time we probed, the system was still very far from its final
state.

##  VI Conclusion and perspectives

### We have numerically solved

the Boltzmann equation of gas kinetics in several test cases, using the Monte-
Carlo method. It is the ordinary Monte-Carlo method, as used in linear
physics, without bias, and where particle paths are computed independently.
The algorithms that we have proposed are based on a second kind Fredholm
writing of the Boltzmann equation, and rely on the following of virtual
particles from their arrivals to their sources. As far as the authors know,
this is the first numerical method in gas kinetics which thus travels along
time, backwards.

This building of particle paths, from the end to the start, brings an
interesting property to the method. No difficulty appears when focusing on low
crowded parts of the phase space. For the same reason we can carry out _probe
calculus_. This ability is widely known in linear transport physics. An
example from radiative transfer is the prediction of the view of a satellite
looking at the Earth: it is commonly known that one can calculate the image by
sampling rays built as incoming to the satellite and directed from the Earth,
and then by evaluating the possible sources of these rays; and that it is
unneccesary to start by computing the whole irradiance field around the Earth.
The paths built from their destination can explore the phase space, and catch
the information strictly necessary to compute the desired result.

### This possibility to use particle paths

built toward sources is usually considered to be linked to a propagative
insight of the phenomenon. This means that what is observed at the probe point
is a sum of the contributions of all possible paths coming from sources. This
insight remains part of our approach, and the non-linearity is taken into
account at the nodes of the henceforth branching paths. So the contributing
paths are now contributing trees. Each tree consists of the path of the
incoming particle, plus the paths of its collision partners, plus the paths of
the collision partners of the aforementioned partners, _etc_ …We have obtained
an arborescent propagation.

This comes at a price. When the mean number of collisions encountered by the
particles inside the simulated space-time increases, the contributing trees
are padded, and the recursivity and estimation variance of our algorithms can
make their calculations unaffordable.

We are convinced that this extension of a propagative and statistical point of
view, performed here to examine gas kinetics, is amenable to aspects of other
non-linear physics problems. Indeed, it has already been applied recently to
some of them Rasulov  _et al._ (2010); Dimov and Gurov (2000); Henry-Labordère
_et al._ (2014).

### Another striking feature

of the computations we have done is the complete absence of any mesh within.
This is again a common property of the Monte-Carlo method in linear transport
physics, but is new in gas kinetics numericals. Here it is pointless to worry
about the good coverage of the space of speeds, about the bias that the mesh
introduces, or about the memory usage of the numerical programme.

If a geometry is present (a possibility not illustrated in this article), it
enters the integral formulation of the Boltzmann equation only through an
intersection calculus, between surfaces and ballistic trajectories. Otherwise,
the boundary conditions in gas kinetics are generally written as integral
expressions of the distribution of the outgoing molecules, directly usable in
our Monte-Carlo method. This leads the authors to state that taking into
account any geometry will bring no difficulty, given that ballistic
trajectories are simple enough – typically, straight lines. If very complex
geometries come into play, we will be able to use the expertise of the image
synthesis community directly De la Torre  _et al._ (2014).

### Finally,

setting aside the theoretical implications this work could have, we have still
obtained a numerical method in gas kinetics that is very complementary to what
is currently available.

It can compute the high energy fractions in the gas very easily, which is
reputed to be a high burden for the DSMC method or for fully discretized
methods in gas kinetics. These fractions can be very important from the
applicative point of view, because some processes happen only here – _eg_
chemical reactions, ionization events, nuclear reactions or degradation
processes.

However, today it works mainly in non-stationary problems, where the initial
condition plays a major role. Going to stationary situations is very
difficult. On the other hand, stationary problems are reputedly easier for the
DSMC or fully discretized methods.

## Acknowledgements

This work has been sponsored by the French government research-program
“Investissements d’avenir” through the ANR program ANR-2012-IS04-0003-03
DEPART, the Laboratories of Excellence ANR-10-LABX-22-01 SOLSTICE and
ANR-10-LABX-16-01 IMobS3, the ATS program ALGUE of the IDEX of Toulouse
ANR-11-IDEX-02 UNITI, by the Occitanie region, and by the LAPLACE laboratory
through a BQR program.

## References

  * Howell (1998) J. R. Howell, Journal of Heat Transfer 120, 547 (1998).
  * Metropolis and Ulam (1949) N. Metropolis and S. Ulam, Journal of the American Statistical Association 44, 335 (1949).
  * Dunn and Shultis (2012) W. L. Dunn and J. K. Shultis,  _Exploring Monte Carlo Methods_ (Elsevier, 2012) p. 391.
  * Dimov (2008) I. T. Dimov,  _Monte Carlo Methods for Applied Scientists_ (2008) p. 308.
  * Karatzas and Shreve (1991) I. Karatzas and S. E. Shreve,  _Acta Applicandae Mathematicae_ , 2nd ed. (Springer-Verlag, 1991) p. 470.
  * Feynman  _et al._ (2010) R. P. Feynman, A. R. Hibbs, and D. F. Styer,  _Quantum mechanics and path integrals_ (Dover Publications, 2010) p. 371.
  * Umrigar  _et al._ (2007) C. J. Umrigar, J. Toulouse, C. Filippi, S. Sorella, and R. G. Hennig, Physical Review Letters 98, 1 (2007), arXiv:0611094 [cond-mat] .
  * Curtiss (1954) J. H. Curtiss, Journal of Mathematics and Physics 32, 209 (1954).
  * Kalos and Whitlock (2008) M. H. Kalos and P. A. Whitlock,  _Monte carlo methods_ (John Wiley & Sons, 2008) p. 203.
  * Ermakov  _et al._ (1976) S. M. Ermakov, V. V. Nekrutkin, A. Y. Proshkin, and A. F. Sizova, Doklady Akademii nauk SSSR 230, 261 (1976).
  * Gurov (1992) T. V. Gurov, Mathematica Balkanica 6, 237 (1992).
  * Dimov and Gurov (2000) I. T. Dimov and T. V. Gurov, Pliska Studia Mathematica Bulgarica 13, 117 (2000).
  * Rasulov  _et al._ (2004) A. Rasulov, A. N. Karaivanova, and M. Mascagni, Monte Carlo Methods and Applications 10, 551 (2004).
  * Rasulov  _et al._ (2010) A. Rasulov, G. Raimova, and M. Mascagni, Mathematics and Computers in Simulation 80, 1118 (2010).
  * Raimova (2016) G. Raimova, Communications in Statistics - Simulation and Computation 45, 2981 (2016).
  * Henry-Labordère  _et al._ (2014) P. Henry-Labordère, X. Tan, and N. Touzi, Stochastic Processes and their Applications 124, 1112 (2014), arXiv:1302.4624 .
  * Dauchet  _et al._ (2016) J. Dauchet, J.-J. Bezian, S. Blanco, C. Caliot, J. Charon, C. Coustet, M. El Hafi, V. Eymet, O. Farges, V. Forest, R. Fournier, M. Galtier, J. Gautrais, A. Khuong, L. Pelissier, B. Piaud, M. Roger, G. Terrée, and S. Weitz, arXiv (2016), arXiv:1610.02684v1 .
  * Dauchet (2012) J. Dauchet,  _Analyse radiative des photobioréacteurs_ , Ph.D. thesis, Université Blaise Pascal - Clermont-Ferrand II (2012).
  * Farges (2014) O. Farges,  _Conception optimale de centrales solaires à concentration : application aux centrales à tour et aux installations « beam down »_ , Ph.D. thesis, Institut National Polytechnique de Toulouse - École Nationale Supérieure des Mines d’Albi-Carmaux (2014).
  * De la Torre  _et al._ (2014) J. De la Torre, G. Baud, J.-J. Bézian, S. Blanco, C. Caliot, J.-F. Cornet, C. Coustet, J. Dauchet, M. El Hafi, V. Eymet, R. Fournier, J. Gautrais, O. Gourmel, D. Joseph, N. Meilhac, A. Pajot, M. Paulin, P. Perez, B. Piaud, M. Roger, J. Y. Rolland, F. Veynandt, and S. Weitz, Solar Energy 103, 653 (2014).
  * Dauchet  _et al._ (2015) J. Dauchet, S. Blanco, J.-F. Cornet, and R. Fournier, Journal of Quantitative Spectroscopy and Radiative Transfer 161, 60 (2015).
  * Charon  _et al._ (2015) J. Charon, S. Blanco, J.-F. Cornet, J. Dauchet, M. El Hafi, R. Fournier, M. K. Abboud, and S. Weitz, Journal of Quantitative Spectroscopy and Radiative Transfer 172, 3 (2015).
  * Galtier  _et al._ (2013) M. Galtier, S. Blanco, C. Caliot, C. Coustet, J. Dauchet, M. El Hafi, V. Eymet, R. Fournier, J. Gautrais, A. Khuong, B. Piaud, and G. Terrée, Journal of Quantitative Spectroscopy and Radiative Transfer 125, 57 (2013).
  * Novák  _et al._ (2014) J. Novák, A. Selle, and W. Jarosz, ACM Transactions on Graphics 33, 11 (2014).
  * Szirmay-Kalos  _et al._ (2017) L. Szirmay-Kalos, I. Georgiev, M. Magdics, B. Molnár, and D. Légrády, Eurographics 36, 9 (2017).
  * Rees (1968) H. D. Rees, Solid State Communications 6, 643 (1968).
  * Woodcock  _et al._ (1965) E. Woodcock, T. Murphy, P. Hemmings, and S. Longworth, in  _Applications of Computing Methods to Reactor Problems_ (1965).
  * Spanier (1966) J. Spanier, SIAM Journal on Applied Mathematics 14, 702 (1966).
  * Coleman (1968) W. A. Coleman, Nuclear Science and Engineering 32, 76 (1968).
  * Skullerud (1968) H. R. Skullerud, Journal of Physics D: Applied Physics 1, 1567 (1968).
  * Longo (2002) S. Longo, Physica A 313, 389 (2002).
  * Bird (1994) G. A. Bird,  _Molecular Gas Dynamics and the Direct Simulation of Gas Flows_ (Clarendon, 1994).
  * Bird (1998) G. A. Bird, Computers & Mathematics with Applications 35, 1 (1998).
  * Oran  _et al._ (1998) E. S. Oran, C. K. Oh, and B. Z. Cybyk, Annual Review of Fluid … 30, 403 (1998).
  * Murrone and Villedieu (2011) A. Murrone and P. Villedieu, AerospaceLab , 1 (2011).
  * Khisamutdinov (2004) A. I. Khisamutdinov, Transport Theory and Statistical Physics 33, 69 (2004).
  * Stefanov (2011) S. K. Stefanov, SIAM Journal on Scientific Computing 33, 677 (2011).
  * Homolle and Hadjiconstantinou (2007) T. M. M. Homolle and N. G. Hadjiconstantinou, Journal of Computational Physics 226, 2341 (2007).
  * Degond  _et al._ (2010) P. Degond, G. Dimarco, and L. Pareschi, International Journal for Numerical Methods in Fluids 67, 189 (2010).
  * Hauge and Præstgaard (1981) E. H. Hauge and E. Præstgaard, Journal of Statistical Physics 24, 21 (1981).
  * Guéry-Odelin  _et al._ (2014) D. Guéry-Odelin, J. G. Muga, M. J. Ruiz-Montero, and E. Trizac, Physical Review Letters 112, 1 (2014).
  * Bobylev (1976) A. V. Bobylev, Doklady Akademii nauk SSSR 231, 571 (1976).
  * Krook and Wu (1976) M. Krook and T. T. Wu, Physical Review Letters 36, 1107 (1976).
  * Note (1) The BKW mode can be generalized to any dimensionality Ernst (1979).
  * Moradi  _et al._ (2014) M. Moradi, Z. Ayati, and M. A. Mirzazadeh, International Journal of Applied Mathematics and Computation 5, 35 (2014).
  * Boeuf and Marode (1982) J.-P. Boeuf and E. Marode, Journal of Physics D: Applied Physics 15, 2169 (1982).
  * Andreucci (1985) N. Andreucci, in  _Lecture Notes in Physics_ (1985) p. 291.
  * Boltzmann (1909) L. Boltzmann,  _Wissenschaftliche Abhandlungen_ , edited by F. Hasenöhrl (Barth, J. A., Leipzig, 1909).
  * Note (2) One could object that the sampling of $t^{\prime}$ is imposed, because it is linked to the raised collision frequency $\cc@accent{"705E}{\nu}$; but $\cc@accent{"705E}{\nu}$ is itself arbitrary.
  * Ernst (1979) M. H. Ernst, Physics Letters 69, 390 (1979).