# Predicting the statistical error of analog particle tracing Monte Carlo

V. Maes  Department of Computer Science, KU Leuven, Belgium  I. Bossuyt
Department of Computer Science, KU Leuven, Belgium  H. Vandecasteele
Department of Computer Science, KU Leuven, Belgium  W. Dekeyser  Department of
Mechanical Engineering, KU Leuven, Belgium  J. Koellermeier  Bernoulli
Institute for Mathematics, Computer Science and Artificial Intelligence,
University of Groningen, Nijenborgh 4, NL-9747 AG Groningen, Netherlands
Groningen Cognitive Systems and Materials Center, University of Groningen,
Nijenborgh 4, NL-9747 AG Groningen, Netherlands  M. Baelmans  Department of
Mechanical Engineering, KU Leuven, Belgium  G. Samaey  Department of Computer
Science, KU Leuven, Belgium

## Abstract

Large particle systems are often described by high-dimensional (linear)
kinetic equations that are simulated using Monte Carlo methods for which the
asymptotic convergence rate is independent of the dimensionality. Even though
the asymptotic convergence rate is known, predicting the actual value of the
statistical error remains a challenging problem. In this paper, we show how
the statistical error of an analog particle tracing Monte Carlo method can be
calculated (expensive) and predicted a priori (cheap) when estimating
quantities of interest (QoI) on a histogram. We consider two types of QoI
estimators: point estimators for which each particle provides one independent
contribution to the QoI estimates, and analog estimators for which each
particle provides multiple correlated contributions to the QoI estimates. The
developed statistical error predictors can be applied to other QoI estimators
and nonanalog simulation routines as well. The error analysis is based on
interpreting the number of particle visits to a histogram bin as the result of
a (correlated) binomial experiment. The resulting expressions can be used to
optimize (non)analog particle tracing Monte Carlo methods and hybrid
simulation methods involving a Monte Carlo component, as well as to select an
optimal particle tracing Monte Carlo method from several available options.
Additionally, the cheap statistical error predictors can be used to determine
a priori the number of particles $N$ that is needed to reach a desired
accuracy. We illustrate the theory using a linear kinetic equation describing
neutral particles in the plasma edge of a fusion device and show numerical
results. The code used to perform the numerical experiments is openly
available.

Keywords: kinetic equation, variance prediction, statistical error, binomial
experiment, correlated samples

##  1 Introduction

The mesoscopic behaviour of large particle systems can be described by kinetic
equations that have a high-dimensional phase space. A common approach to deal
with the high-dimensionality of the problem is to employ Monte Carlo methods,
consisting of a simulation and estimation routine, for which the statistical
error converges with a rate that is independent of the dimensionality. For
linear kinetic equations with a linear collision operator, modelling for
instance collisions of the particles under consideration with a background
medium, particle tracing Monte Carlo methods [41, 30, 24] allow for unbiased
simulation and estimation of quantities of interest (QoI). Such particle
tracing Monte Carlo methods do, however, introduce a statistical error
(noise), which depends on the number of particles $N\in\mathbb{N}$ used in the
simulation. Because of the linearity of the kinetic equation, the $N$
particles can be simulated independently, leading to the well-known result
that the statistical error is inversely proportional to the square root of $N$
[9]. The precise value of the statistical error depends on the kinetic
equation under consideration as well as on the choice of simulation and
estimation method. For the illustration in this paper, we use the following
linear kinetic equation:

$\underbrace{\vphantom{\int}\partial_{t}f(x,v,t)}_{\begin{subarray}{c}\text{transient}\end{subarray}}+\underbrace{\vphantom{\int}v\cdot\partial_{x}f(x,v,t)}_{\begin{subarray}{c}\text{transport}\end{subarray}}=\underbrace{\vphantom{\int}S(x,v,t)}_{\begin{subarray}{c}\text{particle
source}\end{subarray}}-\underbrace{\vphantom{\int}R_{i}(x,t)f(x,v,t)}_{\begin{subarray}{c}\text{ionization
sink}\end{subarray}}+\underbrace{R_{cx}(x,t)\left(M(v\mid x,t)\int
f(x,v^{\prime},t)dv^{\prime}-f(x,v,t)\right)}_{\begin{subarray}{c}\text{charge-
exchange collision operator}\end{subarray}},$ (1)

which governs the neutral particle behaviour in the plasma edge of a fusion
device [35, 17, 30]. The particle velocity distribution $f(x,v,t)$ represents
the density of neutral particles at time $t\geq t_{0}=0$ with a given position
$x\in D\subset\mathbb{R}^{d}$ and velocity $v\in\mathbb{R}^{d}$; $S(x,v,t)$ is
the particle source describing how many neutrals are created at a position $x$
with velocity $v$ per unit of time; the ionization rate $R_{i}(x,t)$ and
charge-exchange rate $R_{cx}(x,t)$ are known positive functions that depend on
the background plasma state; the sum of the two collision rates is called the
total collision rate $R$:

$R(x,t)=R_{i}(x,t)+R_{cx}(x,t);$ (2)

$M(v\mid x,t)$ is the post-collisional velocity distribution describing the
probability of a neutral acquiring velocity $v$ after a charge-exchange
collision with a plasma particle. We here assume $M(v\mid x,t)$ to be a
Maxwellian:

$M(v\mid
x,t)=\frac{1}{(2\pi\sigma_{p}^{2}(x,t))^{d/2}}\exp\left(-\frac{1}{2}\frac{|v-u_{p}(x,t)|^{2}}{\sigma_{p}^{2}(x,t)}\right).$
(3)

The mean velocity $u_{p}(x,t)$ and variance $\sigma_{p}^{2}(x,t)$ (which we
assume to be isotropic, i.e., a scalar) of the post-collisional velocity
distribution are known from the background plasma state, rendering the kinetic
equation (1) linear. The QoI $Q$ that have to be estimated using Monte Carlo
estimation routines are typically moments of the particle velocity
distribution $f(x,v,t)$, i.e., integrals over velocity space of the following
form:

$Q(x,t)=\int q(v)f(x,v,t)dv,$ (4)

where $q(v)$ is a function of velocity. We focus on the neutral particle
density $\rho(x,t)$, momentum density $m(x,t)$, and energy density $E(x,t)$,
the three lowest order moments of the particle velocity distribution, which
are defined as follows:

$\rho(x,t)=\int f(x,v,t)dv,\quad m(x,t)=\int vf(x,v,t)dv,\quad
E(x,t)=\int\frac{|v|^{2}}{2}f(x,v,t)dv.$ (5)

Using particle tracing Monte Carlo methods, we can obtain unbiased estimates
$\hat{Q}$ that approximate the exact QoI $Q$. In this work, we search for
expressions that predict a priori the statistical error on these estimates
$\hat{Q}$ as a function of model parameters such as the collision rates
$R_{i}$, $R_{cx}$ and the post-collisional velocity distribution $M(v\mid
x,t)$. As error measure, we will consider the mean squared error, i.e., the
variance on the QoI estimates (this actually represents the statistical error
squared and therefore has an asymptotic convergence rate of $\frac{1}{N}$
[9]). Predicting the variance allows to determine the number of particles $N$
needed to achieve a desired accuracy a priori and gives a good indication of
the computational cost of a Monte Carlo simulation. Note that there might be
multiple particle sources (so-called strata) in a simulation that each require
a different number of particles [36]. Additionally, the variance predictors
can be used to select the best simulation and estimation combinations for a
given problem [29, 27, 18, 30] and to optimize particle tracing Monte Carlo
methods [37, 19] and hybrid simulation methods, e.g., domain decomposition
methods [7, 42] that encompass a Monte Carlo solver. Finally, there are some
other applications such as in multilevel Monte Carlo [14] where the variance
predictors can be used to determine the required number of particles per
level, and in Bayesian inversion settings [25] where the uncertainty of the
forward Monte Carlo simulations has to be taken into account [4].

In the literature, there is some work on analytical variance calculations,
based on invariant imbedding [40, 5, 30, 31, 32] and on the construction of
so-called history-score moment equations [2, 27, 37, 39]. These analytical
calculations build on derivations that only seem tractable in simple
(academic) test cases, e.g., one-group equations in infinite domains. In these
simple cases, closed form expressions for the variance can be obtained that
allow for selection and optimization of particle tracing Monte Carlo methods.
However, in more complex (realistic) cases such as equation (1), the
expressions are more involved such that they can only be solved numerically.
In those cases, and especially for higher dimensional problems, the numerical
variance calculations rapidly become impractical, which renders these methods
unattractive for a priori error assessment.

Most of the remaining literature focuses on on-the-fly variance estimation,
typically based on the assumption that the number of particles $N$ is high
enough such that asymptotic arguments hold [16] (this is not always easy to
verify). Then the following asymptotic expression for the variance on a QoI is
valid:

$\mathbb{V}[\hat{Q}]=\frac{C}{N},$ (6)

where estimation of the factor $C$ leads to a quantitative variance expression
[23, 13]. Empirical simulation driven variance estimation along these lines is
also performed for DSMC methods [10, 33]. Recently, error analysis has also
been performed for kernel density estimators for particle in cell methods
[12]. For Markov Chain Monte Carlo methods sampling an invariant distribution,
where consecutive samples are correlated, the correlation can be taken into
account in the on-the-fly variance estimation by also estimating the
integrated autocorrelation time [43, 15, 21] to obtain an effective number of
independent samples $N_{eff}$ that can be used in the asymptotic variance
expression (6). The discussed on-the-fly variance estimation techniques cannot
be used for a priori variance prediction and do not provide insight in the
dependence of the variance on model parameters, unless when repeating a lot of
simulations for different parameter settings and comparing the results, but
such an approach quickly becomes an expensive endeavour. Additionally, in
complex codes it might be unclear which error is dominant, e.g., when also a
finite sampling bias and discretization bias are present [13]. For such codes,
the conclusions of on-the-fly variance estimation might be unreliable.

In this paper, we show how the variance on the QoI estimates can be calculated
(expensive) and predicted (cheap) in the setting of an analog particle tracing
Monte Carlo method on a spatial grid, where each grid cell is used as a bin to
construct a histogram for the QoI. This setting is common, e.g., in plasma
edge modelling where equation (1) is solved to model the neutral particle
behaviour, and the grid (histogram) on which the QoI are estimated is also
used to simulate the background plasma using a finite volume method [35, 17,
30]. In Ref. [38], the idea is introduced that building a histogram can be
interpreted as performing a binomial experiment, where for each bin the
Bernoulli trial consists of a particle (or more generally speaking a sample)
being present in that bin (success) or not (failure). In that paper, only
independent samples are considered. However, for particle tracing Monte Carlo
methods a particle can contribute multiple times to the histogram, leading to
multiple correlated samples. Therefore, we extend the ideas introduced in Ref.
[38] to the setting of particle tracing Monte Carlo, where the construction of
a histogram can be interpreted as performing a correlated binomial experiment.

The literature on correlated binomial experiments (binary data) is focussed on
proposing a probability model (e.g., the Bahadur model) for the joint
probability of $N$ Bernoulli trials, and then fitting the parameters in that
model by regression using data (e.g., by using a maximum-likelihood estimator)
[11, 34, 26, 45]. This probability model for the joint probability
distribution then specifies the variance on the outcome of such a binomial
experiment. These probability models are, however, data-driven, which in our
simulation context means that it requires (expensive) simulations to generate
the required data. In this paper, we aim at providing an a priori variance
predictor that does not require (expensive) simulations. As we will see below,
several theory-based variance expressions exist for correlated binomial
experiments, some of which can be used as a cheap a priori variance predictor.
One variance expression (using hidden Markov model theory) allows to calculate
the actual variance, but at a high cost. Another variance expression only
accounts for the Markov dependence between consecutive samples and provides a
cheap approximate predictor for the actual variance. For each variance
expression, we verify the correctness and accuracy.

The rest of this paper is structured as follows. Section 2 introduces the
concepts underlying analog particle tracing Monte Carlo and QoI estimation. In
Section 3, different expressions to describe the variance on (correlated)
binomial experiments in an exact and approximate way are introduced. In
Section 4, we use the obtained expressions to calculate and predict the
variance on an analog particle tracing Monte Carlo method with point
estimators, in which each particle provides one independent contribution to
the QoI estimates. In Section 5, we use the obtained expressions to calculate
and predict the variance on an analog particle tracing Monte Carlo method with
analog estimators, for which each particle provides multiple correlated
contributions to the QoI estimates. Section 6 discusses some important remarks
with respect to the topic of statistical error prediction. Finally, the
conclusions and some pointers for future work are summarized in Section 7. The
code accompanying this paper is openly available in a public repository [1].

##  2 Particle dynamics and QoI estimators

In this section, we describe how kinetic equation (1) can be solved using an
analog particle tracing Monte Carlo scheme and introduce some statistical
assumptions on the particle dynamics that we will use to derive the a priori
variance predictors. Next, we show how QoI can be estimated on a grid and how
this estimation can be interpreted as the execution of a binomial experiment.
Finally, we illustrate how the variance on such a binomial experiment depends
on the model parameters of the kinetic equation.

###  2.1 Analog particle tracing Monte Carlo

Analog particle tracing Monte Carlo builds on the insight that the kinetic
equation (1) can be interpreted as the population dynamics of particles that
execute the following velocity-jump process [24]:

$\begin{split}&x_{n,k+1}=x_{n,k}+v_{n,k}\tau_{n,k},\quad x_{n,k}\in
D\subset\mathbb{R}^{d}\\\ &v_{n,k}\sim M(v_{n,k})=\mathcal{N}(v_{n,k}\mid
u_{p},\sigma_{p}^{2})=\frac{1}{(2\pi\sigma_{p}^{2})^{d/2}}\exp\left(-\frac{1}{2}\frac{|v_{n,k}-u_{p}|^{2}}{\sigma_{p}^{2}}\right),\quad
v_{n,k}\in\mathbb{R}^{d}\\\ &\tau_{n,k}\sim\mathcal{E}(\tau_{n,k}\mid
R)=R\exp(-R\tau_{n,k}),\quad\tau_{n,k}\in\mathbb{R}_{+}.\end{split}$ (7)

In this velocity-jump process, $x_{n,k}$ represents the position where
particle $1\leq n\leq N$ has its $k$-th collision ($x_{n,0}$ represents the
initial particle position at time $t_{0}=0$), $v_{n,k}$ is the velocity of the
particle after that collision, and $\tau_{n,k}$ is the flight time until the
next collision. At each collision point $x_{n,k}$, the trajectory is
terminated with probability

$\mathcal{P}_{t}=\frac{R_{i}}{R_{i}+R_{cx}}=\frac{R_{i}}{R}.$ (8)

The transition kernel corresponding to the velocity-jump process is given by

$\mathcal{K}(x_{n,k}\rightarrow
x_{n,k+1})=\int_{0}^{\infty}\frac{1}{(\tau_{n,k})^{d}}\cdot\mathcal{E}\left(\tau_{n,k}\mid
R\right)\mathcal{N}\left(\frac{x_{n,k+1}-x_{n,k}}{\tau_{n,k}}\mid
u_{p},\sigma_{p}^{2}\right)d\tau_{n,k}.$ (9)

The transition kernel (9) represents the probability density for having a
transition from $x_{n,k}$ to $x_{n,k+1}$ and intuitively follows from
integrating over all combinations of $\tau_{n,k}$ and $v_{n,k}$ that result in
this transition. Mathematically, the transition kernel is obtained by
differentiating the cumulative distribution function

$\begin{split}\mathcal{P}(v_{n,k}\tau_{n,k}\leq
x_{n,k+1}-x_{n,k})&=\int_{0}^{\infty}\mathcal{P}\left(v_{n,k}\leq\frac{x_{n,k+1}-x_{n,k}}{\tau_{n,k}}\mid\tau_{n,k}\right)\mathcal{E}(\tau_{n,k}\mid
R)d\tau_{n,k}\\\
&=\int_{0}^{\infty}\int_{-\infty}^{\frac{x_{n,k+1}-x_{n,k}}{\tau_{n,k}}}\mathcal{N}(v_{n,k}\mid
u_{p},\sigma_{p}^{2})dv_{n,k}\cdot\mathcal{E}(\tau_{n,k}\mid
R)d\tau_{n,k},\end{split}$ (10)

where we used the probability densities defined in (7).

When simulating kinetic equation (1) with velocity-jump process (7), the total
mass that has to be simulated,

$M=\int\int f(x,v,t=0)dxdv+\int\int\int S(x,v,t)dxdvdt,$ (11)

originates from the initial condition and/or the particle source and has to be
represented by $N$ particles. In an analog simulation, each particle gets the
same constant weight:

$w=\frac{M}{N}.$ (12)

There are two ways to treat the particle source $S(x,v,t)$. The first way
works for transient sources and deterministically discretizes the particle
source in time using $I$ time steps of size $\Delta t_{s}$. We then have
$N_{0}$ particles at the start of the simulation (at $t_{0}=0$) and generate
$N_{s,i}$ additional particles in each time step during the simulation. For
variance prediction, we can then simply treat each group (stratum) of
particles ($N_{0}$, $N_{s,i}$ $i=1..I$) separately.

The second way only works if the particle source is stationary
($S(x,v,t)=S(x,v)$) over a time interval $t\in[t_{1},t_{2}]$, meaning that the
times at which particles are generated are distributed uniformly in that time
interval. We can then split up the particles in a group (stratum) of $N_{0}$
particles that is present at the start of the simulation (at $t_{0}=0$) and a
group (stratum) of $N_{s}$ particles that is generated by the particle source
during the simulation. For each particle $n$ of those $N_{s}$ particles, we
then sample a generation time $t_{n}$ at which that particle enters the
simulation. The generation times are uniformly distributed in $[t_{1},t_{2}]$:

$\begin{split}t_{n}&\sim U[t_{1},t_{2}].\end{split}$ (13)

How to take such a stationary particle source into account for variance
prediction will be discussed below in sections 4 and 5.

###  2.2 Statistical assumptions

Note in (7) that a particle $n$ requires no information about other particles,
i.e., each particle is simulated independently from the others. This stems
from the fact that the kinetic equation (1) is linear. For the variance
prediction, we also assume that the simulation parameters $u_{p}$,
$\sigma_{p}^{2}$, $R_{cx}$, and $R_{i}$ are homogeneous (space independent)
and stationary (time independent), such that $\Delta t_{n,k}$, $v_{n,k}$, and
$x_{n,k}$ are independent random variables (probability distributions
$\mathcal{N}(v_{n,k}\mid u_{p},\sigma_{p}^{2})$ and
$\mathcal{E}(\tau_{n,k}\mid R)$ in (7) do not depend on $x$, $v$, $t$). From
these probability distributions, it also follows that both the velocities
$v_{n,k}$ and the flight times $\Delta t_{n,k}$ are independent and
identically distributed (iid). The positions $x_{n,k}$, however, are not iid,
because position $x_{n,k+1}$ depends on $x_{n,k}$. Note that the assumption of
constant simulation parameters is a typical finite volume assumption within
one grid cell of the simulation domain. In realistic simulation cases,
however, the simulation parameters do vary throughout the simulation domain
and in time. We neglect these space and time dependencies below in the
construction of the variance predictors.

###  2.3 QoI estimation on a grid

We estimate the QoI on a grid, which consists of $J$ cells $D_{j}\subset D$,
interpreting each grid cell as a bin of a histogram. We define an indicator
function $I_{j}(x)$ over each of the grid cells that equals 1 for $x\in D_{j}$
and 0 elsewhere. Note that a random evaluation of such an indicator function
is a Bernoulli random variable where 1 can be interpreted as a success and 0
as a failure. We can then define the QoI in cell $j$ as follows:

$Q_{j}(t)=\int_{D_{j}}Q(x,t)dx=\int_{D}Q(x,t)I_{j}(x)dx=\int_{D}\int
q(v)f(x,v,t)I_{j}(x)dvdx.$ (14)

A particle can only contribute to the estimate $\hat{Q}_{j}(t)$ of $Q_{j}(t)$
if its position is in cell $D_{j}$. We consider two QoI estimators: point
estimators and analog estimators. For a point estimator, we are interested in
the QoI at one specific point in time $t=T$. If we have $N$ particles, the
estimator is of the form

$Q_{j}(T)\approx\hat{Q}_{j}(T)=w\sum_{n=1}^{N}q(v_{n}(T))I_{j}(x_{n}(T)),$
(15)

where $x_{n}(T)$, $v_{n}(T)$ denote the position and velocity of particle $n$
at time $T$ without the constraint that it has to be a collision point. Note
that the estimate consists of $N$ independent contributions, as each particle
is simulated independently and only contributes once. For an analog estimator,
we are interested in the average QoI over a time interval $t\in[t_{1},t_{2}]$,
where we count contributions to the QoI at collision points $k\in[0,K(n)]$,
where $k=0$ represents the first collision in the time interval and $K(n)$
represents the last collision of particle $n$ in the time interval. If we have
$N$ particles, the estimator is of the form

$Q_{j}=\frac{1}{\Delta
t}\int_{t_{1}}^{t_{2}}Q_{j}(t)dt\approx\hat{Q}_{j}=\frac{w}{R\Delta
t}\sum_{n=1}^{N}\sum_{k=0}^{K(n)}q(v_{n,k})I_{j}(x_{n,k}),$ (16)

where we represent the duration of the time interval by $\Delta
t=t_{2}-t_{1}$. In this estimator, a particle $n$ delivers $K(n)$
contributions to the estimate, where the particle dynamics introduce
correlations between the positions $x_{n,k}$ in these contributions.

###  2.4 (Correlated) binomial experiments

In both estimators (15) and (16), we have a sum of indicator function
evaluations. We now write a generic sum of $L$ indicator function evaluations
as

$S_{L,j}=\sum_{l=1}^{L}I_{j}(x_{l}).$ (17)

This sum can be interpreted as a binomial experiment, in which each evaluation
of the indicator function $I_{j}(x_{l})$ represents a Bernoulli trial. We
denote the probability for success in Bernoulli trial $l$ by $p_{l,j}$. The
probability for success can be related to the expected value $\langle
I_{j}(x_{l})\rangle$ of the indicator function evaluation $I_{j}(x_{l})$ by
using the fact that an expected value of a discrete random variable can be
written as the sum of each possible outcome times the probability for that
outcome:

$p_{l,j}=\mathcal{P}(I_{j}(x_{l})=1)=1\cdot\mathcal{P}(I_{j}(x_{l})=1)+0\cdot\mathcal{P}(I_{j}(x_{l})=0)=\langle
I_{j}(x_{l})\rangle.$ (18)

The variance on the sum $S_{L,j}$ can be written as

$\begin{split}\mathbb{V}[S_{L,j}]=\mathbb{V}[\sum_{l=1}^{L}I_{j}(x_{l})]&=\langle(\sum_{l=1}^{L}I_{j}(x_{l})-\langle\sum_{l=1}^{L}I_{j}(x_{l})\rangle)^{2}\rangle\\\
&=\langle(\sum_{l=1}^{L}I_{j}(x_{l}))^{2}\rangle-\langle\sum_{l=1}^{L}I_{j}(x_{l})\rangle^{2}\\\
&=\langle\sum_{l=1}^{L}I_{j}(x_{l})^{2}+2\sum_{l=1}^{L}\sum_{m=1}^{l-1}I_{j}(x_{l})I_{j}(x_{m})\rangle-\langle\sum_{l=1}^{L}I_{j}(x_{l})\rangle^{2}\\\
&=\sum_{l=1}^{L}\langle
I_{j}(x_{l})^{2}\rangle+2\sum_{l=1}^{L}\sum_{m=1}^{l-1}\langle
I_{j}(x_{l})I_{j}(x_{m})\rangle-(\sum_{l=1}^{L}\langle
I_{j}(x_{l})\rangle)^{2}\\\ &=\sum_{l=1}^{L}\langle
I_{j}(x_{l})\rangle+2\sum_{l=1}^{L}\sum_{m=1}^{l-1}\langle
I_{j}(x_{l})I_{j}(x_{m})\rangle-(\sum_{l=1}^{L}\langle
I_{j}(x_{l})\rangle)^{2}\\\
&=\sum_{l=1}^{L}p_{l,j}+2\sum_{l=1}^{L}\sum_{m=1}^{l-1}\langle
I_{j}(x_{l})I_{j}(x_{m})\rangle-(\sum_{l=1}^{L}p_{l,j})^{2},\end{split}$ (19)

where in the penultimate equality, we used the following identity for
indicator functions:

$I_{j}(x_{l})^{2}=I_{j}(x_{l}).$ (20)

The goal of Section 3 is to find useful expressions for the variance
$\mathbb{V}[S_{L,j}]$, which we can then use to derive variance predictors for
the QoI estimators (15) and (16) in Section 4 and Section 5, respectively.

###  2.5 Variance as a function of model parameters

The variance on the outcome of a correlated binomial experiment representing
QoI estimation on a grid depends on model parameters such as $u_{p}$,
$\sigma_{p}^{2}$, $R_{i}$, and $R_{cx}$. To illustrate this, we compare three
realizations of two simulation set-ups with different simulation parameters in
Figure 1. The simulation parameters are chosen to be homogeneous and
stationary such that particles are distributed uniformly in space and time
throughout the simulation. In each realization, we build a histogram that
counts the number of collisions happening inside each bin for 1 simulated
particle, which resembles the result of an analog estimator (16) with $N=1$.
Each contribution to the histogram is a location where that particle has a
collision. The histograms have $J=10$ bins and the particle is allowed to
execute $L=100$ collisions. For both simulation set-ups, the expected number
of collisions in each bin is the same, namely $\frac{L}{J}=10$.

The simulation set-up in the left panel has a low value for the charge-
exchange rate $R_{cx}$, meaning that particles can travel a long distance in
between collisions. As a result, the histograms look as if the contributions
(the collision points) are independent of each other, leading to a small
variability between different realizations, which in turn leads to a low
variance on the realized histograms. The simulation set-up in the right panel
has a high value for the charge-exchange rate $R_{cx}$, meaning that particles
can only travel a short distance in between collisions. As a result, the
histograms are tightly centered around their initial collision point, because
consecutive collision points are strongly correlated. The histograms therefore
depend heavily on the initial position of the particle. The variability
between different realizations, with different initial particle positions, is
large, which in turn leads to a high variance on the realized histograms.

The experiment supports the following intuition about the variance of
histogram based QoI estimators in which each particle can deliver an equal
number of correlated contributions to an estimate: (1) model parameters that
result in particles only being able to travel short distances lead to a high
variance; (2) model parameters that result in particles being able to travel
long distances lead to a low variance. This intuition is confirmed
mathematically in the subsequent sections of this paper.

Figure 1: Particle paths (bottom) and corresponding histograms counting the
number of collisions happening inside each bin (top). Left: 3 realizations
with $u_{p}=0$, $\sigma_{p}^{2}=1$, $R_{i}=0$, $R_{cx}=1$. Right: 3
realizations with $u_{p}=0$, $\sigma_{p}^{2}=1$, $R_{i}=0$, $R_{cx}=25$. The
red dots indicate the initial particle positions. All simulations are
terminated after 100 collisions.

##  3 Variance expressions for (correlated) binomial experiments

In this section, we derive four expressions for the variance
$\mathbb{V}[S_{L,j}]$ (19) on a (correlated) binomial experiment consisting of
$L$ Bernoulli trials. First, we derive an upper bound on $\mathbb{V}[S_{L,j}]$
by searching for the worst-case relation between consecutive Bernoulli trials
$x_{l}$; this upper bound is independent of any model parameters of the
kinetic equation (1). Second, we derive an expression for a binomial
experiment in which consecutive Bernoulli trials are assumed to be
independent; this expression is also independent of any model parameters.
Third, we assume that consecutive Bernoulli trials have a Markov dependence,
i.e., we assume that there is a 2-state Markov process that determines if a
particle is inside bin $D_{j}$ or not. The probabilities involved in the
assumed Markov process can be related to the bin size and model parameters of
the kinetic equation (1). Finally, we assume that the binomial experiment is
the observed process of a hidden Markov model, where the underlying (hidden)
process is the particle dynamics (7). The involved probabilities can again be
related to the bin size and model parameters of the kinetic equation.

###  3.1 Upper bound

Let us start by providing a simple upper bound for the variance on $S_{L,j}$.
We look for stationary worst-case values $p$, $z$ for the expected values in
the variance expression (19):

$\begin{split}&p=\langle I_{j}(x_{l})\rangle=p_{l,j},\\\ &z=\langle
I_{j}(x_{l})I_{j}(x_{m})\rangle.\end{split}$ (21)

Because the indicator functions can only take values 0 or 1, we have that
$0\leq p,z\leq 1$. We insert $p$, $z$ in (19) and solve the following
optimization problem to find an upper bound for the variance on $S_{L,j}$:

$\max_{p,z}Lp+L(L-1)z-L^{2}p^{2}\quad\text{s.t.}\quad 0\leq p,z\leq 1,$ (22)

where $L\geq 1$ is the number of Bernoulli trials in the binomial experiment.
It can be verified that the solution for this maximization problem is
$p=\frac{1}{2L}$, $z=1$, resulting in the following upper bound for the
variance:

$\mathbb{V}[S_{L,j}]=\mathbb{V}[\sum_{l=1}^{L}I_{j}(x_{l})]\leq
L(L-1)+0.25\sim\mathcal{O}(L^{2}),$ (23)

which is quadratic in the number of Bernoulli trials $L$. As we will see in
the numerical experiments below, the upper bound (23) typically provides a
huge overestimation of the actual variance of a binomial experiment.

#### Remark 1

A quadratic variance $\sim L^{2}$ on a binomial experiment means that for an
empirical mean of the form $\frac{1}{L}\sum_{l=1}^{L}I_{j}(x_{l})$, the
variance is independent of $L$, i.e., the variance does not decrease for
$L\rightarrow\infty$. This corresponds to taking $L$ times the same sample in
each realization of the sum (perfect correlation, i.e., $z=1$).

###  3.2 Binomial experiment with independent Bernoulli trials

Let us assume that the Bernoulli trials are independent and that each trial
occurs with the same (stationary) success probability, meaning that the
samples of the random variable $I_{j}(x_{l})$ are iid:

$\begin{split}&\langle I_{j}(x_{l})\rangle=p_{l,j}=p_{j},\\\ &\langle
I_{j}(x_{l})I_{j}(x_{m})\rangle=\langle I_{j}(x_{l})\rangle\langle
I_{j}(x_{m})\rangle=p_{l,j}p_{m,j}=p_{j}^{2}.\end{split}$ (24)

Variance expression (19) then simplifies to the form

$\mathbb{V}[S_{L,j}]=\mathbb{V}[\sum_{l=1}^{L}I_{j}(x_{l})]=Lp_{j}+L(L-1)p_{j}^{2}-L^{2}p_{j}^{2}=Lp_{j}(1-p_{j}),$
(25)

which is the well known expression for the variance of a binomial experiment
with independent trials. Note that the variance expression for independent
trials (25) is linear in $L$, while the upper bound (23) derived above is
quadratic. Information on how to obtain a value for $p_{j}$ can be found in
Section 6.1.

###  3.3 Binomial experiment driven by a Markov process

In this section, we assume that the binomial experiment for each bin $D_{j}$
is driven by a 2-state Markov process in which the current Bernoulli trial
depends on the outcome of the previous trial (but not on outcomes of older
trials). Binomial experiments with a Markov dependence have been investigated
in a stationary setting in Ref. [22] and in a transient setting in Ref. [44].
The Markov dependence of the current trial on the previous trial is encoded by
conditional probabilities. In the transient setting, the (conditional)
probabilities are allowed to have a different value for each Bernoulli trial,
while in the stationary setting, the same probabilities are assumed for each
trial. The (trial dependent) probability for success is defined above as
$p_{l,j}$ in (18). The probability for failure follows logically as
$\mathcal{P}(I_{j}(x_{l})=0)=1-\mathcal{P}(I_{j}(x_{l})=1)=1-p_{l,j}$. To
fully characterize the Markov dependence, we define the following conditional
probabilities encoding the dependence on the outcome of the previous trial:

$\begin{split}&\mathcal{P}(I_{j}(x_{l})=1\mid
I_{j}(x_{l-1})=1)=\lambda_{l-1,j},\\\ &\mathcal{P}(I_{j}(x_{l})=0\mid
I_{j}(x_{l-1})=1)=1-\mathcal{P}(I_{j}(x_{l})=1\mid
I_{j}(x_{l-1})=1)=1-\lambda_{l-1,j},\\\ &\mathcal{P}(I_{j}(x_{l})=0\mid
I_{j}(x_{l-1})=0)=\nu_{l-1,j},\\\ &\mathcal{P}(I_{j}(x_{l})=1\mid
I_{j}(x_{l-1})=0)=1-\mathcal{P}(I_{j}(x_{l})=0\mid
I_{j}(x_{l-1})=0)=1-\nu_{l-1,j}.\end{split}$ (26)

The conditional probabilities serve as transition probabilities in the Markov
process that drives the binomial experiment. We can construct the transition
matrix $T_{j}(P_{l-1}\rightarrow P_{l})$ as follows:

$T_{j}(P_{l-1}\rightarrow
P_{l})=\begin{bmatrix}\nu_{l-1,j}&1-\lambda_{l-1,j}\\\
1-\nu_{l-1,j}&\lambda_{l-1,j}\end{bmatrix},$ (27)

defining the vector $P_{l}=\begin{bmatrix}1-p_{l,j}\\\ p_{l,j}\end{bmatrix}$
containing the failure and success probability for trial $l$. We then recover
the following Markov process:

$P_{l}=T_{j}(P_{l-1}\rightarrow P_{l})P_{l-1}.$ (28)

Having defined the Markov process and the probabilities involved, we can
elaborate an expression for the variance on $S_{L,j}$ (19). Following Ref.
[44], we observe that in the second term in the variance expression (19), we
can rewrite $\langle I_{j}(x_{l})I_{j}(x_{m})\rangle$ with $m<l$ as follows:

$\begin{split}\langle I_{j}(x_{l})I_{j}(x_{m})\rangle&=\langle
I_{j}(x_{l})I_{j}(x_{m})\mid
I_{j}(x_{l})=0,I_{j}(x_{m})=0\rangle\cdot\mathcal{P}(I_{j}(x_{l})=0,I_{j}(x_{m})=0)\\\
&+\langle I_{j}(x_{l})I_{j}(x_{m})\mid
I_{j}(x_{l})=1,I_{j}(x_{m})=0\rangle\cdot\mathcal{P}(I_{j}(x_{l})=1,I_{j}(x_{m})=0)\\\
&+\langle I_{j}(x_{l})I_{j}(x_{m})\mid
I_{j}(x_{l})=0,I_{j}(x_{m})=1\rangle\cdot\mathcal{P}(I_{j}(x_{l})=0,I_{j}(x_{m})=1)\\\
&+\langle I_{j}(x_{l})I_{j}(x_{m})\mid
I_{j}(x_{l})=1,I_{j}(x_{m})=1\rangle\cdot\mathcal{P}(I_{j}(x_{l})=1,I_{j}(x_{m})=1)\\\
&=0\times\mathcal{P}(I_{j}(x_{l})=0,I_{j}(x_{m})=0)+0\cdot\mathcal{P}(I_{j}(x_{l})=1,I_{j}(x_{m})=0)\\\
&+0\cdot\mathcal{P}(I_{j}(x_{l})=0,I_{j}(x_{m})=1)+1\cdot\mathcal{P}(I_{j}(x_{l})=1,I_{j}(x_{m})=1)\\\
&=\mathcal{P}(I_{j}(x_{l})=1,I_{j}(x_{m})=1)\\\
&=\begin{bmatrix}0&1\end{bmatrix}\left(\prod_{i=m+1}^{l}T_{j}(P_{i-1}\rightarrow
P_{i})\right)\begin{bmatrix}0\\\ p_{m,j}\end{bmatrix}.\end{split}$ (29)

This result can be interpreted as follows (by reading the final expression
from the right to the left). We take the fraction $p_{m,j}$ of the Markov
chains that have success as outcome in trial $m$. Then, we look at how this
fraction redistributes its probability in each of the next $l-m$ trials, to
find which fraction of those chains ends up with a success as outcome in trial
$l$ and which fraction has a failure as outcome in trial $l$. Finally, we are
only interested in the fraction of Markov chains that have a success in trial
$l$. Inserting (29) into the variance expression (19) results in

$\mathbb{V}[S_{L,j}]=\mathbb{V}[\sum_{l=1}^{L}I_{j}(x_{l})]=\sum_{l=1}^{L}p_{l,j}+2\sum_{l=1}^{L}\sum_{m=1}^{l-1}\begin{bmatrix}0&1\end{bmatrix}\left(\prod_{i=m+1}^{l}T_{j}(P_{i-1}\rightarrow
P_{i})\right)\begin{bmatrix}0\\\
p_{m,j}\end{bmatrix}-\left(\sum_{l=1}^{L}p_{l,j}\right)^{2},$ (30)

which can also be written in the following form [44]:

$\begin{split}\mathbb{V}[S_{L,j}]=\mathbb{V}[\sum_{l=1}^{L}I_{j}(x_{l})]&=\sum_{l=1}^{L}p_{l,j}(1-p_{l,j})\\\
&+2\sum_{l=1}^{L}\sum_{m=1}^{l-1}\left(\begin{bmatrix}0&1\end{bmatrix}\left(\prod_{i=m+1}^{l}T_{j}(P_{i-1}\rightarrow
P_{i})\right)\begin{bmatrix}0\\\
p_{m,j}\end{bmatrix}-p_{l,j}p_{m,j}\right).\end{split}$ (31)

In this last expression, we recognize the variance for a transient independent
binomial experiment (compare the first term to the expression for the variance
of an independent stationary binomial experiment (25)). The second term can
therefore be considered as a correction term that corrects for the presence of
a Markov dependence between consecutive trials.

For stationary processes, we can drop the subscript $l$ on the probabilities
$p_{j}$, $\lambda_{j}$ and $\nu_{j}$. As the probabilities should remain the
same in each trial, we obtain the following detailed balance condition:

$\mathcal{P}(I_{j}(x_{l})=1,I_{j}(x_{l-1})=0)=\mathcal{P}(I_{j}(x_{l})=0,I_{j}(x_{l-1})=1),$
(32)

from which it follows that

$1-\nu_{j}=\mathcal{P}(I_{j}(x_{l})=1\mid
I_{j}(x_{l-1})=0)=\frac{\mathcal{P}(I_{j}(x_{l})=0\mid
I_{j}(x_{l-1})=1)\cdot\mathcal{P}(I_{j}(x_{l-1})=1)}{\mathcal{P}(I_{j}(x_{l-1})=0)}=\frac{(1-\lambda_{j})p_{j}}{1-p_{j}},$
(33)

i.e., $\nu_{j}$ becomes a function of $p_{j}$ and $\lambda_{j}$, meaning that
the last two now fully characterize the 2-state Markov process [22].
Information on how to obtain values for $p_{j}$ and $\lambda_{j}$ for a given
problem can be found in Ref. [20] and in Section 6.1. We now obtain the
stationary Markov process $P_{l}=T_{j}P_{l-1}$ with stationary transition
matrix $T_{j}$:

$T_{j}=\begin{bmatrix}1-\frac{(1-\lambda_{j})p_{j}}{1-p_{j}}&1-\lambda_{j}\\\
\frac{(1-\lambda_{j})p_{j}}{1-p_{j}}&\lambda_{j}\end{bmatrix}.$ (34)

The stationary equivalent of (29) is as follows:

$\langle
I_{j}(x_{l})I_{j}(x_{m})\rangle=\begin{bmatrix}0&1\end{bmatrix}T_{j}^{l-m}\begin{bmatrix}0\\\
p_{m,j}\end{bmatrix}=p_{j}\left(p_{j}+(1-p_{j})\left(\frac{\lambda_{j}-p_{j}}{1-p_{j}}\right)^{l-m}\right).$
(35)

The expression for the variance on $S_{L,j}$ then simplifies from (31) to

$\begin{split}\mathbb{V}[S_{L,j}]=\mathbb{V}[\sum_{l=1}^{L}I_{j}(x_{l})]&=Lp_{j}(1-p_{j})\\\
&+\frac{2p_{j}(1-p_{j})(\lambda_{j}-p_{j})}{1-\lambda_{j}}\left(L-\frac{1-p_{j}}{1-\lambda_{j}}\left(1-\left(\frac{\lambda_{j}-p_{j}}{1-p_{j}}\right)^{L}\right)\right).\end{split}$
(36)

Note that this is a slightly different result than expression (5.5) in Ref.
[22], which contained an error. The variance expression once again splits up
into a term representing the variance for an uncorrelated binomial experiment
with independent trials and a correction term that takes the Markov dependence
into account. Expression (36) can be used as a cheap to compute a priori
predictor for the variance on QoI estimates. In the remainder, we refer to
(36) as the Markov process based variance predictor.

We draw attention to three important points in variance expression (36).
First, all the conditional probabilities in the transition matrix (34) should
be values between 0 and 1. This leads to the following constraint:

$\max\left(\frac{2p_{j}-1}{p_{j}},0\right)\leq\lambda_{j}\leq 1,$ (37)

which naturally constitutes a feasibility border that separates feasible and
infeasible combinations of $p_{j}$ and $\lambda_{j}$. Under this constraint it
follows that

$\left|\frac{\lambda_{j}-p_{j}}{1-p_{j}}\right|\leq 1,$ (38)

meaning that the exponential part
$\left(\frac{\lambda_{j}-p_{j}}{1-p_{j}}\right)^{L}$ in expression (36) cannot
blow up. Second, the correction term scales with $\frac{1}{1-\lambda_{j}}$ and
for $\lambda_{j}\rightarrow 1$, with $0<p_{j}<1$, the variance seems to blow
up to infinity, which is not realistic for the variance of a finite sum.
However, working out the limit $\lambda_{j}\rightarrow 1$ using l’Hôpital’s
rule shows that the actual limit of the variance is

$\lim_{\lambda_{j}\rightarrow
1}\mathbb{V}[\sum_{l=1}^{L}I_{j}(x_{l})]=L^{2}p_{j}(1-p_{j}),$ (39)

which is a quadratic expression in $L$, just like the upper bound (23). It can
be verified that expression (39) respects the upper bound (23) for all $L\geq
1$. The limit $\lambda_{j}\rightarrow 1$ corresponds to transitions between
success and failure becoming rare events, leading to a sharp increase in the
variance. Third, for $\lambda_{j}<1$, we observe that the expression for the
variance grows linearly with $L$, just as the variance of a binomial
experiment with independent trials (25) (the quadratic behaviour (39) is only
present for $\lambda_{j}\rightarrow 1$). For large values of $L$ and
$\lambda_{j}<1$, the variance expression (36) simplifies approximately to

$\mathbb{V}[\sum_{l=1}^{L}I_{j}(x_{l})]\approx
L\cdot\left(p_{j}(1-p_{j})+\frac{2p_{j}(1-p_{j})(\lambda_{j}-p_{j})}{1-\lambda_{j}}\right).$
(40)

In Figure 2, we plot the variance expression (40) divided by $L$ for
$\lambda_{j}\leq 0.9$. If $\lambda_{j}=p_{j}$, then the variance expression
reduces to the expression for independent trials. This is logical, because it
means that knowledge of the previous outcome does not affect the probability
for success. If $\lambda_{j}<p_{j}$, which corresponds to more mixing in the
Markov chain compared to when doing independent trials, then the correction
term is negative, meaning that the Markov dependence lowers the variance on
the binomial experiment. If $\lambda_{j}>p_{j}$, which corresponds to more
clustering in the Markov chain compared to when doing independent trials, then
the correction term is positive, meaning that the Markov dependence increases
the variance on the binomial experiment. As discussed above, the variance
increases sharply for $\lambda_{j}\rightarrow 1$ and $0<p_{j}<1$. E.g., for
$\lambda_{j}=0.9$, $p_{j}=0.5$, we see that the variance per trial is inflated
from $0.25$ for the independent case to about $2.5$ for the Markov dependence
case, showing a variance inflation of about a factor of 10.

#### Remark 2

The conditional probability $\lambda_{j}$ of staying within bin $j$ increases
if particles move over shorter distances between collisions [20]. The observed
variance increase for increasing $\lambda_{j}$ therefore corresponds to the
experimental observations from Section 2.5.

Figure 2: Markov dependence based variance expression for large $L$ (40)
divided by $L$ for $0\leq p_{j}\leq 1$ and $0\leq\lambda_{j}\leq 0.9$. For
large values of $\lambda_{j}$, transitions between success and failure become
rare events, inflating the variance. The blue line represents the variance of
a binomial experiment with $\lambda_{j}=p_{j}$, i.e., independent Bernoulli
trials (25) divided by $L$. The red line represents the feasibility border
(37). Left: top view. Right: 3D view.

###  3.4 Binomial experiment driven by a hidden Markov process

In this section, we assume that the binomial experiment is driven by a hidden
Markov process, which corresponds exactly to the setting of estimating QoI on
a histogram. When we build estimators on a histogram, we are actually coarse-
graining a (hidden) continuous Markov process for the position of a particle
$x_{l}$ to a binomial experiment monitoring the presence of that particle in a
bin $j$, represented by random variable $I_{j}(x_{l})$ (an evaluation of an
indicator function). This corresponds to the description of a hidden Markov
model, where the presence of particles in a given bin, i.e., the binomial
experiment, can be interpreted as the observed process, and the stochastic
particle dynamics can be interpreted as the underlying hidden Markov process.
In this case, however, the hidden Markov process is not really ‘hidden’ as we
know the particle dynamics. If the particle dynamics are governed by the
transition kernel $\mathcal{K}(x_{l-1}\rightarrow x_{l})$ (e.g., (9) for
particle dynamics (7)), then this transition kernel can be identified to be
the conditional probability density $\mathcal{P}(x_{l}\mid x_{l-1})$. If we
then have the probability density $P_{l}(x_{l})=\mathcal{P}(x_{l}=x)$, we can
write the following continuous Markov process to go from a probability density
on state (particle position) $x_{l-1}$ to a probability density on state
(particle position) $x_{l}$ [41]:

$P_{l}(x_{l})=\int\mathcal{K}(x_{l-1}\rightarrow x_{l})P_{l-1}(x_{l-1})dx.$
(41)

This Markov process (a recursion relation that only depends on state
$x_{l-1}$, not on older states) is a continuous equivalent of the discrete
Markov process (28). We can use expression (41) to calculate the different
terms in the expression for the variance of $S_{L,j}$ (19). The first term in
that expression becomes

$\begin{split}\sum_{l=1}^{L}\langle
I_{j}(x_{l})\rangle&=\sum_{l=1}^{L}\mathcal{P}(I_{j}(x_{l})=1)\\\
&=\sum_{l=1}^{L}\int\ldots\int I_{j}(x_{l})\mathcal{K}(x_{l-1}\rightarrow
x_{l})\ldots\mathcal{K}(x_{1}\rightarrow x_{2})P_{1}(x_{1})dx_{1}\ldots
dx_{l},\end{split}$ (42)

i.e., for each value of $l$, we take the probability density $P_{1}(x_{1})$ on
the initial state, propagate it using the transition kernels to the
probability density on the $l$-th state, and then integrate that probability
density over bin $j$, indicated by indicator function $I_{j}(x_{l})$, to count
the expected fraction of Markov chains that have state $x_{l}$ in bin $j$.
Following a similar reasoning, we find the analog for (29) to be

$\begin{split}\langle
I_{j}(x_{l})I_{j}(x_{m})\rangle&=\mathcal{P}(I_{j}(x_{l})=1,I_{j}(x_{m})=1)\\\
&=\int\ldots\int I_{j}(x_{l})I_{j}(x_{m})\mathcal{K}(x_{l-1}\rightarrow
x_{l})\ldots\mathcal{K}(x_{1}\rightarrow x_{2})P_{1}(x_{1})dx_{1}\ldots
dx_{l},\end{split}$ (43)

with $m<l$. The main difference between (29) and (43) is that in the former,
because of the Markov dependence, the probabilities for the current state only
depend on the previous state, while in the latter, the coarse-grained (hidden)
Markov process, the random variables $I_{j}(x_{l})$ do not necessarily have
the Markov property anymore, meaning that we need to take all the transitions
of the (hidden) continuous Markov process into account starting from the
initial state [28].

The high dimensional integrals (42), (43) are typically not computable
analytically. Instead, we discretize the (hidden) continuous Markov process by
subdividing the continuous phase space into small cells and using these cells
as discrete states. We assume that particles are distributed uniformly in each
of the cells, which enables the computation of transition probabilities
between the different cells, as explained in Ref. [20]. The indicator function
$I_{j}(x_{l})$ then becomes a diagonal matrix $\mathbb{I}_{l,j}$ with the
value 1 on the diagonal elements corresponding to the cells that belong to bin
$j$, and the value 0 on the diagonal elements corresponding to the cells that
are outside bin $j$. The probability density $P_{1}(x)$ becomes a vector
$P_{1}$ with at each entry the probability of the initial state being in the
cell corresponding to that element. We obtain a discrete Markov process of the
form (28). The discretized versions of expressions (42), (43) are as follows
[28]:

$\begin{split}\sum_{l=1}^{L}\langle
I_{j}(x_{l})\rangle&=\sum_{l=1}^{L}\mathcal{P}(I_{j}(x_{l})=1)\\\
&=\sum_{l=1}^{L}\left(\boldsymbol{1}^{\prime}\mathbb{I}_{l,j}\left(\prod_{i=2}^{l}T_{j}(P_{i-1}\rightarrow
P_{i})\right)P_{1}\right),\end{split}$ (44)

$\begin{split}\langle
I_{j}(x_{l})I_{j}(x_{m})\rangle&=\mathcal{P}(I_{j}(x_{l})=1,I_{j}(x_{m})=1)\\\
&=\boldsymbol{1}^{\prime}\mathbb{I}_{l,j}\left(\prod_{i=m+1}^{l}T_{j}(P_{i-1}\rightarrow
P_{i})\right)\mathbb{I}_{l,j}\left(\prod_{i=2}^{m}T_{j}(P_{i-1}\rightarrow
P_{i})\right)P_{1},\end{split}$ (45)

where $\boldsymbol{1}$ is a vector of ones. To fully appreciate the difference
with the Markov dependence case from Section 3.3, equation (45) should be
compared to (29).

The high dimensional integrals (42), (43) and their discretizations (44), (45)
are expensive to compute and therefore unfit for a priori variance prediction.
There is, however, a way to obtain a cheap approximate expression for the
variance on $S_{L,j}$. A typical simplification for a coarse-grained (hidden)
Markov process, is to only consider the Markovian component of the observed
dynamics [3]. We then retrieve the case of a binomial experiment driven by a
Markov process discussed in Section 3.3, which means that we can use the
expression (36) as a cheap predictor of the variance on $S_{L,j}$, provided
that we have expressions for $p_{j}$ and $\lambda_{j}$.

###  3.5 Numerical experiments

We test the derived expressions for a range of values for $R$ and
$\sigma_{p}^{2}$, while setting $u_{p},R_{i}=0$. For each set of parameters,
we simulate a single particle for $L=10^{3}$ collisions and count how many
collision points (Bernoulli trials) are in a bin $j$, which corresponds to the
setting of an analog estimator with $N=1$, see (16). The test cases are all
executed on a 1D ($d=1$) domain $D=[0,1]$ with length 1 and periodic boundary
conditions. All particles are launched at time $t=0$, with an initial position
distributed uniformly over the domain and a velocity drawn from the post-
collisional velocity distribution given in (7). The space is subdivided in
$J=10$ bins, such that the probability for success equals
$p_{j}=\frac{1}{J}=0.1$. The subdomain $D_{j}=[0,0.1]$ is considered as the
bin of interest. We compare the upper bound (23) for the variance, the
variance expression for independent trials (25), the Markov process based
variance expression (36), and the discretized hidden Markov process based
variance expression. The hidden Markov process is discretized on an
equidistant grid consisting of 100 cells. The conditional probability
$\lambda_{j}$ for the Markov process and transition probabilities for the
discretized hidden Markov process can be calculated as a function of the model
parameters as described in Ref. [20]. The code is openly available in a public
repository [1].

We consider three experiments: the first experiment is in the so-called
hydrodynamic scaling (vary $R$ over orders of magnitude while keeping
$\sigma_{p}^{2}\sim\mathcal{O}(1)$); the second experiment is in the so-called
temperature scaling (vary $\sigma_{p}^{2}$ over orders of magnitude while
keeping $R\sim\mathcal{O}(1)$); and the third experiment is in the so-called
diffusive scaling (choose $\sigma_{p}^{2}=R$ and vary over orders of
magnitude). As a reference solution, we determine for each set of parameters
the empirical variance based on $10^{4}$ realizations of the experiment. We
plot the different variance predictors divided by $L$ for these three scalings
in Figure 3, where we each time use $\sigma_{p}^{2}=1$ in the hydrodynamic
scaling and $R=1$ in the temperature scaling. The upper bound for the variance
is omitted in the temperature scaling plot, to zoom in on the other variance
expressions.

As expected, parameter settings that result in particles being able to travel
longer distances (low $R$, high $\sigma_{p}^{2}$) correspond to lower
variances, while parameter settings that result in particles being able to
travel shorter distances (high $R$, low $\sigma_{p}^{2}$) correspond to higher
variances. The discretized hidden Markov process based expression correctly
predicts the empirically variance, except for a small discretization error. In
each plot, we can see that using the Markov process based expression for the
variance (36) is an enormous improvement over the assumption of independent
trials, except in the temperature scaling. The peculiar behaviour in the
temperature scaling can be explained by a boundary condition effect as
follows. Due to increasing $\sigma_{p}^{2}$ for constant $R$, particles can
travel further and further in one step. Because of the periodic boundary
conditions, this means that at some point particles can fly around the whole
domain returning to the bin where they left. In the limit for very large step
sizes, the periodic boundary conditions make that the variance on the binomial
experiment evolves to the case with independent samples, i.e., to
$\lambda_{j}=p_{j}$. The conditional probability $\lambda_{j}$, however, is
calculated as described in Ref. [20] and does not take the periodic boundary
conditions into account. The calculated $\lambda_{j}$ therefore corresponds to
the case in which a particle that escapes from a cell cannot end up in that
cell again, leading to $\lambda_{j}<p_{j}$ and an underestimation of the
variance. Taking boundary conditions into account in the variance predictors
is left for future work.

Figure 3: Variance divided by $L$ for $L=10^{3}$ Bernoulli trials. Left:
hydrodynamic scaling. Middle: temperature scaling. Right: diffusive scaling.
The empirical variance (blue) is calculated based on $10^{4}$ realizations.
The upper bound (red) is given by (23). The expression for independent
Bernoulli trials (yellow) is given by (25). The used Markov dependence
approximation (green) is given by (36). The variance expression based on a
discretized hidden Markov process (purple) is explained in Section 3.4.

###  3.6 Towards variance expressions for QoI estimators

The expressions obtained above deal with the variance on sums $S_{L,j}$ of the
form (17). The QoI estimators (15) and (16), however, contain more general
sums of the form

$S_{L,x,j}=w\sum_{l=1}^{L}q_{x}(v_{l})I_{j}(x_{l}),$ (46)

where for the lowest order moments (5), we have $q_{\rho}(v_{l})=1$,
$q_{m}(v_{l})=v_{l}$, and $q_{E}(v_{l})=|v_{l}|^{2}/2$. We compute the
variance on such a sum for a fixed value of $L$ by performing a similar
derivation as in (19). Using the fact that the velocities $v_{l}$ are iid, and
that $v_{l}$ and $x_{l}$ are independent and therefore functions of $x_{l}$
and functions of $v_{l}$ are also independent, we obtain

$\mathbb{V}[S_{L,x,j}]=w^{2}\mathbb{V}[q_{x}(v_{l})]\langle
S_{L,j}\rangle+w^{2}\langle q_{x}(v_{l})\rangle^{2}\mathbb{V}[S_{L,j}].$ (47)

The statistics on $q_{x}(v_{l})$ can be calculated given the probability
density $M(v_{l})$ in (7). For the three lowest order moments (5), we obtain

$\begin{split}&\langle q_{\rho}(v_{l})\rangle=1,\quad\hskip 5.69054pt\langle
q_{m}(v_{l})\rangle=|u_{p}|,\quad\hskip 5.69054pt\langle
q_{E}(v_{l})\rangle=\frac{1}{2}(|u_{p}|^{2}+\sigma_{p}^{2}),\\\
&\mathbb{V}[q_{\rho}(v_{l})]=0,\quad\mathbb{V}[q_{m}(v_{l})]=\sigma_{p}^{2},\quad\mathbb{V}[q_{E}(v_{l})]=|u_{p}|^{2}\sigma_{p}^{2}+\frac{1}{2}\sigma_{p}^{4}.\end{split}$
(48)

From above, we have an upper bound (23), an expression for independent trials
(25), Markov process based expressions, and hidden Markov process based
expressions for the variance $\mathbb{V}[S_{L,j}]$. We therefore have all the
required ingredients to construct variance predictors for analog particle
tracing Monte Carlo methods with point estimators and analog estimators.

##  4 Predicting the variance of analog particle tracing Monte Carlo with
point estimators

For point estimators (15), the contributions (Bernoulli trials) are
independent. The variance on the point estimator $\hat{Q}_{j}(T)$ can be
written in the form (47) as follows:

$\begin{split}\mathbb{V}[\hat{Q}_{j}(T)]&=w^{2}\mathbb{V}[q_{x}(v_{n}(T))]\langle\sum_{n=1}^{N}I_{j}(x_{n}(T))\rangle+w^{2}\langle
q_{x}(v_{n}(T))\rangle^{2}\mathbb{V}[\sum_{n=1}^{N}I_{j}(x_{n}(T))]\\\
&=w^{2}Np_{j}\mathbb{V}[q_{x}(v_{n}(T))]+w^{2}\langle
q_{x}(v_{n}(T))\rangle^{2}\mathbb{V}[\sum_{n=1}^{N}I_{j}(x_{n}(T))],\end{split}$
(49)

where the last equality only holds for a stationary probability for success
$p_{j}$. Because of the independent contributions, we can use expression (25)
for the variance on the binomial experiment with independent trials:

$\mathbb{V}[\hat{Q}_{j}(T)]=w^{2}Np_{j}\mathbb{V}[q_{x}(v_{n}(T))]+w^{2}Np_{j}(1-p_{j})\langle
q_{x}(v_{n}(T))\rangle^{2}.$ (50)

Working out this expression explicitly for the three lowest order moments (5)
by inserting (48) results in

$\begin{split}\mathbb{V}[\hat{\rho}_{j}(T)]&=w^{2}Np_{j}(1-p_{j}),\\\
\mathbb{V}[\hat{m}_{j}(T)]&=w^{2}Np_{j}\sigma_{p}^{2}+w^{2}Np_{j}(1-p_{j})|u_{p}|^{2},\\\
\mathbb{V}[\hat{E}_{j}(T)]&=w^{2}Np_{j}(|u_{p}|^{2}\sigma_{p}^{2}+\frac{1}{2}\sigma_{p}^{4})+\frac{1}{4}w^{2}Np_{j}(1-p_{j})(|u_{p}|^{2}+\sigma_{p}^{2})^{2}.\end{split}$
(51)

Inserting (12) for the particle weight $w$ shows that the variances scale with
$\frac{1}{N}$, as they should.

###  4.1 Influence of sources and sinks on probability for success $p_{j}$

If there are no sinks ($R_{i}=0$), then we can select a probability for
success $p_{j}=\bar{p}_{j}$, e.g., by assuming that a success in each bin $j$
occurs with the same probability: $\bar{p}_{j}=\frac{1}{J}$, such that
$\sum_{j=1}^{J}p_{j}=\sum_{j=1}^{J}\bar{p}_{j}=1$. If there is an ionization
sink, then we have that $\sum_{j=1}^{J}p_{j}<1$, because a particle can
disappear in the sink before reaching time $T$. We can incorporate the sink in
the probability for success as follows:

$\begin{split}p_{j}&=\mathcal{P}(\text{particle ends up in cell $j$ and is not
ionized})\\\ &=\mathcal{P}(\text{particle ends up in cell
$j$})\cdot\mathcal{P}(\text{particle is not ionized})\\\
&=\bar{p}_{j}\cdot\exp(-R_{i}T).\end{split}$ (52)

Note that each particle still has the same probability for success $p_{j}$.
When the particles originate from a stationary particle source, as described
in Section 2.1, then a particle $n$ that survives (does not get ionized) up to
time $T$ spends a time $\Delta t_{n}\sim U[0,T]$ in the simulation. In that
case, the probability for success $p_{n,j}=\bar{p}_{j}\cdot\exp(-R_{i}\Delta
t_{n})$ becomes particle dependent, because each particle has to survive for a
different duration $\Delta t_{n}$, which means that they have a different
probability of being ionized before time $T$. As the probability distribution
for $\Delta t_{n}$ is known, we can calculate the following stationary
probability for success:

$\begin{split}p_{j}=\langle
p_{n,j}\rangle&=\int_{0}^{T}\bar{p}_{j}\cdot\exp(-R_{i}\delta)\cdot\mathcal{P}(\Delta
t_{n}=\delta)d\delta\\\
&=\int_{0}^{T}\bar{p}_{j}\cdot\exp(-R_{i}\delta)\cdot\frac{1}{T}d\delta\\\
&=\bar{p}_{j}\cdot\frac{1-\exp(-R_{i}T)}{R_{i}T}.\end{split}$ (53)

This stationary probability for success $p_{j}$ can be used in variance
expression (50) for simulations including sources and sinks.

###  4.2 Numerical experiments

We test the point estimator variance predictors for the three lowest order
moments (51) for a range of values for $R$ and $\sigma_{p}^{2}$. For each set
of parameters, we simulate $N=10^{2}$ particles and count one contribution to
each of the point estimators for each particle (Bernoulli trial) that ends up
in a bin $j$ at a time $T$, where $T$ is chosen such that $RT=10$. We repeat
the experiments once with $R_{i}=0$, and once with $R_{i}=1$. The test cases
are all executed on a 1D ($d=1$) domain $D=[0,1]$ with length 1 and periodic
boundary conditions. The particles are launched from a stationary particle
source, with an initial position distributed uniformly over the domain and a
velocity drawn from the post-collisional velocity distribution as given in
(7). The space is subdivided in $J=10$ bins, such that the probability for
success (53) equals

$p_{j}=\frac{1}{J}\cdot\frac{1-\exp(-R_{i}T)}{R_{i}T}=0.1\cdot\frac{1-\exp(-R_{i}T)}{R_{i}T}.$
(54)

The subdomain $D_{j}=[0,0.1]$ is considered as the bin of interest. The code
is openly available in a public repository [1].

We consider the same three experiments as in Section 3.5: the hydrodynamic
scaling experiment, the temperature scaling experiment, and the diffusive
scaling experiment. In each scaling, we estimate the variance on the 3 lowest
order moments: particle density, momentum, and energy. As a reference
solution, we determine for each set of parameters the empirical variance based
on $10^{4}$ realizations of the experiment. We plot the variance predictors
divided by $N$ for these three scalings in Figure 4 (with $R_{i}=0$) and
Figure 5 (with $R_{i}=1$), where we each time choose $u_{p}=0$, and use
$\sigma_{p}^{2}=1$ in the hydrodynamic scaling and $R=1$ in the temperature
scaling. In the plots where the variance remains constant for different values
of $R$ and $\sigma_{p}^{2}$, the statistical noise on the empirical variance
is visible. In each plot, the variance predictor corresponds well to the
empirical values for each QoI.

Figure 4: Analog particle tracing Monte Carlo with point estimator. Particles
are generated from a random source that is uniform in time, there is no
ionization sink ($R_{i}=0$). Top to bottom: variance on density, momentum and
energy. Left to right: hydrodynamic, temperature and diffusive scaling.

Figure 5: Analog particle tracing Monte Carlo with point estimator. Particles
are generated from a random source that is uniform in time, there is an
ionization sink with $R_{i}=1$. Top to bottom: variance on density, momentum
and energy. Left to right: hydrodynamic, temperature and diffusive scaling.

##  5 Predicting the variance of analog particle tracing Monte Carlo with
analog estimators

For analog estimators (16), the contributions (Bernoulli trials) delivered by
each particle are correlated. Using the knowledge that the correlation follows
from the positions $x_{n,k}$ (the current position $x_{n,k}$ depends on
previous positions $x_{n,m}$ with $0\leq m<k$), while the velocities $v_{n,k}$
are iid and position independent, the variance on an analog estimator can be
written in a similar form as (47) as follows:

$\mathbb{V}[\hat{Q}_{j}]=\frac{w^{2}}{(R\Delta
t)^{2}}\mathbb{V}[q_{x}(v_{n,k})]\langle\sum_{n=1}^{N}\sum_{k=0}^{K(n)}I_{j}(x_{n,k})\rangle+\frac{w^{2}}{(R\Delta
t)^{2}}\langle
q_{x}(v_{n,k})\rangle^{2}\mathbb{V}[\sum_{n=1}^{N}\sum_{k=0}^{K(n)}I_{j}(x_{n,k})].$
(55)

Note that $N$ and $K(n)$ are random variables, as the number of particles $N$
that reaches time $t_{1}$, the time at which the analog estimator starts
recording, is random, just as the number of collisions $K(n)$ during the time
interval of duration $\Delta t=t_{2}-t_{1}$. In the following, we first assume
that $N$ is known, such that only $K(n)$ is a random variable, in the case
that there is no stationary particle source. Next, we adjust the derivations
for the case where there is an active stationary particle source. Finally, we
let $N$ be a random variable as well, where we make the distinction between
having no stationary particle source and having an active stationary particle
source.

###  5.1 Known number of particles $N$

Let us first assume that the number of particles $N$ at time $t_{1}$ is known,
e.g., because $t_{1}=t_{0}=0$ such that the time interval $\Delta t$ starts at
the known initial condition or because there are no sources or sinks. Then
(55) becomes:

$\mathbb{V}[\hat{Q}_{j}]=\frac{w^{2}N}{(R\Delta
t)^{2}}\mathbb{V}[q_{x}(v_{n,k})]\langle\sum_{k=0}^{K(n)}I_{j}(x_{n,k})\rangle+\frac{w^{2}N}{(R\Delta
t)^{2}}\langle
q_{x}(v_{n,k})\rangle^{2}\mathbb{V}[\sum_{k=0}^{K(n)}I_{j}(x_{n,k})].$ (56)

Assuming that $N$ is known allows us to focus on the number of collisions
occurring during the time interval $\Delta t$. The collisions happen at
exponentially distributed times, so the number of collision events in time
interval $\Delta t$ is Poisson distributed. A collision is either a charge-
exchange collision or an ionization collision. If there is an ionization sink
($R_{i}>0$), then we have the following probability density for having
$K(n)=\kappa$ collisions:

$\begin{split}\mathcal{P}(K(n)=0)&=\exp(-R\Delta t),\\\
\mathcal{P}(K(n)=\kappa>0)&=\left(\frac{R_{cx}}{R}\right)^{\kappa-1}\cdot\frac{(R\Delta
t)^{\kappa}\exp(-R\Delta t)}{\kappa!}\\\
&+\sum_{l=\kappa+1}^{\infty}\left(\frac{R_{cx}}{R}\right)^{\kappa-1}\left(1-\frac{R_{cx}}{R}\right)\cdot\frac{(R\Delta
t)^{l}\exp(-R\Delta t)}{l!}.\end{split}$ (57)

For $K(n)=0$, we have the probability that no collision happens during time
interval $\Delta t$. For $K(n)>0$, the probability density (57) has two
contributions. First, a sequence of $\kappa$ events counts as a sequence of
$\kappa$ collisions if the first $\kappa-1$ events are charge-exchange
collisions. Second, a sequence with $l>\kappa$ events counts as a sequence of
$\kappa$ collisions if the first $\kappa-1$ events are charge-exchange
collisions and the $\kappa$-th event is an ionization collision. The expected
value and variance of distribution (57) are

$\begin{split}\langle K(n)\rangle&=\frac{R}{R_{i}}\left(1-\exp(-R_{i}\Delta
t)\right),\\\ \mathbb{V}[K(n)]&=\frac{2R^{2}-R_{i}R-(2R_{i}R_{cx}R\Delta
t+2R^{2}-R_{i}R)\exp(-R_{i}\Delta t)}{R_{i}^{2}}\\\
&-\left(\frac{R}{R_{i}}(1-\exp(-R_{i}\Delta t))\right)^{2}.\end{split}$ (58)

A derivation of these results is provided in Appendix A. Taking the limit
$R_{i}\rightarrow 0$ results in the standard Poisson distribution

$\lim_{R_{i}\rightarrow 0}\mathcal{P}(K(n)=\kappa)=\frac{(R\Delta
t)^{\kappa}\exp(-R\Delta t)}{\kappa!},$ (59)

with as expected value and variance

$\begin{split}\lim_{R_{i}\rightarrow 0}\langle K(n)\rangle&=R\Delta t,\\\
\lim_{R_{i}\rightarrow 0}\mathbb{V}[K(n)]&=R\Delta t.\end{split}$ (60)

Using the law of total expectation (see Appendix B), we can write the expected
value of a sum of indicator function evaluations with a random number of terms
$K(n)$ as follows:

$\langle{\sum_{k=0}^{K(n)}I_{j}(x_{n,k})}\rangle=\sum_{\kappa=0}^{\infty}\langle\sum_{k=0}^{\kappa}I_{j}(x_{n,k})\rangle\cdot\mathcal{P}(K(n)=\kappa).$
(61)

Using the law of total variance (see Appendix B), we can write the variance of
such a sum as

$\mathbb{V}[\sum_{k=0}^{K(n)}I_{j}(x_{n,k})]=\sum_{\kappa=0}^{\infty}\left(\mathbb{V}[\sum_{k=0}^{\kappa}I_{j}(x_{n,k})]+\left(\langle\sum_{k=0}^{\kappa}I_{j}(x_{n,k})\rangle-\langle\sum_{k=0}^{K(n)}I_{j}(x_{n,k})\rangle\right)^{2}\right)\cdot\mathcal{P}(K(n)=\kappa).$
(62)

Expression (62) is completely general and any of the expressions for binomial
experiments with a fixed number of trials $\kappa$, e.g., the upper bound
(23), the expression for independent trials (25), the Markov process based
expressions, or the hidden Markov process based expressions, can be inserted
to obtain an (approximate) expression for the variance on the binomial
experiment with a random number of trials.

As an example, we insert the Markov process based expression for large $K(n)$
(40) for the variance on the binomial experiment into (61) and (62) to obtain

$\begin{split}\langle{\sum_{k=0}^{K(n)}I_{j}(x_{n,k})}\rangle&\approx\sum_{\kappa=0}^{\infty}\kappa
p_{j}\cdot\mathcal{P}(K(n)=\kappa)=\langle K(n)\rangle p_{j},\\\
\mathbb{V}[\sum_{k=0}^{K(n)}I_{j}(x_{n,k})]&\approx\langle
K(n)\rangle\left(p_{j}(1-p_{j})+\frac{2p_{j}(1-p_{j})(\lambda_{j}-p_{j})}{1-\lambda_{j}}\right)+\mathbb{V}[K(n)]p_{j}^{2}.\end{split}$
(63)

Inserting these expressions in (56) results in:

$\begin{split}\mathbb{V}[\hat{Q}_{j}]&\approx\frac{w^{2}N}{(R\Delta
t)^{2}}\mathbb{V}[q_{x}(v_{n,k})]\langle K(n)\rangle p_{j}\\\
&+\frac{w^{2}N}{(R\Delta t)^{2}}\langle q_{x}(v_{n,k})\rangle^{2}\left(\langle
K(n)\rangle\left(p_{j}(1-p_{j})+\frac{2p_{j}(1-p_{j})(\lambda_{j}-p_{j})}{1-\lambda_{j}}\right)+\mathbb{V}[K(n)]p_{j}^{2}\right),\end{split}$
(64)

where we can substitute (48) to obtain variance predictors for the three
lowest order moments.

###  5.2 Treatment of a stationary particle source

When there is a stationary source as described in Section 2.1, the time that a
particle $n$ can effectively spend in the time interval $[t_{1},t_{2}]$ is a
uniformly distributed random variable $\Delta t_{n}\sim U[0,t_{2}-t_{1}]$, see
(13). This directly affects the probability distribution of the number of
collisions (Bernoulli trials) that the particle contributes to the analog
estimator. We can write the probability of having an effective time interval
of length $\delta$ and having $\kappa$ collisions in that time interval as
follows:

$\begin{split}\mathcal{P}(\Delta
t_{n}=\delta,K(n)=\kappa)&=\mathcal{P}(K(n)=\kappa\mid\Delta
t_{n}=\delta)\cdot\mathcal{P}(\Delta t_{n}=\delta)\\\
&=\mathcal{P}(K(n)=\kappa\mid\Delta
t_{n}=\delta)\cdot\frac{1}{(t_{2}-t_{1})},\end{split}$ (65)

where $\mathcal{P}(K(n)=\kappa\mid\Delta t_{n}=\delta)$ is given by (57) with
$\Delta t=\delta$. In variance expression (62), we need the marginal
probability distribution $\mathcal{P}(K(n)=\kappa)$. Integrating out the
dependence on $\delta$ in (65), we obtain the following expression for the
marginal probability distribution for having $\kappa$ collisions:

$\begin{split}\mathcal{P}(K(n)=\kappa)=&\int_{0}^{(t_{2}-t_{1})}\mathcal{P}(\Delta
t_{n}=\delta,K(n)=\kappa)d\delta\\\
=&\left(\frac{R_{cx}}{R}\right)^{\kappa-1}\Bigg{[}\frac{1}{R(t_{2}-t_{1})}-\frac{\Gamma(\kappa+1,R(t_{2}-t_{1}))}{R(t_{2}-t_{1})\cdot\kappa!}+\\\
&\left(1-\frac{R_{cx}}{R}\right)\left(1-\sum_{l=0}^{\kappa}\left(\frac{1}{R(t_{2}-t_{1})}-\frac{\Gamma(l+1,R(t_{2}-t_{1}))}{R(t_{2}-t_{1})\cdot
l!}\right)\right)\Bigg{]},\end{split}$ (66)

where we make use of incomplete Gamma functions (see Appendix C).
Additionally, we can derive the following expression for the expected number
of collisions:

$\begin{split}\langle
K(n)\rangle&=\int_{0}^{(t_{2}-t_{1})}\sum_{\kappa=0}^{\infty}\kappa\cdot\mathcal{P}(\Delta
t_{n}=\delta,K(n)=\kappa)d\delta\\\
&=\frac{R(R_{i}(t_{2}-t_{1})+\exp(-R_{i}(t_{2}-t_{1}))-1)}{R_{i}^{2}(t_{2}-t_{1})}.\end{split}$
(67)

Using l’Hôpital’s rule, we can verify that in the limit $R_{i}\rightarrow 0$,
we obtain

$\lim_{R_{i}\rightarrow 0}\langle K(n)\rangle=\frac{R(t_{2}-t_{1})}{2}.$ (68)

In this limit, the expected number of collisions is the expected number of
collisions for a fixed time interval of duration $\frac{(t_{2}-t_{1})}{2}$,
i.e., for the expected value of the effective time interval $\Delta t_{n}$.

The variance on the number of collisions is

$\begin{split}\mathbb{V}[K(n)]&=\int_{0}^{(t_{2}-t_{1})}\sum_{\kappa=0}^{\infty}(\kappa-\langle
K(n)\rangle)^{2}\cdot\mathcal{P}(\Delta t_{n}=\delta,K(n)=\kappa)d\delta\\\
&=\frac{2R^{2}\exp(-R_{i}(t_{2}-t_{1}))(R_{i}(t_{2}-t_{1})+\exp(R_{i}(t_{2}-t_{1}))(R_{i}(t_{2}-t_{1})-2)+2)}{R_{i}^{3}(t_{2}-t_{1})}\\\
&-\frac{RR_{i}\exp(-R_{i}(t_{2}-t_{1}))(2R_{i}(t_{2}-t_{1})+\exp(R_{i}(t_{2}-t_{1}))(R_{i}(t_{2}-t_{1})-3)+3)}{R_{i}^{3}(t_{2}-t_{1})}\\\
&-\left(\frac{R(R_{i}(t_{2}-t_{1})+\exp(-R_{i}(t_{2}-t_{1}))-1)}{R_{i}^{2}(t_{2}-t_{1})}\right)^{2}.\end{split}$
(69)

Again, using l’Hôpital’s rule, we can verify that in the limit
$R_{i}\rightarrow 0$ we obtain

$\lim_{R_{i}\rightarrow
0}\mathbb{V}[K(n)]=\frac{R(t_{2}-t_{1})}{2}+\frac{(R(t_{2}-t_{1}))^{2}}{12}.$
(70)

Note in this limit of the variance expression that
$\frac{(t_{2}-t_{1})^{2}}{12}$ is the variance of the uniform distribution
from which we draw $\Delta t_{n}$. So the new variance is the variance on the
number of collisions in the expected effective time interval (first term),
augmented with a term that scales with the variance on the effective time
interval duration.

To account for the stationary particle source, we can insert the marginal
probability distribution (66) in the variance expression (62). The expressions
for the mean and variance on $K(n)$ can, e.g., be used in the variance
predictor for the analog estimator based on the approximate Markov process
expression for large $K(n)$ (64).

###  5.3 Random number of particles $N$

Let us now consider the case of equation (55) where the number of particles
$N$ at time $t_{1}$, the time at which the analog estimator starts recording,
is uncertain. We denote the number of particles in the initial condition by
$N(t=t_{0}=0)=N_{0}$, the particles generated by a stationary particle source
by $N_{S}$, and the particles present at time $t_{1}$ by $N(t=t_{1})=N$.

If we know the probability distribution on $N$, then it follows from the law
of total expectation (see Appendix B) that

$\langle\sum_{n=1}^{N}\sum_{k=0}^{K(n)}I_{j}(x_{n,k})\rangle=\sum_{m=1}^{N_{i}}m\langle\sum_{k=0}^{K(n)}I_{j}(x_{n,k})\rangle\cdot\mathcal{P}(N=m)=\langle
N\rangle\langle\sum_{k=0}^{K(n)}I_{j}(x_{n,k})\rangle,$ (71)

where $N_{i}=N_{0}$ for particles originating from the initial condition and
$N_{i}=N_{S}$ for particles originating from the stationary particle source
that is active during the time interval $[t_{0},t_{1}]$. Additionally, we can
use the law of total variance (see Appendix B) in the following way:

$\begin{split}\mathbb{V}[\sum_{n=1}^{N}\sum_{k=0}^{K(n)}I_{j}(x_{n,k})]&=\mathbb{V}[\langle\sum_{n=1}^{N}\sum_{k=0}^{K(n)}I_{j}(x_{n,k})\mid
N\rangle]+\langle\mathbb{V}[\sum_{n=1}^{N}\sum_{k=0}^{K(n)}I_{j}(x_{n,k})\mid
N]\rangle\\\
&=\mathbb{V}[N\langle\sum_{k=0}^{K(n)}I_{j}(x_{n,k})\rangle]+\langle
N\mathbb{V}[\sum_{k=0}^{K(n)}I_{j}(x_{n,k})]\rangle\\\
&=\mathbb{V}[N]\langle\sum_{k=0}^{K(n)}I_{j}(x_{n,k})\rangle^{2}+\langle
N\rangle\mathbb{V}[\sum_{k=0}^{K(n)}I_{j}(x_{n,k})].\end{split}$ (72)

Expressions for the statistics of the sum over $k$ have been derived above in
sections 5.1 and 5.2. Therefore, the variance expressions (72) and (55) are
completely determined if we derive expressions for $\mathbb{V}[N]$ and
$\langle N\rangle$. The randomness in $N$ is introduced by the particle source
and ionization sink. We first consider the case where there is an ionization
sink, but no particle source. Then, we will look at the case where there is
also a stationary particle source.

Case 1: no stationary particle source
We are interested in the probability $\mathcal{P}(N=m)$, i.e., we want to know
how many particles there are at time $t_{1}$, the start of the estimation time
interval. Each of the $N_{0}$ particles at the start of the simulation ionizes
with probability $P_{ion}=1-\exp(-R_{i}t_{1})$ or survive with probability
$P_{survive}=1-P_{ion}=\exp(-R_{i}t_{1})$. If we then count all the different
realizations in which $m$ particles survive up to time $t_{1}$, we get the
following binomial probability distribution:

$\begin{split}\mathcal{P}(N=m)&=\frac{N_{0}!}{m!(N_{0}-m)!}P_{survive}^{m}P_{ion}^{N_{0}-m}\\\
&=\frac{N_{0}!}{m!(N_{0}-m)!}\exp(-R_{i}t_{1})^{m}(1-\exp(-R_{i}t_{1}))^{N_{0}-m}.\end{split}$
(73)

Recall that the binomial distribution describes binomial experiments with
independent trials. Using the properties of this distribution, we find the
following expressions for the expected number of particles $N$ and the
variance on $N$:

$\begin{split}\langle N\rangle&=N_{0}\exp(-R_{i}t_{1}),\\\
\mathbb{V}[N]&=N_{0}\exp(-R_{i}t_{1})(1-\exp(-R_{i}t_{1})).\end{split}$ (74)

Note that we can also find these expressions by defining a binomial experiment
$N=\sum_{n=1}^{N_{0}}I_{n}$, where this time $I_{n}$ represents the Bernoulli
trial of a particle surviving the time interval $[0,t_{1}]$ without being
ionized, which happens with success probability $p_{n}=p=\exp(-R_{i}t_{1})$,
and then considering the expected value and variance of this sum under the
assumption of independent trials. Also note that in the limit of
$R_{i}\rightarrow 0$ we obtain $\langle N\rangle=N_{0}$ and $\mathbb{V}[N]=0$,
which both are to be expected as the initial particles then are certain to
reach time $t_{1}$.

Case 2: active stationary particle source
When $N_{S}$ particles originate from a stationary particle source, the
survival probability becomes particle dependent: $p_{n}=\exp(-R_{i}\Delta
t_{n})$ with $\Delta t_{n}\sim\mathcal{P}(\Delta t_{n}=\delta)=U[0,t_{1}]$.
This time, obtaining an expression for $\mathcal{P}(N=m)$ is not trivial. We
therefore immediately write the binomial experiment
$N=\sum_{n=1}^{N_{S}}I_{n}$ with Bernoulli trials $I_{n}$ as explained in the
previous paragraph. The expected value of a binomial experiment with
independent trials, but trial dependent success probabilities $p_{n}$, is as
follows:

$\begin{split}\langle N\rangle&=\sum_{n=1}^{N_{S}}\langle
I_{n}\rangle=\sum_{n=1}^{N_{S}}\int_{0}^{t_{1}}(p_{n}\cdot 1+(1-p_{n})\cdot
0)\cdot\frac{1}{t_{1}}d\delta=N_{S}\int_{0}^{t_{1}}\exp(-R_{i}\delta)\cdot\frac{1}{t_{1}}d\delta\\\
&=N_{S}\cdot\frac{1-\exp(-R_{i}t_{1})}{R_{i}t_{1}}.\end{split}$ (75)

In the limit for $R_{i}\rightarrow 0$, this expression reduces using
l’Hôpital’s rule to $\langle N\rangle=N_{S}$. For the variance, we find

$\begin{split}\mathbb{V}[N]&=\mathbb{V}[\sum_{n=1}^{N_{S}}I_{n}]=\sum_{n=1}^{N_{S}}\langle
I_{n}\rangle(1-\langle I_{n}\rangle)\\\
&=N_{S}\cdot\left(\frac{1-\exp(-R_{i}t_{1})}{R_{i}t_{1}}\right)\left(1-\left(\frac{1-\exp(-R_{i}t_{1})}{R_{i}t_{1}}\right)\right).\end{split}$
(76)

In the limit for $R_{i}\rightarrow 0$, this expression reduces to a variance
of $0$, as expected because then we have $N\equiv N_{S}$.

###  5.4 Numerical experiments

We test the analog estimator variance predictors of the form (55) for the
three lowest order moments (48) for a range of values for $R$ and
$\sigma_{p}^{2}$. For each set of parameters, we simulate a single particle
and count one contribution to each of the analog estimators for each collision
point (Bernoulli trial) that is inside a bin $j$ during a time interval
$[t_{1},t_{2}]$. All the experiments use $t_{1}=0$ and $t_{2}=\frac{100}{R}$
such that the expected number of collisions during the time interval equals
100. We repeat the experiments two times: once with particles being launched
from an initial condition without sink ($R_{i}=0$); and once with particles
being launched from a stationary particle source with an active sink
($R_{i}=1$). The test cases are all executed on a 1D ($d=1$) domain $D=[0,1]$
with length 1 and periodic boundary conditions. The particles have an initial
position distributed uniformly over the domain and a velocity drawn from the
post-collisional velocity distribution as given in (7). The space is
subdivided in $J=10$ bins, such that the probability for success equals
$p_{j}=\frac{1}{J}=0.1$. The subdomain $D_{j}=[0,0.1]$ is considered as the
bin of interest. The code is openly available in a public repository [1].

We consider the same three experiments as in Section 3.5: the hydrodynamic
scaling experiment, the temperature scaling experiment, and the diffusive
scaling experiment. In each scaling, we estimate the variance on the 3 lowest
order moments: particle density, momentum, and energy. For the variance on the
binomial experiments, we insert the upper bound (23), the variance expression
for independent trials (25), the Markov process based variance expression
(36), and the discretized hidden Markov process based variance expression. The
hidden Markov process is discretized on an equidistant grid consisting of 100
cells. The conditional probability $\lambda_{j}$ for the Markov process and
transition probabilities for the discretized hidden Markov process are
calculated as a function of the model parameters as described in Ref. [20]. As
a reference solution, we determine for each set of parameters the empirical
variance based on $10^{4}$ realizations of the experiment. We plot the
different variance predictors divided by $Rt_{2}=100$ for these three scalings
in Figure 6 (no stationary particle source or ionization sink) and Figure 7
(active stationary particle source and ionization sink), where we each time
choose $u_{p}=0$, and use $\sigma_{p}^{2}=1$ in the hydrodynamic scaling and
$R=1$ in the temperature scaling. The upper bound for the variance is omitted
in the temperature scaling plots, to zoom in on the other variance
expressions. In all cases, the discretized hidden Markov process based
predictor corresponds well to the measured empirical values. The Markov
process based variance predictor achieves the correct order of magnitude,
which is valuable for error prediction where knowledge of the order of
magnitude typically is sufficient.

Figure 6: Analog particle tracing Monte Carlo with analog estimator. Particles
originate from an uniform initial condition, there is no ionization sink. Top
to bottom: variance on density, momentum and energy. Left to right:
hydrodynamic, temperature and diffusive scaling.

Figure 7: Analog particle tracing Monte Carlo with analog estimator. Particles
are generated from a random source that is uniform in time, there is an
ionization sink with $R_{i}=1$. Top to bottom: variance on density, momentum
and energy. Left to right: hydrodynamic, temperature and diffusive scaling.

##  6 Important remarks

The focus of this work is on predicting the variance on QoI estimators. For
point estimators, the expressions based on independent contributions can
easily be adopted. For analog estimators, the Markov process based expressions
constitute easy to compute a priori variance predictors. This section deals
with several aspects of statistical error prediction that need some more
elaboration. The first aspect is how to choose the success probability $p_{j}$
and the conditional probability $\lambda_{j}$ in the Markov process based
variance expressions. The second aspect is the performance of the different
variance expressions in higher dimensions. Finally, the variance on a QoI
estimator only gives a view of the absolute squared error. A better
representation of the accuracy of the QoI estimates is given by the relative
statistical error.

###  6.1 Choice of $p_{j}$ and $\lambda_{j}$

The probability for success $p_{j}$ scales with the particle density
$\rho_{j}$. E.g., for the density point estimator, we have that

$\rho_{j}(T)=\langle\hat{\rho}_{j}(T)\rangle=wNp_{j},$ (77)

meaning that $p_{j}=\frac{\rho_{j}(T)}{wN}$, i.e., the probability for success
depends on the solution itself. The solution, of course, is not available a
priori. To get a value for $p_{j}$, one can use a cheap (stationary)
approximative solution or simply assume that $p_{j}$ is uniform:
$p_{j}=\frac{1}{J}$.

The conditional probability $\lambda_{j}$, introduced in (26), depends on
local model parameters in cell $D_{j}$ and the size of that cell, but unlike
the probability for success $p_{j}$, it does not depend on the particle
density $\rho_{j}$. Ref. [20] describes how $\lambda_{j}$ can be calculated
when assuming an infinite domain without boundary conditions. As the
calculation only involves local model parameters, it can be used to obtain
local error predictors for each cell separately. Alternatively, one can only
consider the cell with the worst-case model parameters to predict an upper
bound on the variance. The effect of boundary conditions and non-homogeneities
between different cells in a real simulation still has to be investigated. For
boundary cells, it might be possible to take the boundary conditions into
account in the calculation of $\lambda_{j}$.

###  6.2 Additional experiments and higher dimensions

In addition to the experiments shown in this work, there are experiments with
different parameter settings and in higher dimensions (2D and 3D) available in
a public repository [1]. The higher dimensional experiments are
computationally expensive for the discretized hidden Markov process based
variance predictors. The Markov process based expressions remain cheap and
easily generalize to higher dimensions, because the generalization only
requires the calculation of the conditional probability $\lambda_{j}$ for a
higher dimensional bin of interest, which can be done as explained in Ref.
[20]. The variance predictors also perform satisfactorily in these higher
dimensional test cases.

###  6.3 Relative statistical errors

Using the predicted variances on QoI estimates $\mathbb{V}[\hat{Q}_{j}]$, we
can also predict relative statistical errors, which we define as follows:

$e_{j,x}=\frac{\sqrt{\mathbb{V}[\hat{Q}_{j}]}}{\langle\hat{Q}_{j}\rangle}=\sqrt{\frac{\mathbb{V}[\hat{Q}_{j}]}{\langle\hat{Q}_{j}\rangle^{2}}}.$
(78)

As an illustration, we will elaborate (78) for the point estimator expression
(50). Using (77), we obtain

$\begin{split}e_{j,x}&=\sqrt{\frac{w^{2}Np_{j}\mathbb{V}[q_{x}(v_{n}(T))]+w^{2}Np_{j}(1-p_{j})\langle
q_{x}(v_{n}(T))\rangle^{2}}{w^{2}\langle
q_{x}(v_{n}(T))\rangle^{2}N^{2}p_{j}^{2}}}=\sqrt{\frac{\mathbb{V}[q_{x}(v_{n}(T))]+(1-p_{j})\langle
q_{x}(v_{n}(T))\rangle^{2}}{\langle q_{x}(v_{n}(T))\rangle^{2}Np_{j}}}\\\
&=\sqrt{\frac{\mathbb{V}[q_{x}(v_{n}(T))]+(1-p_{j})\langle
q_{x}(v_{n}(T))\rangle^{2}}{\langle
q_{x}(v_{n}(T))\rangle^{2}}}\cdot\sqrt{\frac{w}{\rho_{j}(T)}}.\end{split}$
(79)

If we now, for example, write out the relative statistical error on the
density, we obtain

$e_{j,\rho}=\sqrt{\frac{1-p_{j}}{Np_{j}}}=\sqrt{1-p_{j}}\cdot\sqrt{\frac{w}{\rho_{j}(T)}}.$
(80)

This expression can also be found in Section 3.4.9 of the book [41] and a more
generalized version for the so-called delta-f method can be found in Ref. [8].
Note how the relative statistical errors scale with the square root of
$\frac{w}{\rho_{j}(T)}$. If the particle weight $w$ is small compared to the
density $\rho_{j}(T)$, then the relative statistical error will be small. If
the particle weight $w$ is large compared to the density $\rho_{j}(T)$, then
the relative statistical error is high. Increasing the number of particles
$N$, results in smaller particle weights $w$, resulting in smaller relative
statistical errors.

#### Remark 3

For the momentum estimator, we have that
$\mathbb{V}[q_{m}(v_{n}(T))]=\sigma_{p}^{2}$ and $\langle
q_{m}(v_{n}(T))\rangle=|u_{p}|$. If $|u_{p}|^{2}\ll\sigma_{p}^{2}$, then the
momentum estimator has a large relative statistical error (see (79)), making
it difficult to get an accurate QoI estimate.

##  7 Conclusions

In this work, variance predictors are devised for analog particle tracing
Monte Carlo methods that simulate linear kinetic equations. The expressions
are illustrated in 1D test cases but can easily be applied to higher
dimensional problems as well [20]. The variance predictors are built on the
idea that constructing a histogram can be interpreted as performing a
(correlated) binomial experiment. Using this insight, we first derived an
upper bound, an independent trials based, Markov process based, and hidden
Markov process based variance expressions for a (correlated) binomial
experiment. Next, these expressions were used to construct variance predictors
for different quantity of interest estimators. For point estimators, the
independent trials based estimators suffice. For analog estimators, the hidden
Markov process based expressions capture the variance exactly, but at a high
cost, making them unattractive for a priori variance prediction. The Markov
process based expressions are cheap to compute and typically capture the order
of magnitude of the variance correctly, which is satisfactory from an error
prediction point of view. The derived expressions give a prediction of the
variance in a given bin of the histogram, based on the bin size, the used
number of particles $N$, and the local model parameters. One can then choose
to either predict the variance for all bins, or to select the bin with worst-
case parameters and predict a worst-case upper bound for the variance over the
whole simulation domain. The resulting variance predictors can be used in
multiple ways: (1) to predict the relative statistical error, i.e., the
accuracy, on an estimate of a quantity of interest; (2) to characterize the
performance of a Monte Carlo method as a function of model parameters; (3) for
Monte Carlo method selection and optimization, also, e.g., if the Monte Carlo
method is a component of a more complex hybrid code.

In future work, the variance predictors can perhaps be extended to nonlinear
(or coupled) Monte Carlo methods such as DSMC and PIC, or other methods with
interacting particles, see, e.g., Ref. [6] and references therein. They can be
applied to different kinetic equations and other (nonanalog) simulation and
estimation procedures, enabling selection and optimization of nonanalog
particle tracing Monte Carlo methods. The error analysis can be used in the
design of hybrid simulation methods that contain a Monte Carlo component. The
theory has to be tested on real application cases to assess the influence of
the neglected $x,t$ dependencies of the model parameters. Finally, some work
can still be done on taking boundary conditions into account, e.g., in the
calculation of the conditional probability $\lambda_{j}$ in the Markov
dependence approximation.

## Acknowledgements

The authors would like to acknowledge the financial support of the CogniGron
research center and the Ubbo Emmius Funds (University of Groningen). The first
author is a SB PhD fellow of the Research Foundation Flanders (FWO), funded by
grant 1S64723N. The second author has received funding from the European High-
Performance Computing Joint Undertaking (JU) under grant agreement No 955701.
The JU receives support from the European Union’s Horizon 2020 research and
innovation programme, and Belgium, France, Germany, and Switzerland. Parts of
the work are supported by the Research Foundation Flanders (FWO) under grant
1179820N for fundamental research and project grant G085922N. This work has
been carried out within the framework of the EUROfusion Consortium, funded by
the European Union via the Euratom Research and Training Programme (Grant
Agreement No 101052200 — EUROfusion). Views and opinions expressed are,
however, those of the author(s) only and do not necessarily reflect those of
the European Union or the European Commission. Neither the European Union nor
the European Commission can be held responsible for them.

##  Appendix A Derivation of mean and variance of $\kappa$ events happening in
a time interval with sinks

In this appendix, we derive the mean and variance (58) of $K(n)$ with respect
to the probability distribution $\mathcal{P}(K(n)=\kappa)$ given by (57).

###  A.1 The mean

We start by inserting (57) in the definition of the mean:

$\begin{split}\langle
K(n)\rangle&=\sum_{\kappa=0}^{\infty}\kappa\cdot\mathcal{P}(K(n)=\kappa)\\\
&=\sum_{\kappa=1}^{\infty}\kappa\cdot\biggl{(}\left(\frac{R_{cx}}{R}\right)^{\kappa-1}\cdot\frac{(R\Delta
t)^{\kappa}\exp(-R\Delta t)}{\kappa!}\\\
&+\sum_{l=\kappa+1}^{\infty}\left(\frac{R_{cx}}{R}\right)^{\kappa-1}\left(1-\frac{R_{cx}}{R}\right)\cdot\frac{(R\Delta
t)^{l}\exp(-R\Delta t)}{l!}\biggr{)}.\\\ \end{split}$ (81)

We tackle this big equation in multiple parts. For the first term we get

$\begin{split}&\sum_{\kappa=1}^{\infty}\kappa\cdot\left(\frac{R_{cx}}{R}\right)^{\kappa-1}\cdot\frac{(R\Delta
t)^{\kappa}\exp(-R\Delta t)}{\kappa!}\\\ &=R\Delta t\exp(-R\Delta
t)\sum_{\kappa=1}^{\infty}\frac{(R_{cx}\Delta t)^{\kappa-1}}{(\kappa-1)!}\\\
&=R\Delta t\exp(-R\Delta t)\exp(R_{cx}\Delta t)\\\ &=R\Delta
t\exp(-R_{i}\Delta t).\end{split}$ (82)

For the remaining term of (81), we use (in order) property (100) and the
definition (96) of the incomplete Gamma function (see Appendix C), and the
Taylor expansion of the exponential function:

$\begin{split}&\sum_{\kappa=1}^{\infty}\kappa\cdot\sum_{l=\kappa+1}^{\infty}\left(\frac{R_{cx}}{R}\right)^{\kappa-1}\left(1-\frac{R_{cx}}{R}\right)\cdot\frac{(R\Delta
t)^{l}\exp(-R\Delta t)}{l!}\\\
&=\sum_{\kappa=1}^{\infty}\kappa\cdot\left(\frac{R_{cx}}{R}\right)^{\kappa-1}\left(1-\frac{R_{cx}}{R}\right)\cdot\exp(-R\Delta
t)\left(\exp(R\Delta t)-\sum_{l=0}^{\kappa}\frac{(R\Delta
t)^{l}}{l!}\right)\\\
&=\sum_{\kappa=1}^{\infty}\kappa\cdot\left(\frac{R_{cx}}{R}\right)^{\kappa-1}\left(1-\frac{R_{cx}}{R}\right)\cdot\left(1-\frac{\Gamma(\kappa+1,R\Delta
t)}{\kappa!}\right)\\\
&=\left(\sum_{\kappa=1}^{\infty}\kappa\cdot\left(\frac{R_{cx}}{R}\right)^{\kappa-1}\left(1-\frac{R_{cx}}{R}\right)\right)-\left(\sum_{\kappa=1}^{\infty}\kappa\cdot\left(\frac{R_{cx}}{R}\right)^{\kappa-1}\left(1-\frac{R_{cx}}{R}\right)\cdot\frac{\Gamma(\kappa+1,R\Delta
t)}{\kappa!}\right)\\\
&=\frac{R}{R_{i}}-\left(1-\frac{R_{cx}}{R}\right)\sum_{\kappa=1}^{\infty}\left(\frac{R_{cx}}{R}\right)^{\kappa-1}\cdot\frac{\Gamma(\kappa+1,R\Delta
t)}{(\kappa-1)!}\\\
&=\frac{R}{R_{i}}-\left(1-\frac{R_{cx}}{R}\right)\int_{R\Delta
t}^{\infty}t\exp(-t)\sum_{\kappa=1}^{\infty}\frac{\left(\frac{R_{cx}t}{R}\right)^{\kappa-1}}{(\kappa-1)!}dt\\\
&=\frac{R}{R_{i}}-\left(1-\frac{R_{cx}}{R}\right)\int_{R\Delta
t}^{\infty}t\exp(-t)\exp\left(\frac{R_{cx}t}{R}\right)dt\\\
&=\frac{R}{R_{i}}-\frac{R(R_{i}\Delta t+1)}{R_{i}}\exp(-R_{i}\Delta
t).\end{split}$ (83)

Combining the two parts, we finally obtain

$\begin{split}\langle K(n)\rangle&=R\Delta t\exp(-R_{i}\Delta
t)+\frac{R}{R_{i}}-\frac{R(R_{i}\Delta t+1)}{R_{i}}\exp(-R_{i}\Delta t)\\\
&=\frac{R}{R_{i}}(1-\exp(-R_{i}\Delta t)).\end{split}$ (84)

###  A.2 The variance

We start from the definition of the variance:

$\begin{split}\mathbb{V}[K(n)]&=\sum_{\kappa=0}^{\infty}(\kappa-\langle
K(n)\rangle)^{2}\cdot\mathcal{P}(K(n)=k)\\\
&=\left(\sum_{\kappa=1}^{\infty}\kappa^{2}\cdot\mathcal{P}(K(n)=\kappa)\right)-\langle
K(n)\rangle^{2}.\end{split}$ (85)

As the second term is known from above, we focus on the first term and insert
(57):

$\begin{split}&\sum_{\kappa=1}^{\infty}\kappa^{2}\cdot\mathcal{P}(K(n)=\kappa)\\\
&=\sum_{\kappa=1}^{\infty}\kappa^{2}\cdot\biggl{(}\left(\frac{R_{cx}}{R}\right)^{\kappa-1}\cdot\frac{(R\Delta
t)^{\kappa}\exp(-R\Delta t)}{\kappa!}\\\
&+\sum_{l=\kappa+1}^{\infty}\left(\frac{R_{cx}}{R}\right)^{\kappa-1}\left(1-\frac{R_{cx}}{R}\right)\cdot\frac{(R\Delta
t)^{l}\exp(-R\Delta t)}{l!}\biggr{)}.\end{split}$ (86)

We can work this out part by part as we did for the expected value of $K(n)$.
We provide some intermediate results. For the first part of the expression, we
obtain:

$\begin{split}&\sum_{\kappa=1}^{\infty}\kappa^{2}\cdot\left(\frac{R_{cx}}{R}\right)^{\kappa-1}\cdot\frac{(R\Delta
t)^{\kappa}\exp(-R\Delta t)}{\kappa!}\\\
&=\left(\sum_{\kappa=1}^{\infty}(\kappa-1)\left(\frac{R_{cx}}{R}\right)^{\kappa-1}\cdot\frac{(R\Delta
t)^{\kappa}\exp(-R\Delta
t)}{(\kappa-1)!}\right)+\left(\sum_{\kappa=1}^{\infty}\left(\frac{R_{cx}}{R}\right)^{\kappa-1}\frac{(R\Delta
t)^{\kappa}\exp(-R\Delta t)}{(\kappa-1)!}\right)\\\ &=(R\Delta
t)^{2}\left(\frac{R_{cx}}{R}\right)\exp(-R\Delta t)\exp(R_{cx}\Delta
t)+R\Delta t\exp(-R\Delta t)\exp(R_{cx}\Delta t)\\\ &=R\Delta t(1+R_{cx}\Delta
t)\exp(-R_{i}\Delta t).\end{split}$ (87)

Again using property (100) and the definition (96) of the incomplete Gamma
function, we can write the second part of the expression as follows:

$\begin{split}&\sum_{\kappa=1}^{\infty}\kappa^{2}\cdot\sum_{l=\kappa+1}^{\infty}\left(\frac{R_{cx}}{R}\right)^{\kappa-1}\left(1-\frac{R_{cx}}{R}\right)\cdot\frac{(R\Delta
t)^{l}\exp(-R\Delta t)}{l!}\\\
&=\sum_{\kappa=1}^{\infty}\kappa^{2}\left(\frac{R_{cx}}{R}\right)^{\kappa-1}\left(1-\frac{R_{cx}}{R}\right)\left(1-\frac{\Gamma(\kappa+1,R\Delta
t)}{\kappa!}\right)\\\
&=\frac{R(2R-R_{i})}{R_{i}^{2}}-\sum_{\kappa=1}^{\infty}\kappa^{2}\left(\frac{R_{cx}}{R}\right)^{\kappa-1}\left(1-\frac{R_{cx}}{R}\right)\frac{\Gamma(\kappa+1,R\Delta
t)}{\kappa!}\\\ &=\frac{R(2R-R_{i})}{R_{i}^{2}}\\\
&-\left(1-\frac{R_{cx}}{R}\right)\int_{R\Delta
t}^{\infty}\exp(-t)\left(t^{2}\left(\frac{R_{cx}}{R}\right)\exp\left(\left(1-\frac{R_{i}}{R}\right)t\right)+t\exp\left(\left(1-\frac{R_{i}}{R}\right)t\right)\right)dt\\\
&=\frac{R(2R-R_{i})}{R_{i}^{2}}\\\
&-\frac{R_{i}}{R}\left(\frac{R^{3}\exp(-R_{i}\Delta t)\left((R_{i}\Delta
t)^{2}+2R_{i}\Delta
t+2\right)}{R_{i}^{3}}\left(1-\frac{R_{i}}{R}\right)+\frac{R^{2}\exp(-R_{i}\Delta
t)(R_{i}\Delta t+1)}{R_{i}^{2}}\right).\end{split}$ (88)

Combining the two parts, we finally obtain

$\begin{split}&\sum_{\kappa=1}^{\infty}\kappa^{2}\cdot\mathcal{P}(K(n)=\kappa)=\frac{2R^{2}-R_{i}R-(2R_{i}R_{cx}R\Delta
t+2R^{2}-R_{i}R)\exp(-R_{i}\Delta t)}{R_{i}^{2}}.\end{split}$ (89)

Inserting this result in the expression for the variance (85), we get the
following result:

$\begin{split}\mathbb{V}[K(n)]&=\frac{2R^{2}-R_{i}R-(2R_{i}R_{cx}R\Delta
t+2R^{2}-R_{i}R)\exp(-R_{i}\Delta
t)}{R_{i}^{2}}-\left(\frac{R}{R_{i}}(1-\exp(-R_{i}\Delta
t))\right)^{2}.\end{split}$ (90)

##  Appendix B The law of total expectation and variance

The law of total expectation is as follows:

$\langle x\rangle=\langle\langle x\mid y\rangle\rangle.$ (91)

The inner expectation is over the random variable $x$ conditioned on $y$. The
outer expectation is over the random variable $y$. For example, if
$x=\sum_{l=1}^{L}z_{l}$, with a random number of terms $L$, and $y=L$, then we
have

$\begin{split}\langle\sum_{l=1}^{L}z_{l}\rangle&=\langle\langle\sum_{l=1}^{L}z_{l}\mid
L=\lambda\rangle\rangle\\\
&=\sum_{\lambda=0}^{\infty}\langle\sum_{l=1}^{L}z_{l}\mid
L=\lambda\rangle\cdot\mathcal{P}(L=\lambda)\\\
&=\sum_{\lambda=0}^{\infty}E(L=\lambda)\cdot\mathcal{P}(L=\lambda)\\\
&=E(E(L=\lambda)),\end{split}$ (92)

where we define $E(L)$ as the expected value of the sum $\sum_{l=1}^{L}z_{l}$
for a fixed value of $L$. This expression is used, e.g., in (61) and (71).

The law of total variance is as follows:

$\mathbb{V}[x]=\langle\mathbb{V}[x\mid y]\rangle+\mathbb{V}[\langle x\mid
y\rangle].$ (93)

In the first term, the variance is over $x$ conditioned on $y$ and the
expected value is over $y$. In the second term, the expected value is over $x$
conditioned on $y$ and the variance is over $y$. Considering again the example
with $x=\sum_{l=1}^{L}z_{l}$ and $y=L$, where for a fixed value of $L$ we have
$\mathbb{V}[\sum_{l=1}^{L}z_{l}]=\mathcal{V}(L)$. We find

$\begin{split}\mathbb{V}[\sum_{l=1}^{L}z_{l}]&=\langle\mathbb{V}[\sum_{l=1}^{L}z_{l}\mid
L=\lambda]\rangle+\mathbb{V}[\langle\sum_{l=1}^{L}z_{l}\mid
L=\lambda\rangle]\\\
&=\langle\mathcal{V}(L=\lambda)\rangle+\mathbb{V}[E(L=\lambda)]\\\
&=\sum_{\lambda=0}^{\infty}\mathcal{V}(L=\lambda)\cdot\mathcal{P}(L=\lambda)+\sum_{\lambda=0}^{\infty}(E(L=\lambda)-E(E(L=\lambda)))^{2}\cdot\mathcal{P}(L=\lambda),\end{split}$
(94)

where we defined $E(E(L=\lambda))$ in (92). This expression is used, e.g., in
(62) and (72).

##  Appendix C The incomplete Gamma function

The Gamma function is defined as

$\Gamma(z)=\int_{0}^{\infty}t^{z-1}\exp(-t)dt,$ (95)

with $\text{Re}(z)>0$. The incomplete Gamma function is defined as

$\Gamma(z,x)=\int_{x}^{\infty}t^{z-1}\exp(-t)dt,$ (96)

with $\text{Re}(z)>0$ and $x>0$.

For the Poisson distributions used throughout this text, we are interested in
the (incomplete) Gamma function where $z$ represents a number of collisions,
i.e., an integer $z\in\mathbb{N}$. In that case, the Gamma function can be
simplified using integration by parts as follows:

$\begin{split}\Gamma(z)&=\int_{0}^{\infty}t^{z-1}\exp(-t)dt\\\
&=(z-1)\int_{0}^{\infty}t^{z-2}\exp(-t)dt\\\
&=(z-1)(z-2)\int_{0}^{\infty}t^{z-3}\exp(-t)dt\\\ &\vdots\\\
&=(z-1)!\int_{0}^{\infty}\exp(-t)dt\\\ &=(z-1)!.\end{split}$ (97)

Analogously, we can simplify the incomplete Gamma function using integration
by parts as follows:

$\begin{split}\Gamma(z,x)&=\int_{x}^{\infty}t^{z-1}\exp(-t)dt\\\
&=x^{z-1}\exp(-x)+(z-1)\int_{x}^{\infty}t^{z-2}\exp(-t)dt\\\
&=x^{z-1}\exp(-x)+(z-1)x^{z-2}\exp(-x)+(z-1)(z-2)\int_{x}^{\infty}t^{z-3}\exp(-t)dt\\\
&\vdots\\\
&=\exp(-x)(z-1)!\sum_{k=1}^{z-1}\left(\frac{x^{k}}{k!}\right)+(z-1)!\int_{x}^{\infty}\exp(-t)dt\\\
&=\exp(-x)(z-1)!\sum_{k=1}^{z-1}\left(\frac{x^{k}}{k!}\right)+(z-1)!\exp(-x),\\\
&=\exp(-x)(z-1)!\left(\sum_{k=1}^{z-1}\left(\frac{x^{k}}{k!}\right)+1\right),\\\
\end{split}$ (98)

from which it follows that

$\sum_{k=1}^{z-1}\left(\frac{x^{k}}{k!}\right)=\frac{\exp(x)\Gamma(z,x)}{(z-1)!}-1=\frac{\exp(x)\Gamma(z,x)}{\Gamma(z)}-1$
(99)

and

$\sum_{k=0}^{z-1}\left(\frac{x^{k}}{k!}\right)=\frac{\exp(x)\Gamma(z,x)}{(z-1)!}=\frac{\exp(x)\Gamma(z,x)}{\Gamma(z)}.$
(100)

These expressions turn out to be useful when doing calculations with Poisson
distributions where we can then replace sums by Gamma functions (integrals)
that can be solved analytically, see, e.g., Appendix A.

## References

  * [1] Code accompanying the paper: ‘Predicting the statistical error of analog particle tracing Monte Carlo’.  https://gitlab.kuleuven.be/numa/public/statistical-errors-analog-ptmc.  Accessed: 2024-01-23.
  * [2] H. J. Amster and M. J. Djomehri.  Prediction of Statistical Error in Monte Carlo Transport Calculations.  Nuclear Science and Engineering, 60(2):131–142, June 1976.
  * [3] D. Andrieux.  Bounding the coarse graining error in hidden Markov dynamics.  Applied Mathematics Letters, 25(11):1734–1739, Nov. 2012.
  * [4] G. Bal, I. Langmore, and Y. Marzouk.  Bayesian inverse problems with Monte Carlo forward models.  Inverse Problems and Imaging, 7(1):81–105, 2013.
  * [5] R. Bellman and G. M. Wing.  An introduction to invariant imbedding.  Pure and applied mathematics. Wiley, New York, 1975.
  * [6] A. Bouillon, T. Ingelaere, and G. Samaey.  Single-ensemble multilevel Monte Carlo for discrete interacting-particle methods.  In preparation.
  * [7] I. D. Boyd and T. R. Deschenes.  Hybrid Particle-Continuum Numerical Methods for Aerospace Applications, 2011.
  * [8] S. Brunner, E. Valeo, and J. A. Krommes.  Collisional delta-f scheme with evolving background for transport time scale simulations.  Phys. Plasmas, 6(12):19, 1999.
  * [9] R. E. Caflisch.  Monte Carlo and quasi-Monte Carlo methods.  Acta Numerica, 7:1–49, Jan. 1998.
  * [10] G. Chen and I. D. Boyd.  Statistical Error Analysis for the Direct Simulation Monte Carlo Technique.  Journal of Computational Physics, 126(2):434–448, July 1996.
  * [11] D. R. Cox.  The Analysis of Multivariate Binary Data.  Applied Statistics, 21(2):113, 1972.
  * [12] E. G. Evstatiev, J. M. Finn, B. A. Shadwick, and N. Hengartner.  Noise and error analysis and optimization in particle-based kinetic plasma simulations.  Journal of Computational Physics, 440:110394, Sept. 2021.  arXiv:2102.02377 [physics].
  * [13] K. Ghoos, W. Dekeyser, G. Samaey, P. Börner, and M. Baelmans.  Accuracy and convergence of coupled finite-volume/Monte Carlo codes for plasma edge simulations of nuclear fusion reactors.  Journal of Computational Physics, 322:162–182, Oct. 2016.
  * [14] M. B. Giles.  Multilevel Monte Carlo methods.  Acta Numerica, 24:259–328, May 2015.
  * [15] J. Goodman and J. Weare.  Ensemble samplers with affine invariance.  Communications in Applied Mathematics and Computational Science, 5(1):65–80, Jan. 2010.
  * [16] C. Graham and D. Talay.  Stochastic Simulation and Monte Carlo Methods: Mathematical Foundations of Stochastic Simulation, volume 68.  01 2013.
  * [17] N. Horsten.  Fluid and Hybrid Fluid-Kinetic Models for the Neutral Particles in the Plasma Edge of Nuclear Fusion Devices.  Ph.D., KU Leuven, 2019.
  * [18] R. Indira.  Analytical study of leakage estimators in Monte-Carlo simulation of radiation transport.  Annals of Nuclear Energy, 15(5):261–269, Jan. 1988.
  * [19] R. Indira.  Optimization of biasing parameters in particle transport simulations.  Annals of Nuclear Energy, 16(11):571–576, Jan. 1989.
  * [20] T. Ingelaere, V. Maes, and G. Samaey.  Cell Escape Probabilities for Markov Processes on a Grid.  Submitted.
  * [21] M. A. Kalugin, D. S. Oleynik, and E. A. Sukhino-Khomenko.  Statistical error of reactor calculations by the Monte Carlo method.  Physics of Atomic Nuclei, 74(14):1871–1877, Dec. 2011.
  * [22] J. Klotz.  Statistical Inference in Bernoulli Trials with Dependence.  The Annals of Statistics, 1(2), Mar. 1973.
  * [23] E. Koehler, E. Brown, and S. J.-P. A. Haneuse.  On the Assessment of Monte Carlo Error in Simulation-Based Statistical Analyses.  The American Statistician, 63(2):155–162, May 2009.
  * [24] B. Lapeyre, E. Pardoux, and R. Sentis.  Introduction to Monte-Carlo Methods for Transport and Diffusion Equations.  Oxford University Press, 2003.
  * [25] K. Law, A. Stuart, and K. Zygalakis.  Data Assimilation: A Mathematical Introduction, volume 62 of Texts in Applied Mathematics.  Springer International Publishing, Cham, 2015.
  * [26] S. R. Lipsitz, G. M. Fitzmaurice, L. Sleeper, and L. P. Zhao.  Estimation Methods for the Joint Distribution of Repeated Binary Observations.  Biometrics, 51(2):562, June 1995.
  * [27] I. Lux.  Systematic Study of Some Standard Variance Reduction Techniques.  Nuclear Science and Engineering, 67(3):317–325, Sept. 1978.
  * [28] I. L. MacDonald.  A coarse-grained Markov chain is a hidden Markov model.  Physica A: Statistical Mechanics and its Applications, 541:123661, Mar. 2020.
  * [29] D. B. Macmillan.  Comparison of Statistical Estimators for Neutron Monte Carlo Calculations.  Nuclear Science and Engineering, 26(3):366–372, Nov. 1966.
  * [30] B. Mortier.  Advanced Monte Carlo simulation and estimation for kinetic neutral particles in the plasma edge of fusion reactors.  Ph.D., KU Leuven, 2020.
  * [31] B. Mortier, M. Baelmans, and G. Samaey.  A study of source term estimators in coupled finite-volume/Monte-Carlo methods with applications to plasma edge simulations in nuclear fusion: analog and collision estimators.  Accepted.
  * [32] B. Mortier, M. Baelmans, and G. Samaey.  Invariant imbedding applied to source term estimation procedures in a simplified 1d0d slab case.  Technical report, KU Leuven, 2020.
  * [33] R. S. Myong, A. Karchani, and O. Ejtehadi.  A review and perspective on a convergence analysis of the direct simulation Monte Carlo and solution verification.  Physics of Fluids, 31(6):066101, June 2019.
  * [34] J. F. Pendergast, S. J. Gange, M. A. Newton, M. J. Lindstrom, M. Palta, and M. R. Fisher.  A Survey of Methods for Analyzing Clustered Binary Response Data.  International Statistical Review / Revue Internationale de Statistique, 64(1):89, Apr. 1996.
  * [35] D. Reiter, M. Baelmans, and P. Börner.  The EIRENE and B2-EIRENE Codes.  Fusion Science and Technology, 47(2):172–186, Feb. 2005.
  * [36] D. Reiter, C. May, D. Coster, and R. Schneider.  Time dependent neutral gas transport in tokamak edge plasmas.  Journal of Nuclear Materials, 220-222:987–992, Apr. 1995.
  * [37] P. K. Sarkar and M. A. Prasad.  Prediction of Statistical Error and Optimization of Biased Monte Carlo Transport Calculations.  Nuclear Science and Engineering, 70(3):243–261, June 1979.
  * [38] D. W. Scott.  On optimal and data-based histograms.  Biometrika, 66(3):605–610, 1979.
  * [39] C. J. Solomon, A. Sood, and T. E. Booth.  Verification of the history-score moment equations for weight-window variance reduction.  Conference INIS-BR-17055, Brazil, 2011.
  * [40] J. Spanier.  An Analytic Approach to Variance Reduction.  SIAM Journal on Applied Mathematics, 18(1):172–190, Jan. 1970.
  * [41] J. Spanier and E. M. Gelbard.  Monte Carlo Principles and Neutron Transport Problems.  Addison-Wesley Publishing Company, 1969.
  * [42] M. Valentinuzzi, Y. Marandet, H. Bufferand, G. Ciraolo, and P. Tamain.  Two-phases hybrid model for neutrals.  Nuclear Materials and Energy, 18:41–45, Jan. 2019.
  * [43] U. Wolff.  Monte Carlo errors with less errors.  Computer Physics Communications, 176(5):383, Mar. 2007.  arXiv:hep-lat/0306017.
  * [44] E. Xekalaki and J. Panaretos.  A Binomial Distribution with Dependent Trials and Its Use in Stochastic Model Evaluation.  Communications in Statistics - Theory and Methods, 33(12):3043–3058, Jan. 2004.
  * [45] N. F. Zhang.  The Use of Correlated Binomial Distribution in Estimating Error Rates for Firearm Evidence Identification.  Journal of Research of the National Institute of Standards and Technology, 124:124026, Oct. 2019.