# On the Unbiasedness and Optimization of Regional Weight Cancellation

Hunter Belanger<EMAIL_ADDRESS>Davide Mancusi<EMAIL_ADDRESS>Andrea Zoia<EMAIL_ADDRESS>Université Paris-Saclay, CEA, Service d’Études
des Réacteurs et de Mathématiques Appliquées, 91191, Gif-sur-Yvette, France

###### Abstract

The Monte Carlo method is often used to simulate systems which can be modeled
by random walks. In order to calculate observables, in many implementations
the “walkers” carry a statistical weight which is generally assumed to be
positive. Some random walk simulations, however, may require walkers to have
positive or negative weights: it has been shown that the presence of a mixture
of positive and negative weights can impede the statistical convergence, and
special weight-cancellation techniques must be adopted in order to overcome
these issues. In a recent work we demonstrated the usefulness of one such
method, exact regional weight cancellation, to solve eigenvalue problems in
nuclear reactor physics in three spatial dimensions. The method previously
exhibited had several limitations (including multi-group transport and
isotropic scattering) and needed homogeneous cuboid cancellation regions. In
this paper we lift the previous limitations, in view of applying exact
regional cancellation to more realistic continuous-energy neutron transport
problems. This extended regional cancellation framework is used to optimize
the efficiency of the weight cancellation. Our findings are illustrated on a
benchmark configuration for reactor physics.

##  I Introduction

For day-to-day industrial needs in the field of nuclear reactor physics,
deterministic codes are used to solve the neutron transport equation,
estimating the reaction rates and the power distribution in the reactor core
[1, 2, 3]. Deterministic methods have the advantage of running fast, but this
speed comes at the cost of accuracy: many approximations are introduced,
discretizing the phase space (position, direction and energy) and thus leading
to a bias in the results. The gold-standard in reactor physics for solving the
neutron transport equation is the Monte Carlo method, which does not need to
introduce any discretization of the phase space, and is therefore free of any
bias [4]. This high-fidelity simulation method comes at the cost of requiring
extensive computer resources. Because of this computational cost, multi-
physics simulations of a full-core nuclear reactor model, using Monte Carlo
neutronics codes coupled with other state-of-the-art thermal-hydraulics and
fuel performance codes, have become possible only very recently, mainly thanks
to the increase in available computer power and to the development of
efficient variance-reduction techniques [5, 6].

In these Monte Carlo simulations, the particles being simulated (typically
neutrons or photons) carry a statistical weight, which is used to estimate
observable quantities such as reaction rates and power distributions within
the core of the nuclear reactor. For most applications involved in nuclear
reactor physics or radiation shielding problems, these statistical weights are
always positive. However, several types of Monte Carlo neutronics simulations
require that the particles being sampled carry negative statistical weights
(or complex weights, where each component is allowed to be negative). Such
types of calculations include the evaluation of the second harmonic of the
flux, critical buckling, neutron noise, and special rejection sampling methods
for spatially continuous material properties [7, 8, 9, 10, 11]. Random walk
problems using positive and negative statistical weights emerge more broadly
in many applications outside of nuclear reactor physics, e.g. in quantum
diffusion Monte Carlo [12]. Such simulations can be particularly challenging,
as the summing of positive and negative contributions to estimate the
observable quantities leads to very large variances in these tallies: it is
often recognized that weight cancellation is mandatory to ensure convergence
[8, 9, 13].

In a recent work, we have focused on the case of spatially continuous material
properties for particle transport applications emerging in reactor physics.
Material cross sections for neutron transport depend on the energy of the
incident particle, as well as on the temperature and density of the material.
Traditional neutronics codes (both Monte Carlo and deterministic) make the
approximation that each material region in the reactor model has a constant
temperature and density [14, 15, 16, 17, 1, 2, 3]. In a real nuclear reactor,
however, this is certainly not the case, as the temperature and density will
depend continuously on position. To improve the fidelity of our Monte Carlo
codes for the advent of the impending era of high-performance computing (HPC),
in a previous work we have examined the possibility of treating spatially-
continuous material temperatures and densities in Monte Carlo simulations
[11]. In particular, we have focused on assessing which particle-tracking
methods might be best suitable to treat spatially-continuous cross sections
for fixed-source transport problems, typically occurring in radiation
shielding applications [11]: among the possible choices, the delta-tracking
[18, 19] and negative-weighted delta-tracking [20, 21] sampling strategies
were deemed particularly attractive. Negative-weighted delta-tracking,
although beneficial for dealing with spatially-continuous cross sections, has
a potential drawback due to the statistical weights of the particles being
allowed to become negative: in a subsequent study concerning $k$-eigenvalue
problems, we have shown that the coupling of positive and negative particle
weights prevents convergence of the power iteration method to the fundamental
mode of the physical system being studied [22].

To overcome these issues, an exact regional weight cancellation method,
originally proposed by Booth and Gubernatis in a 1D context [23], was extended
to 3D and was shown to allow the convergence of power iteration with negative-
weighted delta-tracking in a multi-group reactor physics benchmark [22]. Such
a weight cancellation technique might be useful to improve the simulation
methods of the other previously mentioned problems which have particles with
negative statistical weights. While our previous work in Ref. 22 demonstrated
potential for the method, many questions were left unanswered: under what
conditions is regional cancellation unbiased? How might one extend
cancellation from multi-group to continuous-energy material cross sections? Is
it possible to maximize the efficiency of weight cancellation, for a given set
of particles in a cancellation region? Our goal in this paper is to build upon
our previous results in Ref. 22 and to start addressing these very questions.

Our manuscript is organised as follows. In Sec. II, we provide a brief summary
of the exact regional cancellation technique which we presented in our
previous paper. Section III will develop the the mathematical theory behind
the family of techniques for regional cancellation, for the most general
continuous-energy case. We also examine what conditions must be satisfied to
ensure that a regional cancellation method is unbiased. The theory presented
in this section elucidates the concepts which are integral to the technique
(for both multi-group and continuous-energy calculations), and blazes the
trail for implementing exact regional cancellation in continuous-energy
problems. The question of optimizing cancellation is then treated in Sec. IV,
where two candidate optimization methods are proposed. Section V discusses the
implementation of the two optimization strategies in our Monte Carlo code, and
discusses how these two strategies allow us to deal with heterogeneous
cancellation regions. The different optimization strategies are compared in
Sec. VI, and we also assess the performances of our methods on a reactor
physics benchmark with heterogeneous cancellation regions. Some concluding
thoughts and remarks are provided in Sec. VII.

##  II Regional Weight Cancellation

Previously, we have extended the 1D exact regional cancellation scheme of
Booth and Gubernatis [23] to work in 3D multi-group neutron transport problems
[22]. Here, we shall briefly outline the mechanics of this method, in a
general continuous-energy framework. For the case of $k$-eigenvalue problems,
the fundamental mode and eigenvalue are sought by Monte Carlo methods using
power iteration, which basically consists of following the neutron histories
over fission generations [4]. When negative-weighted delta-tracking is used to
sample particle flights, the transported neutron will have positive and
negative weights, and weight cancellation will be mandatory to ensure the
convergence of power iteration [22]. In this context, the regional
cancellation operation is applied to neutrons born from fission. The fission
particles are first sorted into user-defined cancellation regions, based on
their position. Once all of the particles have been sorted into their
cancellation regions, we may then consider each cancellation region
independently for the cancellation procedure. In our previous work, a simple
rectilinear mesh was imposed on top of the problem geometry.

Consider cancellation region $\mathcal{R}$ (which is assumed to be composed of
only one fissile material), containing fission neutrons which have already
been sampled. In addition to storing its own position ($\bm{r}$), energy
($E$), and direction ($\bm{\hat{\Omega}}$), each fission particle also stores
its parent’s energy ($E^{\prime}$), the position of the previous collision
($\bm{r}^{\prime}$), and the direction of the parent’s penultimate flight
($\bm{\hat{\Omega}}^{\prime\prime}$)111The direction of the parent’s last
flight is not explicitly stored, as it can be calculated as
$\bm{\hat{\Omega}}^{\prime}=\frac{\bm{r}-\bm{r}^{\prime}}{\absolutevalue{\bm{r}-\bm{r}^{\prime}}}$..
From this information, we can calculate the “fission density function”, i.e.
the expected fission density at $\bm{r}$ due to a collision at
$\bm{r}^{\prime}$ coming from direction $\bm{\hat{\Omega}}^{\prime\prime}$,
and a subsequent flight from $\bm{r}^{\prime}$ to $\bm{r}$ at energy
$E^{\prime}$; this is a key ingredient for the weight cancellation procedure.
The exact form the fission density function depends on the particle tracking
method being used. For the case of negative-weighted delta-tracking222We will
only mention negative-weighted delta-tracking in the text, since that was the
focus of our previous work, but Eq. (1) is also valid for regular delta-
tracking, where $\Sigma_{\text{smp}}$ would be the majorant cross section.
This could be of use for neutron noise or critical buckling problems, which
would not necessarily require the use of negative-weighted delta-tracking but
nonetheless require weight cancellation., as examined in our previous work,
the fission density function was taken to be

$\zeta(\bm{r}|\bm{r}^{\prime},\bm{\hat{\Omega}}^{\prime},E^{\prime})=\\\
\frac{\mathcal{P}\left(\frac{\bm{r}-\bm{r}^{\prime}}{\absolutevalue{\bm{r}-\bm{r}^{\prime}}}\cdot\bm{\hat{\Omega}}^{\prime}\right)\Sigma_{f}(\bm{r},E^{\prime}|\bm{r}^{\prime},E^{\prime})}{2\pi\absolutevalue{\bm{r}-\bm{r}^{\prime}}^{2}}e^{-\Sigma_{\text{smp}}(E^{\prime})\absolutevalue{\bm{r}-\bm{r}^{\prime}}}\text{.}$
(1)

In this notation, $\Sigma_{\text{smp}}$ is the sampling cross-section required
for negative-weighted delta tracking, $\Sigma_{f}$ is the fission cross-
section, and $\mathcal{P}$ is the probability density function for the cosine
of the scattering angle for the previous collision.333While the symbol $f$ was
used for the fission density function in Ref. 22, we have instead chosen to
use $\zeta$ in this paper, to avoid any confusion with other subsequent
symbols. Based on $\zeta$, we are able to split each fission particle in
$\mathcal{R}$ into two components: a point-wise component with weight $w_{p}$,
and a uniform component of weight $w_{u}$. The point-wise portion, $w_{p}$,
keeps the phase space coordinates $(\bm{r},\bm{\hat{\Omega}},E)$ of the split
fission particle. The uniform component, $w_{u}$, is spread uniformly over the
region $\mathcal{R}$. To calculate the point-wise and uniform weights, we use

$\displaystyle
w_{p}=\frac{\zeta(\bm{r}|\bm{r}^{\prime},\bm{\hat{\Omega}}^{\prime},E^{\prime})-\beta}{\zeta(\bm{r}|\bm{r}^{\prime},\bm{\hat{\Omega}}^{\prime},E^{\prime})}w$
(2) $\displaystyle
w_{u}=\frac{\beta}{\zeta(\bm{r}|\bm{r}^{\prime},\bm{\hat{\Omega}}^{\prime},E^{\prime})}w\text{,}$
(3)

respectively, with $w$ being the weight of the original fission particle [23].
Note that $w_{p}+w_{u}=w$, so that the net weight is conserved. The free
parameter $\beta$ can take any value, and in general is chosen independently
for each particle in $\mathcal{R}$. Our previous work followed the
recommendation of Booth and Gubernatis, and always took $\beta$ to be the
minimum value of
$\zeta(\bm{r}^{\prime\prime}|\bm{r}^{\prime},\bm{\hat{\Omega}}^{\prime},E^{\prime})$
over all possible $\bm{r}^{\prime\prime}\in\mathcal{R}$, for the particle of
interest. We demonstrated that, for the case of isotropic scattering and
cuboid cancellation regions, one only needs to evaluate
$\zeta(\bm{r}^{\prime\prime}|\bm{r}^{\prime},\bm{\hat{\Omega}}^{\prime},E^{\prime})$
for the eight corners of the cuboid (taking $\bm{r}^{\prime\prime}$ to be the
corner positions) to find the minimum value within $\mathcal{R}$.

With $w_{p}$ and $w_{u}$ having been calculated for each fission particle in
$\mathcal{R}$, we then take the sum of all the uniform weight components

$U=\sum_{i=1}^{N}w_{u,i}\text{,}$ (4)

where the extra subscript $i$ indicates the fission particle. This operation
is effectively where the cancellation occurs: depending on the initial weights
$w_{i}$ of the fission particles, the individual uniform components $w_{u,i}$
will be positive or negative, and taking their sum cancels some of the
positive and negative weight which was in the region $\mathcal{R}$. The
uniform weight $U$ must be distributed uniformly within $\mathcal{R}$. To do
this, $n=\lceil\absolutevalue{U}\rceil$ new fission particles are sampled
within $\mathcal{R}$, each having a weight of $U/n$. The positions of the $n$
uniform particles are sampled uniformly in $\mathcal{R}$. In our previous
work, the direction was sampled from an isotropic distribution, as fission was
assumed to be perfectly isotropic, and the energy was sampled from the fission
spectrum of the material in $\mathcal{R}$, as it was assumed that the fission
spectrum had no dependence on incident neutron energy. These $n$ new uniform
fission particles must be added to the fission bank, and will then be
transported along with the other fission particles during the next fission
generation.

The method proposed in Ref. 22 that we have recalled here has several
limitations. First, cancellation regions must be homogeneous, containing only
a single fissile material. Second, fission must always be isotropic, and the
fission spectrum must be independent of the incident energy. In general, even
in continuous-energy transport, fission is almost always represented as
isotropic, so this is not necessarily a large inconvenience. However, the
fission energy spectrum is generally assumed to be dependent on the incident
neutron energy. Furthermore, while Booth and Gubernatis argue that
cancellation is unbiased for any value of the parameter $\beta$, the amount of
canceled weight (and thus the efficiency of the method) clearly does depend on
$\beta$. Using the minimum value of the fission density as $\beta$ is not
necessarily the most efficient choice for achieving the highest amount of
cancellation. Nonetheless, taking $\beta$ to be the minimum within
$\mathcal{R}$ guarantees that both $w_{p}$ and $w_{u}$ have the same sign as
$w$: when $\beta$ is larger than the minimum, the point-wise portion, $w_{p}$,
can change sign, potentially leading to even more positive and negative weight
in the region than there was initially. The cancellation operation does not
change the net weight $W_{\text{net}}$ in the bin, as

$W_{\text{net}}=\sum_{i}^{N}w_{i}=\sum_{i}^{N}w_{u,i}+\sum_{i}^{N}w_{p,i}=U+\sum_{i}^{N}w_{p,i}$
(5)

will still be located in the bin. However, cancellation does change the total
weight, $W_{\text{tot}}$, defined as the sum of the absolute values of all
weights. The total weight before cancellation is

$W_{\text{tot}}=\sum_{i}^{N}\absolutevalue{w_{i}}\text{,}$ (6)

while the post-cancellation total weight is

$W_{\text{tot,post}}=\sum_{i}^{N}\absolutevalue{w_{p,i}}+\absolutevalue{U}=\sum_{i}^{N}\absolutevalue{w_{p,i}}+\absolutevalue{\sum_{i}^{N}w_{u,i}}\text{.}$
(7)

By using the triangle inequality, it is possible to show that

$W_{\text{tot,post}}\geq|W_{\text{net}}|\text{.}$

The more efficient cancellation is, the closer $W_{\text{tot,post}}$ will be
to $W_{\text{net}}$, with 100% cancellation efficiency corresponding to
$W_{\text{tot,post}}=W_{\text{net}}$ (i.e. all negative weight is removed).
The optimal choice for $\beta$ will maximize the cancellation efficiency, and
therefore minimize $W_{\text{tot,post}}$. This optimal choice of $\beta$ is
clearly dependent on the other particles in the bin, and is hitherto unknown.

##  III Unbiasedness of Cancellation

In this section, several different methods for performing exact regional
cancellation shall be examined, and it will be demonstrated under what
conditions such schemes lead to an unbiased fission source. For this purpose,
it is mathematically beneficial to use the integral form of the transport
equation, as opposed to the integro-differential form adopted in our previous
work. We will begin by presenting the integral transport form for the
eigenvalue transport problem in Sec. III.1. Section III.2 is divided into
several portions, each with a particular goal. Section III.2.1 make a first
attempt at developing an estimator for the fission emission density in a
region, which averages over all possible collisions and subsequent flights
which induce the fission. While this exact estimator is likely of little use
to a practical application, we are able to use it to examine what requirements
must be observed in order to have an unbiased fission emission density
estimator. Section III.2.2 discusses how far back in a particle’s history one
must look, so that it is possible for it to have contributed to the fission
emission density everywhere within the cancellation region. Section III.2.3
uses the ideas from Section III.2.2 to decompose the collision operator as is
done in most Monte Carlo codes, to produce a fission emission density
estimator which could potentially be used in an industrial code to achieve
exact regional cancellation. Section III.2.4 outlines the possibility of
distributing some of the the fission emission density within the region
according to a generic function, instead of distributing it uniformly.
Finally, Section III.2.5 examines why delta-tracking algorithms are more
suited to exact regional cancellation, and the peculiarities which can arise
from delta-scatters.

###  III.1 Integral Formulation of the Transport Equation

We will start with the $k$-eigenvalue Boltzmann transport equation in integral
form. Let $P=(\bm{r},\bm{\hat{\Omega}},E)$ denote the coordinates of a point
in phase space. The collision density
$\psi(P)=\Sigma_{t}(\bm{r},E)\varphi(\bm{r},\bm{\hat{\Omega}},E)$ and the
emission density $\chi(P)$ are related by [4]:

$\displaystyle\psi(P)$ $\displaystyle{}=\mathbb{T}\chi(P)$ (8)
$\displaystyle\chi(P)$
$\displaystyle{}=\left[\mathbb{C}_{s}+\frac{1}{k}\mathbb{C}_{f}\right]\psi(P)\text{.}$
(9)

In this notation, $\mathbb{T}$ is the flight operator, defined as

$\mathbb{T}g(P)=\int T(P^{\prime}\to P)g(P^{\prime})\differential
P^{\prime}\text{,}$ (10)

where we have made use of the flight kernel

$T(P^{\prime}\to P)=\\\
\frac{\Sigma_{t}(\bm{r},E^{\prime})}{\absolutevalue{\bm{r}-\bm{r}^{\prime}}^{2}}\exp\left(-\int\limits_{0}^{\absolutevalue{\bm{r}-\bm{r}^{\prime}}}\Sigma_{t}(\bm{r}^{\prime}+u\bm{\hat{\Omega}}^{\prime},E^{\prime})\differential
u\right)\\\
\delta\left(\bm{\hat{\Omega}}-\frac{\bm{r}-\bm{r}^{\prime}}{\absolutevalue{\bm{r}-\bm{r}^{\prime}}}\right)\delta\left(\bm{\hat{\Omega}}^{\prime}-\bm{\hat{\Omega}}\right)\delta\left(E-E^{\prime}\right)\text{.}$
(11)

We note that the flight kernel $T(P^{\prime}\to P)$ is normalized, and can be
interpreted as the probability density function (PDF) for a particle having a
flight and landing at the phase space coordinate $P$, conditioning on its
initial phase space coordinate being $P^{\prime}$.

The scattering operator $\mathbb{C}_{s}$ in Eq. (9) is defined as

$\mathbb{C}_{s}g(P)=\int C_{s}(P^{\prime}\to P)g(P^{\prime})\differential
P^{\prime}\text{,}$ (12)

with the scattering kernel $C_{s}(P^{\prime}\to P)$ being

$C_{s}(P^{\prime}\to
P)=\frac{\nu_{s}(\bm{r}^{\prime},E^{\prime})\Sigma_{s}(\bm{r}^{\prime},E^{\prime})}{\Sigma_{t}(\bm{r}^{\prime},E^{\prime})}\\\
f_{\text{scat}}\left(\bm{\hat{\Omega}},E|\bm{r}^{\prime},\bm{\hat{\Omega}}^{\prime},E^{\prime}\right)\delta\left(\bm{r}-\bm{r}^{\prime}\right)\text{.}$
(13)

The fission operator $\mathbb{C}_{f}$ is similar to the scattering operator in
Eq. (12), but instead uses a fission kernel $C_{f}(P^{\prime}\to P)$, defined
as

$C_{f}(P^{\prime}\to
P)=\frac{\nu_{f}(\bm{r}^{\prime},E^{\prime})\Sigma_{f}(\bm{r}^{\prime},E^{\prime})}{\Sigma_{t}(\bm{r}^{\prime},E^{\prime})}\\\
f_{\text{fiss}}\left(\bm{\hat{\Omega}},E|\bm{r}^{\prime},\bm{\hat{\Omega}}^{\prime},E^{\prime}\right)\delta\left(\bm{r}-\bm{r}^{\prime}\right)\text{.}$
(14)

The scattering and fission operators may be combined into a collision operator

$\chi(P)=\mathbb{C}\psi(P)\text{,}$ (15)

which has a corresponding collision kernel

$C(P^{\prime}\to P)=C_{s}(P^{\prime}\to P)+\frac{1}{k}C_{f}(P^{\prime}\to
P)\text{.}$ (16)

Here $C(P^{\prime}\to P)$ can be interpreted as the average number of of
particles produced about the phase space coordinate $P$, from a collision
induced by a particle at $P^{\prime}$. Given this interpretation, it is also
possible to rewrite $C(P^{\prime}\to P)$ in a more concise form, using an
average yield $\bar{\nu}(\bm{r}^{\prime},E^{\prime})$, and an average transfer
function
$\bar{f}(\bm{\hat{\Omega}},E|\bm{r}^{\prime},\bm{\hat{\Omega}}^{\prime},E^{\prime})$:

$C(P^{\prime}\to
P)=\bar{\nu}(\bm{r}^{\prime},E^{\prime})\bar{f}\left(\bm{\hat{\Omega}},E|\bm{r}^{\prime},\bm{\hat{\Omega}}^{\prime},E^{\prime}\right)\delta(\bm{r}-\bm{r}^{\prime})\text{.}$
(17)

It is clear that Eq. (17) is true if

$\bar{\nu}(\bm{r}^{\prime},E^{\prime})=\frac{\nu_{s}(\bm{r}^{\prime},E^{\prime})\Sigma_{s}(\bm{r}^{\prime},E^{\prime})}{\Sigma_{t}(\bm{r}^{\prime},E^{\prime})}+\frac{\nu_{f}(\bm{r}^{\prime},E^{\prime})\Sigma_{f}(\bm{r}^{\prime},E^{\prime})}{k\Sigma_{t}(\bm{r}^{\prime},E^{\prime})}$
(18)

and

$\bar{f}\left(\bm{\hat{\Omega}},E|\bm{r}^{\prime},\bm{\hat{\Omega}}^{\prime},E^{\prime}\right)=\\\
\frac{\nu_{s}(\bm{r}^{\prime},E^{\prime})\Sigma_{s}(\bm{r}^{\prime},E^{\prime})}{\bar{\nu}(\bm{r}^{\prime},E^{\prime})\Sigma_{t}(\bm{r}^{\prime},E^{\prime})}f_{\text{scat}}\left(\bm{\hat{\Omega}},E|\bm{r}^{\prime},\bm{\hat{\Omega}}^{\prime},E^{\prime}\right)+\\\
\frac{\nu_{f}(\bm{r}^{\prime},E^{\prime})\Sigma_{f}(\bm{r}^{\prime},E^{\prime})}{k\bar{\nu}(\bm{r}^{\prime},E^{\prime})\Sigma_{t}(\bm{r}^{\prime},E^{\prime})}f_{\text{fiss}}\left(\bm{\hat{\Omega}},E|\bm{r}^{\prime},\bm{\hat{\Omega}}^{\prime},E^{\prime}\right)\text{.}$
(19)

###  III.2 Unbiased Estimators of the Fission Emission Density

####  III.2.1 Averaging over all Scattering Events

Figure 1: Depicted here is the relationship between phase space points
$P_{1}$, $P_{2}$, $P_{3}$, and $P_{4}$. Points connected by a flight operator
(solid line) share the same direction and energy, and are only discontinuous
in position. Points connected by a collision or fission operator (dotted and
dashed lines respectively) share the same position, but are generally
discontinuous in direction and energy. Any branches which might result from
the application of the collision operator from $P_{1}$ to $P_{2}$ are not
depicted.

Consider the following particle history. A neutron enters a collision at
$P_{1}$, and then leaves that collision at $P_{2}$. The particle then
undergoes a flight and experiences a fission at $P_{3}$. The fission at
$P_{3}$ then contributes to the fission emission density at $P_{4}$. It is
assumed that $P_{4}$ is located within the generalized phase space region
$\mathcal{R}$. This partial particle history is depicted in Fig. 1. Despite
the fact that $\bm{r}_{4}=\bm{r}_{3}$, $P_{3}$ is not, in general, located in
the cancellation region $\mathcal{R}$, as $\bm{\hat{\Omega}}_{3}$ and $E_{3}$
may not be within the domain of $\mathcal{R}$. In order to examine the fission
emission density at point $P_{4}$, we must first determine the collision
density $\psi(P_{3})$, for a given collision at $P_{1}$. From Eq. (8) and Eq.
(9), it follows that

$\psi(P_{3})=\mathbb{T}\mathbb{C}\psi(P_{3})=\int T(P_{2}\to P_{3})\int
C(P_{1}\to P_{2})\psi(P_{1})\differential P_{1}\differential P_{2}$ (20)
$=\iiint\differential\bm{r}_{1}\differential\bm{\hat{\Omega}}_{1}\differential
E_{1}\iiint\differential\bm{r}_{2}\differential\bm{\hat{\Omega}}_{2}\differential
E_{2}\psi\left(\bm{r}_{1},\bm{\hat{\Omega}}_{1},E_{1}\right)\bar{\nu}(\bm{r}_{1},E_{1})\bar{f}\left(\bm{\hat{\Omega}}_{2},E_{2}|\bm{r}_{1},\bm{\hat{\Omega}}_{1},E_{1}\right)\delta\left(\bm{r}_{2}-\bm{r}_{1}\right)\\\
\Sigma_{t}\left(\bm{r}_{3},E_{3}\right)\exp\left(-\int\limits_{0}^{\absolutevalue{\bm{r}_{3}-\bm{r}_{2}}}\Sigma_{t}\left(\bm{r}_{2}+u\bm{\hat{\Omega}}_{2},E_{2}\right)\differential
u\right)\frac{\delta\left(\bm{\hat{\Omega}}_{3}-\frac{\bm{r}_{3}-\bm{r}_{2}}{\absolutevalue{\bm{r}_{3}-\bm{r}_{2}}}\right)\delta\left(\bm{\hat{\Omega}}_{3}-\bm{\hat{\Omega}}_{2}\right)\delta(E_{3}-E_{2})}{\absolutevalue{\bm{r}_{3}-\bm{r}_{2}}^{2}}$
(21)
$=\iiint\differential\bm{r}_{1}\differential\bm{\hat{\Omega}}_{1}\differential
E_{1}\psi\left(\bm{r}_{1},\bm{\hat{\Omega}}_{1},E_{1}\right)\bar{\nu}(\bm{r}_{1},E_{1})\\\
\delta\left(\bm{\hat{\Omega}}_{3}-\frac{\bm{r}_{3}-\bm{r}_{1}}{\absolutevalue{\bm{r}_{3}-\bm{r}_{1}}}\right)\frac{\bar{f}\left(\bm{\hat{\Omega}}_{3},E_{3}|\bm{r}_{1},\bm{\hat{\Omega}}_{1},E_{1}\right)\Sigma_{t}(\bm{r}_{3},E_{3})}{\absolutevalue{\bm{r}_{3}-\bm{r}_{1}}^{2}}\exp\left(-\int\limits_{0}^{\absolutevalue{\bm{r}_{3}-\bm{r}_{1}}}\Sigma_{t}\left(\bm{r}_{1}+u\bm{\hat{\Omega}}_{3},E_{3}\right)\differential
u\right)\text{.}$ (22)

The integral over $P_{1}$ in Eq. (22) indicates that $\psi(P_{3})$ is a sum of
contributions from all possible initial phase space points $P_{1}$ for which
$\bm{\hat{\Omega}}_{3}=\frac{\bm{r}_{3}-\bm{r}_{1}}{\absolutevalue{\bm{r}_{3}-\bm{r}_{1}}}$.

The fission emission density $\chi_{f}(P_{4})$ is defined as

$\chi_{f}(P_{4})=\frac{1}{k}\mathbb{C}_{f}\psi(P_{4})\text{.}$ (23)

Combining Eq. (14) and Eq. (22), we obtain

$\chi_{f}(P_{4})=\iiint\differential\bm{r}_{1}\differential\bm{\hat{\Omega}}_{1}\differential
E_{1}\psi\left(\bm{r}_{1},\bm{\hat{\Omega}}_{1},E_{1}\right)\bar{\nu}(\bm{r}_{1},E_{1})\\\
\int\differential
E_{3}\frac{\nu_{f}(\bm{r}_{4},E_{3})\Sigma_{f}(\bm{r}_{4},E_{3})f_{\text{fiss}}\left(\bm{\hat{\Omega}}_{4},E_{4}|\bm{r}_{4},\frac{\bm{r}_{4}-\bm{r}_{1}}{\absolutevalue{\bm{r}_{4}-\bm{r}_{1}}},E_{3}\right)\bar{f}\left(\frac{\bm{r}_{4}-\bm{r}_{1}}{\absolutevalue{\bm{r}_{4}-\bm{r}_{1}}},E_{3}|\bm{r}_{1},\bm{\hat{\Omega}}_{1},E_{1}\right)}{k\absolutevalue{\bm{r}_{4}-\bm{r}_{1}}^{2}}\\\
\exp\left(-\int\limits_{0}^{\absolutevalue{\bm{r}_{4}-\bm{r}_{1}}}\Sigma_{t}\left(\bm{r}_{1}+u\frac{\bm{r}_{4}-\bm{r}_{1}}{\absolutevalue{\bm{r}_{4}-\bm{r}_{1}}},E_{3}\right)\differential
u\right)\\\
=\iiint\psi\left(\bm{r}_{1},\bm{\hat{\Omega}}_{1},E_{1}\right)\bar{\nu}(\bm{r}_{1},E_{1})\zeta(P_{1}\to
P_{4})\differential\bm{r}_{1}\differential\bm{\hat{\Omega}}_{1}\differential
E_{1}\text{.}$ (24)

In the last step we have introduced the function $\zeta$:

$\zeta(P_{1}\to P_{4})=\int\differential
E_{3}\frac{\nu_{f}(\bm{r}_{4},E_{3})\Sigma_{f}(\bm{r}_{4},E_{3})f_{\text{fiss}}\left(\bm{\hat{\Omega}}_{4},E_{4}|\bm{r}_{4},\frac{\bm{r}_{4}-\bm{r}_{1}}{\absolutevalue{\bm{r}_{4}-\bm{r}_{1}}},E_{3}\right)\bar{f}\left(\frac{\bm{r}_{4}-\bm{r}_{1}}{\absolutevalue{\bm{r}_{4}-\bm{r}_{1}}},E_{3}|\bm{r}_{1},\bm{\hat{\Omega}}_{1},E_{1}\right)}{k\absolutevalue{\bm{r}_{4}-\bm{r}_{1}}^{2}}\\\
\exp\left(-\int\limits_{0}^{\absolutevalue{\bm{r}_{4}-\bm{r}_{1}}}\Sigma_{t}\left(\bm{r}_{1}+u\frac{\bm{r}_{4}-\bm{r}_{1}}{\absolutevalue{\bm{r}_{4}-\bm{r}_{1}}},E_{3}\right)\differential
u\right)$ (25)

Here $\zeta(P_{1}\to P_{4})$ is the transition kernel for a particle starting
at $P_{1}$, undergoing a collision, then a flight, and then producing fission
particles at $P_{4}$.

We now wish to construct an estimator for the expected fission emission
density at a point $Q\in\mathcal{R}$. Our estimator operates on events where a
fission particle is emitted at $P_{4}\in\mathcal{R}$, from a particle
originally entering a collision at $P_{1}$. To be unbiased, our estimator
$\vartheta(P_{1}\to P_{4}|\mathcal{R},Q)$ for the fission emission density
must have the property [4]

$\int\zeta(P_{1}\to P_{4})\vartheta(P_{1}\to P_{4}|\mathcal{R},Q)\differential
P_{4}=\zeta(P_{1}\to Q)\text{.}$ (26)

With regional cancellation, we would effectively like to define an estimator
$\vartheta$ for the fission emission density at $Q$ where a portion of the
fission emission density is located exactly at $Q$, and the remaining portion
is uniformly distributed within the phase space region $\mathcal{R}$. We shall
define this estimator to have the form

$\vartheta_{\eta}(P_{1}\to
P_{4}|\mathcal{R},Q)=(1-\eta)\delta(Q-P_{4})+\frac{\eta}{\mathcal{V}_{\mathcal{R}}}\text{.}$
(27)

Here, $\mathcal{V}_{\mathcal{R}}$ is the generalized phase space volume
occupied by $\mathcal{R}$, and $\eta$ is the portion of the fission emission
density that we wish to uniformly distribute within $\mathcal{R}$.444Note that
the parameter $\eta$ may take any value (real or complex); in particular, it
is not required to lie in the $[0,1]$ interval. If $\eta$ is taken to be a
constant with respect to $P_{4}$, then, upon evaluation of the left-hand side
of Eq. (26), using Eq. (27), we obtain:

$\int\zeta(P_{1}\to P_{4})\vartheta_{\eta}(P_{1}\to
P_{4}|\mathcal{R},Q)\differential P_{4}=\\\ \int\zeta(P_{1}\to
P_{4})\left[(1-\eta)\delta(Q-P_{4})+\frac{\eta}{\mathcal{V}_{\mathcal{R}}}\right]\differential
P_{4}=\\\ (1-\eta)\zeta(P_{1}\to
Q)+\frac{\eta}{\mathcal{V}_{\mathcal{R}}}\int\zeta(P_{1}\to
P_{4})\differential P_{4}\text{.}$ (28)

Comparing Eq. (28) and Eq. (26), it is clear that the only unbiased option is
$\eta=0$, corresponding to no cancellation. It is permissible however to allow
$\eta=\eta(P_{1},P_{4})$ to be both a function of $P_{1}$ and $P_{4}$, as
$\vartheta$ is already a function of these parameters. Using the ansatz

$\eta(P_{1},P_{4})=\frac{\beta}{\zeta(P_{1}\to P_{4})}\text{,}$ (29)

inspired by Eq. (3), we see that

$\int\zeta(P_{1}\to P_{4})\vartheta_{\eta(P_{1},P_{4})}(P_{1}\to
P_{4}|\mathcal{R},Q)\differential P_{4}=\\\ \int\zeta(P_{1}\to
P_{4})\left[\left(1-\frac{\beta}{\zeta(P_{1}\to
P_{4})}\right)\delta(Q-P_{4})+\frac{\beta}{\mathcal{V}_{\mathcal{R}}\zeta(P_{1}\to
P_{4})}\right]\differential P_{4}=\\\ \zeta(P_{1}\to
Q)-\beta+\beta=\zeta(P_{1}\to Q)\text{,}$ (30)

which, compared with right-hand side of Eq. (26), shows that this choice leads
to an unbiased estimator.

Equation (30) indicates that we are allowed to distribute a factor
$\beta/\zeta(P_{1}\to P_{4})$ of the fission particle uniformly within
$\mathcal{R}$, so long as $\beta$ has no functional dependence on $P_{4}$.
This requirement on $\beta$ is essential to ensure that, after integrating
over $P_{4}$, both $\beta$ terms will cancel; note however that $\beta$ is
allowed to depend on $P_{1}$. We are therefore allowed to pick $\beta=0$
whenever it is convenient, so long as information from $P_{4}$ is not used to
make this choice.

Additionally, Eq. (30) indicates that we must require $\zeta(P_{1}\to
P_{4})\neq 0\text{ }\forall P_{4}\in\mathcal{R}$. If this is not the case,
then $\eta$ is undefined. In particular, this implies that we require
$\Sigma_{f}(\bm{r}_{4},E_{3})>0$ everywhere within our cancellation region. We
must also require
$f_{\text{fiss}}(\bm{\hat{\Omega}}_{4},E_{4}|\bm{r}_{4},\frac{\bm{r}_{4}-\bm{r}_{1}}{\absolutevalue{\bm{r}_{4}-\bm{r}_{1}}},E_{3})>0\text{
}\forall P_{4}\in\mathcal{R}$; as fission is nearly perfectly isotropic, the
angular component is not problematic, but the energy component could indeed be
zero for very low energies, and care must therefore be taken when selecting
the energy bounds for $\mathcal{R}$. Despite these restrictions, we are given
some liberty as to the definition of $\mathcal{R}$, as it is allowed to be
non-convex.

The ansatz of Eq. (29) has the following remarkable property. Suppose that
$\zeta(P_{1}\to P_{4})$ has the structure

$\zeta(P_{1}\to P_{4})=h_{0}(P_{1})\zeta_{0}(P_{1}\to P_{4})\text{.}$ (31)

Consider the estimators

$\displaystyle\vartheta_{\zeta}$ $\displaystyle=\vartheta\left(P_{1}\to
P_{4}|\mathcal{R},Q,\eta=\frac{\beta}{\zeta}\right)$ (32)
$\displaystyle\vartheta_{\zeta_{0}}$ $\displaystyle=\vartheta\left(P_{1}\to
P_{4}|\mathcal{R},Q,\eta=\frac{\beta_{0}}{\zeta_{0}}\right)\text{.}$ (33)

The two estimators are actually identical for $\beta_{0}=\beta/h_{0}$. In
other words, any factor in $\zeta(P_{1}\to P_{4})$ that is independent of
$P_{4}$ can be pulled out of the definition of $\zeta$ and still yield an
unbiased estimator for the fission emission density.

Finally, note that the integral over $E_{3}$ in the definition of $\zeta$
could be somewhat problematic and/or expensive to compute on the fly in a
continuous-energy Monte Carlo code. It effectively corresponds to averaging
over all possible nuclides, reaction channels, and energies, which could have
been sampled in determining $P_{3}$, and leading to a fission particle at
$P_{4}$. We therefore would like to determine if there is an alternative,
simpler, unbiased option.

####  III.2.2 Expected-value estimators for cancellation

We have seen that, in order to be able to distribute a portion of the fission
emission density uniformly over $\mathcal{R}$, we need to replace the analogue
estimator of the fission emission density with a suitable form of expected-
value estimator [24]. The idea of expected-value estimators is to replace a
sampled quantity with the expected value of the same quantity over its
probability distribution. Next-event estimators are particular forms of
expected-value estimators that average the sampled quantity over the following
event in the stochastic process.

Let us now evaluate if a next-event estimator is suitable for the purpose of
cancellation. Consider a next-fission estimator for the fission emission
density; such an estimator is applied to particles undergoing a collision at
$P_{3}=(\bm{r}_{3},\bm{\hat{\Omega}}_{3},E_{3})$ and yields the expected
fission emission density at a generic point
$Q=(\bm{r},\bm{\hat{\Omega}},E)\in\mathcal{R}$. Since fission does not change
the position of particles, the contribution of the next-fission estimator
vanishes everywhere except for $\bm{r}=\bm{r}_{3}$. Therefore, a next-fission
estimator is not able to yield a non-vanishing contribution at all the points
in cancellation region $\mathcal{R}$.

Thus, in order for cancellation to be possible, we need to include more than
one event in our expected-value estimator, i.e. we need to use at least a
next-next-event estimator, or possibly an estimator of even higher order. It
is then crucial to determine the number of events that our estimator needs to
look ahead and average over, in order to yield a non-vanishing contribution to
all the phase space points in the cancellation region. Indeed, we want to
_minimize_ the number of look-ahead events, because the evaluation of
expected-value estimators becomes more and more cumbersome as the number of
look-ahead events increases.

Consider now the possibility of a next-flight-fission estimator. In our
notation, such an estimator acts on particles emitted at $P_{2}$ and yields
the expected fission emission density at $Q$, averaged over all possible
flights from $P_{2}$ and all possible fission events. Since the flight
operator does not modify the particle direction (see Eq. (11)), the expected
fission density vanishes everywhere except at positions reachable from
$\bm{r}_{2}$ with direction $\bm{\hat{\Omega}}_{2}$. In general, this does not
cover the whole cancellation region, except in the one-dimensional case [23,
22].

It is now probably clear that a next-collision-flight-fission estimator should
_in general_ yield a non-vanishing contribution to the fission emission
density at all phase space points within $\mathcal{R}$. In other words, given
a particle undergoing a collision at $P_{1}$, the expected fission emission
density (averaged over the next collision, flight, and fission) should not
vanish anywhere within $\mathcal{R}$. This corresponds to Eq. (25) above and
justifies the construction of the previous section.

Two remarks are in order here. First, there are cases where even a next
collision-flight-fission estimator is not sufficient to achieve a non-
vanishing expected fission emission density at all points within
$\mathcal{R}$. Indeed, the collision between $P_{1}$ and $P_{2}$ may be
subject to kinematic constraints, and Eq. (25) shows that the expected fission
emission density vanishes if
$\bar{f}\left(\frac{\bm{r}-\bm{r}_{1}}{\absolutevalue{\bm{r}-\bm{r}_{1}}},E_{3}|\bm{r}_{1},\bm{\hat{\Omega}}_{1},E_{1}\right)=0$
for some $\bm{r}$ in cancellation region $\mathcal{R}$. Second, an estimator
based on Eq. (25) would require the evaluation of the integral over $E_{3}$ at
every collision, which is impractical. In fact, regional cancellation attempts
to perform the cancellation algorithm a posteriori, after flights have already
been sampled, and fission particles have been produced. Since we have already
sampled a Monte Carlo history from $P_{1}$ through $P_{2}$, $P_{3}$, and
$P_{4}$, we would like to reuse as much information as possible from the
sampled history to remove part of the fission density from $P_{4}$ and
redistribute it uniformly within $\mathcal{R}$.

Thus, our expected-value estimator needs to average over sufficiently many
event samplings to be able to “see” the whole region $\mathcal{R}$; at the
same time, we want our estimator to average over the strict minimum number of
samplings. Each additional real variable that we average over introduces an
extra integration in the expression of the expected fission emission density
and reduces the usefulness of the $P_{1}\to P_{4}$ history that we have
already sampled.

####  III.2.3 Conditioning on Sampled Scatter

For the subsequent analysis, it is useful to consider a different form of the
collision kernel, more aligned with how most continuous-energy Monte Carlo
codes sample a collision event. While Eq. (17) presents the collision kernel
in terms of the averaged macroscopic cross sections and yields, most
continuous-energy Monte Carlo codes do not handle collisions in such a manner.
In production-level codes, microscopic cross sections are tabulated for
different nuclides and different reaction channels (elastic, level inelastic,
etc.) [17, 14, 15, 16]. Each combination of nuclide and reaction channel has
an independent transfer function for each type of non-capture collision. The
concentration $N_{i}(\bm{r}^{\prime})$ of nuclide $i$ is a function of
position, and the total microscopic cross section
$\sigma_{i}(\bm{r}^{\prime},E^{\prime})$ is a function of position and
energy.555The microscopic cross section is typically given as a function of
temperature and energy. However, since the temperature is a function of
position, we have chosen to present the microscopic cross section as a
function of position and energy, to avoid the introduction of a superfluous
variable. The total macroscopic cross section is

$\Sigma_{t}(\bm{r}^{\prime},E^{\prime})=\sum_{i}N_{i}(\bm{r}^{\prime})\sigma_{i}(\bm{r}^{\prime},E^{\prime})\text{.}$
(34)

At a collision site, we select the nuclide with which our particle will
undergo a collision: nuclide $i$ is chosen with probability
$N_{i}(\bm{r}^{\prime})\sigma_{i}(\bm{r}^{\prime},E^{\prime})/\Sigma_{t}(\bm{r}^{\prime},E^{\prime})$.
With nuclide $i$ having been sampled, a reaction channel $m$ must next be
sampled. If we let $\sigma_{i,m}(\bm{r}^{\prime},E^{\prime})$ be the partial
microscopic cross section for channel $m$, then the total microscopic cross
section is

$\sigma_{i}(\bm{r}^{\prime},E^{\prime})=\sum_{m}\sigma_{i,m}(\bm{r}^{\prime},E^{\prime})\text{,}$
(35)

and channel $m$ will be selected with probability
$\sigma_{i,m}(\bm{r}^{\prime},E^{\prime})/\sigma_{i}(\bm{r}^{\prime},E^{\prime})$.
This channel has an associated yield of $\nu_{i,m}(E^{\prime})$, and transfer
function
$f_{i,m}\left(\bm{\hat{\Omega}},E|\bm{\hat{\Omega}}^{\prime},E^{\prime}\right)$.
Continuous-energy nuclear data files typically give $f_{i,m}$ as a product of
a marginal PDF in energy and a conditional PDF in direction:

$f_{i,m}\left(\bm{\hat{\Omega}},E|\bm{\hat{\Omega}}^{\prime},E^{\prime}\right)=\\\
f_{i,m}\left(E|\bm{\hat{\Omega}}^{\prime},E^{\prime}\right)f_{i,m}\left(\bm{\hat{\Omega}}|\bm{\hat{\Omega}}^{\prime},E^{\prime},E\right)\text{.}$
(36)

When this is the case, the energy $E$ is first sampled from the marginal PDF,
and the direction is subsequently sampled from the conditional PDF. With these
provisions, it is then possible to write the collision kernels, Eqs. (13) and
(14), as

$\displaystyle C_{s}(P^{\prime}\to
P)=\frac{\delta\left(\bm{r}-\bm{r}^{\prime}\right)}{\Sigma_{t}(\bm{r}^{\prime},E^{\prime})}\sum_{i}N_{i}(\bm{r}^{\prime})\sum_{\begin{subarray}{c}m\\\
m\not=\text{fiss}\end{subarray}}\nu_{i,m}(E^{\prime})\sigma_{i,m}(\bm{r}^{\prime},E^{\prime})f_{i,m}\left(E|\bm{\hat{\Omega}}^{\prime},E^{\prime}\right)f_{i,m}\left(\bm{\hat{\Omega}}|\bm{\hat{\Omega}}^{\prime},E^{\prime},E\right)$
(37a) $\displaystyle C_{f}(P^{\prime}\to
P)=\frac{\delta\left(\bm{r}-\bm{r}^{\prime}\right)}{\Sigma_{t}(\bm{r}^{\prime},E^{\prime})}\sum_{i}N_{i}(\bm{r}^{\prime})\nu_{i,\text{fiss}}(E^{\prime})\sigma_{i,\text{fiss}}(\bm{r}^{\prime},E^{\prime})f_{i,\text{fiss}}\left(E|\bm{\hat{\Omega}}^{\prime},E^{\prime}\right)f_{i,\text{fiss}}\left(\bm{\hat{\Omega}}|\bm{\hat{\Omega}}^{\prime},E^{\prime},E\right)\text{.}$
(37b)

Based on the form of Eqs. (37), we introduce then the concept of a “partial
collision point”, indicating that the required pieces of information are
sampled incrementally when performing a collision. Examples of partial
collision points would be the state where we have sampled only the nuclide, or
the nuclide and the channel, or the nuclide, reaction channel, and energy.
With the concept of a partial collision, it is then evident that there is a
partial collision point between $P_{1}$ and $P_{2}$, where the particle has
selected an isotope to collide with ($i$), a reaction channel ($m$), and even
an outgoing energy ($E_{2}$), but has yet to select a direction
$\bm{\hat{\Omega}}_{2}$ out of the collision. In general, this partial
collision point between $P_{1}$ and $P_{2}$ is the strict minimum number of
steps we must look back in a particle’s history, in order to see a nonzero
fission emission density everywhere within $\mathcal{R}$ (assuming that is is
possible to scatter into all directions subtended by $\mathcal{R}$). This
state is accessible in a Monte Carlo simulation, as the nuclide, reaction
channel, and energy $E_{3}$ were all sampled when producing the fission
particle at $P_{4}$, and this information can be stored with the particle. The
transition kernel from $P_{1}$ to $P_{4}$, given a collision with nuclide $i$
in reaction channel $m$ and outgoing energy $E_{3}$, is then

$\zeta(P_{1}\to
P_{4}|i,m,E_{3})=\frac{\nu_{f}(\bm{r}_{4},E_{3})\Sigma_{f}(\bm{r}_{4},E_{3})f_{\text{fiss}}\left(\bm{\hat{\Omega}}_{4},E_{4}|\bm{r}_{4},\frac{\bm{r}_{4}-\bm{r}_{1}}{\absolutevalue{\bm{r}_{4}-\bm{r}_{1}}},E_{3}\right)f_{i,m}\left(\frac{\bm{r}_{4}-\bm{r}_{1}}{\absolutevalue{\bm{r}_{4}-\bm{r}_{1}}}|\bm{\hat{\Omega}}_{1},E_{1},E_{3}\right)}{k\absolutevalue{\bm{r}_{4}-\bm{r}_{1}}^{2}}\times\\\
\exp\left(-\displaystyle\int\limits_{0}^{\absolutevalue{\bm{r}_{4}-\bm{r}_{1}}}\Sigma_{t}\left(\bm{r}_{1}+u\frac{\bm{r}_{4}-\bm{r}_{1}}{\absolutevalue{\bm{r}_{4}-\bm{r}_{1}}},E_{3}\right)\differential
u\right)\text{.}$ (38)

This is now quite reminiscent of the fission density function which we used in
our previous work, as summarized in Sec. II [22]. It is worth stressing that
Eq. (38) uses the macroscopic fission cross section $\Sigma_{f}$ and the
average fission transfer function $f_{\text{fiss}}$, which are averaged over
all fissile nuclides at $\bm{r}_{4}$. In general, these quantities might vary
within $\mathcal{R}$, due to spatial dependence in the nuclide concentrations
and temperature. Examining under what circumstances the estimator
$\vartheta(P_{1}\to P_{4}|\mathcal{R},Q,\eta)$ is unbiased for the transition
kernel presented in Eq. (38), it is straightforward to observe that this
condition is met for the choice of $\eta=\beta/\zeta\left(P_{1}\to
P_{4}|i,m,E_{3}\right)$.

Equation (38) (and its associated estimator) are subject to the same
constraints as Eq. (25), as discussed in Sec. III.2.2: namely, the expected
fission emission density must be non-zero at all the points in the
cancellation region. In particular, it is required that
$f_{i,m}\left(\frac{\bm{r}_{4}-\bm{r}_{1}}{\absolutevalue{\bm{r}_{4}-\bm{r}_{1}}}|\bm{\hat{\Omega}}_{1},E_{1},E_{3}\right)>0\text{
}\forall\text{ }\bm{r}_{4}\in\mathcal{R}$. A special case arises when the
reaction channel $m$ uses a delta distribution for either the energy or
direction (such as in level inelastic scattering). If such a channel was
selected during the last collision, then the value of $\zeta(P_{1}\to P_{4})$
is infinite at the collision point (as we are evaluating the delta
distribution at the singularity), and it vanishes almost everywhere within
$\mathcal{R}$. Thus, such channels do not generally partake in cancellation,
as the uniform portion would then necessarily be zero according to Eq. (29).
Finally, not all nuclear data facilitates the decomposition provided by Eq.
(36). Sometimes the joint PDF might be provided as a marginal PDF in
direction, and a conditional PDF in energy. If this is the case, one must go
back to the partial collision before having sampled the direction, in order to
see the entire cancellation region.

If the region $\mathcal{R}$ contains only one material, which is completely
homogeneous in nuclide concentrations and temperature, then $\nu_{f}$ and
$\Sigma_{f}$ are independent of $P_{4}$. As we have discussed in Sec. III.2.1,
these factors may be removed from the the definition of $\zeta$, without
compromising the unbiasedness of the method. In addition, if fission is
assumed to be perfectly isotropic (a frequent assumption), and if the fission
energy $E_{4}$ is completely independent of the incident energy and direction,
then we do not actually need to perform cancellation on the fission emission
density, but only on the collision density at $\bm{r}_{4}$. For the multi-
group benchmark in Ref. 22 in which we previously demonstrated exact regional
cancellation in 3D, we used homogeneous cancellation regions, where both
fission and scattering were isotropic, and the fission energy was also assumed
to be independent of the incident energy $E_{3}$. This indicates that $\zeta$
could be simplified to

$\zeta(P_{1}\to P_{4}|E_{3})=\\\
\frac{\exp\left(-\displaystyle\int\limits_{0}^{\absolutevalue{\bm{r}_{4}-\bm{r}_{1}}}\Sigma_{t}\left(\bm{r}_{1}+u\frac{\bm{r}_{4}-\bm{r}_{1}}{\absolutevalue{\bm{r}_{4}-\bm{r}_{1}}},E_{3}\right)\differential
u\right)}{\absolutevalue{\bm{r}_{4}-\bm{r}_{1}}^{2}}$ (39)

and still result in an unbiased cancellation method.

####  III.2.4 Non-Uniform Cancellation

Suppose now that, instead of distributing some fission emission density
uniformly over $\mathcal{R}$, we would like to distribute fission emission
density according to an arbitrary function, $D(Q)$. In this case, our
estimator must be modified as

$\vartheta_{D}(P_{1}\to P_{4}|\mathcal{R},Q,\eta)=\\\ \frac{\zeta(P_{1}\to
Q|i,m,E_{3})-D(Q)\beta}{\zeta(P_{1}\to Q|i,m,E_{3})}\delta(Q-P_{4})+\\\
\frac{\beta}{\zeta(P_{1}\to
P_{4}|i,m,E_{3})}\frac{D(Q)}{\mathcal{V}_{\mathcal{R}}}\text{.}$ (40)

This can be shown by following the same approach taken in Eq. (30).

Here, we have placed a portion $(\zeta(P_{1}\to
Q|i,m,E_{3})-D(Q)\beta)/\zeta(P_{1}\to Q|i,m,E_{3})$ of the weight at the
sampled point $P_{4}$, and a portion $\beta/\zeta(P_{1}\to Q|i,m,E_{3})$ of
the particle is distributed according to $D(Q)$. It is only assumed here that
$D(Q)$ is dimensionless and that it is piece-wise continuous. In theory, there
is no reason that $D(Q)$ could not be negative, or even complex valued666While
this paper only considers particles with a single real-valued statistical
weight, some transport problems require that particles carry a complex weight
[9, 10].; if $D(Q)$ is negative or complex, it might not be possible to sample
it directly, but such a situation might be treated using e.g. importance
sampling [4]. Finally, we note that for the choice of $D(Q)=1$ the case of
uniform cancellation is retrieved.

####  III.2.5 Fission Emission Density Function with Delta-Tracking Schemes

In the above derivations, we have often made use of the non-homogeneous
exponential distribution

$\Sigma_{t}(\bm{r}+d\bm{\hat{\Omega}},E)\exp\left(-\int\limits_{0}^{d}\Sigma_{t}(\bm{r}+s\bm{\hat{\Omega}},E)\differential
s\right)\text{,}$ (41)

which occurs in the flight kernel $T(P^{\prime}\to P)$ and in the transition
kernel $\zeta(P^{\prime}\to P)$. This distribution is sampled when trying to
determine the distance $d$ a particle will travel from initial position
$\bm{r}$ along direction $\bm{\hat{\Omega}}$, before undergoing a collision.
For the case of piece-wise constant macroscopic cross sections, this
distribution is straightforward to sample for $d$. However, when the
macroscopic cross section is not piece-wise constant, more sophisticated
methods than direct sampling are often employed [11]. Delta-tracking and
negative-weighted delta-tracking are two such methods that sample the distance
to collision using a sampling cross section $\Sigma_{\text{smp}}(E)$, and then
sample whether a real or virtual collision has occurred with a specific
criterion [19, 21, 20]; in the distinct case of negative-weighted delta
tracking, a weight modifier may be additionally applied to the particle’s
weight [21, 20]. In a real collision, the particle undergoes a reaction as
normal. In a virtual collision, the particle’s energy and direction do not
change; this event is usually known as _delta scattering_. The particle
continues to sample new flight distances and to move to the new location,
until a real collision is sampled. Coleman [25] and Legrady et al. [21] have
previously provided evidence as to why such sampling methods are unbiased.

It is possible to include virtual collisions in the transport equations by
modifying the flight kernel to be

$T_{DT}(P^{\prime}\to P)=\\\
\frac{\Sigma_{\text{smp}}(E^{\prime})\exp\left(-\Sigma_{\text{smp}}(E^{\prime})\absolutevalue{\bm{r}-\bm{r}^{\prime}}\right)}{\absolutevalue{\bm{r}-\bm{r}^{\prime}}^{2}}\\\
\delta\left(\bm{\hat{\Omega}}-\frac{\bm{r}-\bm{r}^{\prime}}{\absolutevalue{\bm{r}-\bm{r}^{\prime}}}\right)\delta\left(\bm{\hat{\Omega}}^{\prime}-\bm{\hat{\Omega}}\right)\delta\left(E-E^{\prime}\right)\text{,}$
(42)

and the collision kernel to be

$C_{DT}(P^{\prime}\to P)=\\\
\frac{\Sigma_{t}(\bm{r}^{\prime},E^{\prime})}{\Sigma_{\text{smp}}(E^{\prime})}C(P^{\prime}\to
P)+\\\
\left(1-\frac{\Sigma_{t}(\bm{r}^{\prime},E^{\prime})}{\Sigma_{\text{smp}}(E^{\prime})}\right)\delta(E-E^{\prime})\delta(\bm{\hat{\Omega}}-\bm{\hat{\Omega}}^{\prime})\delta(\bm{r}-\bm{r}^{\prime})\text{.}$
(43)

These equations are valid for both delta tracking and negative-weighted delta
tracking. From Eq. (42), it is noticed that the PDF for leaving a collision
site at $\bm{r}^{\prime}$ and flying to $\bm{r}$ and inducing a fission (given
that we are flying in the direction of $\bm{r}$, i.e.
$\bm{\hat{\Omega}}^{\prime}=\frac{\bm{r}-\bm{r}^{\prime}}{\absolutevalue{\bm{r}-\bm{r}^{\prime}}}$)
is

$\Sigma_{f}(\bm{r},E^{\prime})\exp\left(-\Sigma_{\text{smp}}(E^{\prime})\absolutevalue{\bm{r}-\bm{r}^{\prime}}\right)\text{,}$
(44)

which is exactly the form presented in Eq. (1). Thus, delta-tracking-like
algorithms provide the large advantage of not requiring the integration of the
total cross section along the flight path. This makes them interesting for the
purpose of performing exact regional cancellation.

The form of Eq. (42) is valid regardless of whether the collision at
$\bm{r}^{\prime}$ was real or virtual. Equation (43) shows that the angular
distribution for virtual collisions is singular, because it is described by a
delta distribution. As discussed in Sec. III.2.3, channels with singular
distributions are not allowed to partake in cancellation, i.e. we need to set
$\beta=0$ for all such channels. For the particular case of virtual
collisions, however, another treatment is possible. At the site where the
virtual collision took place, there was a probability that the particle could
have instead undergone a real collision. We can therefore imagine “splitting”
the particle before the collision. A weight
$w(1-\Sigma_{t}(\bm{r},E^{\prime})/\Sigma_{\text{smp}}(E^{\prime}))$ is
considered to undergo a virtual collision, and have its next collision at
$P_{3}$; this virtual collision portion cannot be used in cancellation, as the
angular distribution was a delta distribution, and the uniform component is
then always zero, as explained in Sec. III.2.3. The rest of the particle
weight, namely
$w\Sigma_{t}(\bm{r},E^{\prime})/\Sigma_{\text{smp}}(E^{\prime})$, is
considered to undergo a real collision and have its next collision at $P_{3}$,
like the virtual part. However, this part can also partake in cancellation.
The point-wise fission particle weight which must remain at the sampled
fission particle site is then

$\displaystyle
w\left(1-\frac{\Sigma_{t}(\bm{r},E^{\prime})}{\Sigma_{\text{smp}}(E^{\prime})}\right)+w\frac{\Sigma_{t}(\bm{r},E^{\prime})}{\Sigma_{\text{smp}}(E^{\prime})}\left(1-\frac{\beta}{\zeta(P^{\prime}\to
P)}\right)=$ $\displaystyle w\left(1-\frac{\beta}{\zeta(P^{\prime}\to
P)}\frac{\Sigma_{t}(\bm{r},E^{\prime})}{\Sigma_{\text{smp}}(E^{\prime})}\right)=$
$\displaystyle w\left(1-\frac{\beta^{\prime}}{\zeta(P^{\prime}\to
P)}\right)\text{,}$ (45)

where we have set $\beta^{\prime}=\beta\Sigma_{t}/\Sigma_{\text{smp}}$. Thus,
splitting shows that _all_ collisions can be assumed to partake in
cancellation as if they were real, because the presence of virtual collisions
only affects the choice of $\beta$. Since the estimator is unbiased for any
$\beta$, the factor
$\Sigma_{t}(\bm{r},E^{\prime})/\Sigma_{\text{smp}}(E^{\prime})$ is not
necessary. However, note that this approach is only unbiased so long as at the
virtual collision site the real component of the scattering kernel for forward
scattering with no energy change is not zero (i.e. $C(P\to P)\not=0$).
Otherwise, the real collision component cannot be left at $P_{3}$, as it is
impossible to have had a real collision with forward scattering and no change
in energy. This was possible in our previous multi-group example, because in-
group scattering was always allowed and all scattering was assumed to be
isotropic; however, this might not be as trivial in a continuous-energy
setting.

##  IV Optimization of Cancellation Efficiency

We now turn our attention to the optimal choice of the free parameter $\beta$
of the cancellation estimator, used to calculate the weight that can be
uniformly distributed over the cancellation region. In Booth and Gubernatis’s
seminal paper [23] and in our previous work [22], $\beta$ was chosen to be the
minimum of the expected fission density over the cancellation region. This
choice has the advantage of being relatively easy to evaluate, but it is not
necessarily the most efficient one. In this section we attempt to introduce a
better strategy to determine the cancellation parameters.

In order to optimize for the cancellation efficiency, one must first properly
define the quantity to be optimized. As we mentioned in Sec. II, the maximum
amount of cancellation will occur when the sum of the absolute value of all
the weights remaining after cancellation in the region has been minimized. For
$N$ particles which initially land in a cancellation region, we define the
absolute value of all weight in a region after cancellation as

$\Gamma_{1}=\sum\limits_{k=1}^{N}\absolutevalue{w_{k,p}}+\absolutevalue{\sum\limits_{k=1}^{N}w_{k,u}},$
(46)

with $w_{k,p}$ being the point-wise weight of particle $k$, and $w_{k,u}$ the
uniform weight portion of particle $k$.777While there are $N$ particles in the
cancellation region before the cancellation operations have been carried out,
there will be more than $N$ particles after cancellation, due to the new
uniform particles which are created during the cancellation process. Equation
(46) is the total post-cancellation weight discussed in Sec. II, Eq. (7). As
each particle has a different value for the cancellation parameter $\beta$, we
then may substitute to obtain

$\Gamma_{1}=\sum\limits_{k=1}^{N}\absolutevalue{\frac{\zeta_{k}-\beta_{k}}{\zeta_{k}}w_{k}}+\absolutevalue{\sum\limits_{k=1}^{N}\frac{\beta_{k}}{\zeta_{k}}w_{k}}\text{,}$
(47)

where $w_{k}$ and $\beta_{k}$ are respectively the pre-cancellation weight and
the cancellation parameter of the $k$-th particle, and
$\zeta_{k}=\zeta(P_{k}\rightarrow P^{\prime}_{k})$ is the expected fission
density of the $k$-th particle, which is assumed to have had its previous
collision in $P_{k}$ and its fission event in $P^{\prime}_{k}$ (note that we
have simplified the notation here compared to Sec. III; for a given particle,
$P_{k}$ and $P^{\prime}_{k}$ respectively correspond to $P_{1}$ and $P_{4}$).

Due to the presence of the absolute values, it is quite difficult to optimize
the expression of $\Gamma_{1}$ analytically with respect to $\beta_{k}$. We
instead define a modified quantity $\Gamma_{2}$, which shares a minimum with
$\Gamma_{1}$:

$\Gamma_{2}=\sum\limits_{k=1}^{N}\bigg{(}\frac{\zeta_{k}-\beta_{k}}{\zeta_{k}}w_{k}\bigg{)}^{2}+\bigg{(}\sum\limits_{k=1}^{N}\frac{\beta_{k}}{\zeta_{k}}w_{k}\bigg{)}^{2}.$
(48)

We now wish to obtain the set of optimal values $\beta_{k}$ that minimize
$\Gamma_{2}$. To remain unbiased, we are not allowed to calculate $\beta_{k}$
based on the phase space coordinates $P^{\prime}_{k}$ where the particle $k$
landed in the cancellation region (this was made evident in Eq. (30)). As
$\zeta_{k}=\zeta(P_{k}\rightarrow P^{\prime}_{k})$, we cannot directly
minimize Eq. (48). In the two subsequent sections, we will present two
reasonable options to avoid this problem, both possibly giving way to a method
of optimizing the regional cancellation algorithm. We go through the
optimization for each case, obtaining two different formulations for
calculating the set of optimal values for $\beta_{k}$.

###  IV.1 Replacing $\zeta_{k}$ with $\expectationvalue{\zeta_{k}}$

The first approach consists in averaging $\zeta_{k}$ over the entire phase
space of the region $\mathcal{R}$, such that

$\expectationvalue{\zeta_{k}}=\frac{\displaystyle\int_{\mathcal{R}}\zeta(P_{k}\rightarrow
P^{\prime}_{k})\differential
P^{\prime}_{k}}{\displaystyle\int_{\mathcal{R}}\differential P^{\prime}_{k}}.$
(49)

We may then replace $\zeta_{k}$ with $\expectationvalue{\zeta_{k}}$ in Eq.
(48), and optimize the new approximate form

$\Gamma^{*}_{2}=\sum_{k=1}^{N}\bigg{(}\frac{\expectationvalue{\zeta_{k}}-\beta_{k}}{\expectationvalue{\zeta_{k}}}w_{k}\bigg{)}^{2}+\bigg{(}\sum_{k=1}^{N}\frac{\beta_{k}}{\expectationvalue{\zeta_{k}}}w_{k}\bigg{)}^{2}\text{,}$
(50)

which is now independent of $P^{\prime}_{k}$. The detailed derivation for this
approach is presented in Appendix A, and the resulting equation for the
cancellation parameter $\beta_{k}$ is found to be

$\beta_{k}=\expectationvalue{\zeta_{k}}\left(1-\frac{S^{*}}{w_{k}}\right),$
(51)

where we make use of the definition

$S^{*}=\frac{W}{N+1},$ (52)

and

$W=\sum_{k=1}^{N}w_{k},$ (53)

$W$ being the net weight in the region $\mathcal{R}$ before cancellation. It
should also be mentioned that Eq. (51) would also be obtained if we first
minimized Eq. (48) with respect to $\beta_{k}$ and then averaged over
$P^{\prime}_{k}$.

###  IV.2 Optimization of $\expectationvalue{\Gamma_{2}}$

The second approach consists in averaging $\Gamma_{2}$ over the phase space of
the region $\mathcal{R}$, obtaining

$\expectationvalue{\Gamma_{2}}=\frac{\displaystyle\int_{\mathcal{R}}\Gamma_{2}\prod\limits_{k=1}^{N}\zeta_{k}\differential
P^{\prime}_{k}}{\displaystyle\int_{\mathcal{R}}\prod\limits_{k=1}^{N}\zeta_{k}\differential
P^{\prime}_{k}}.$ (54)

We may then optimize $\expectationvalue{\Gamma_{2}}$ instead of $\Gamma_{2}$.
The complete derivation is provided in Appendix B, and leads to a different
equation for $\beta_{k}$:

$\beta_{k}=\expectationvalue{\zeta_{k}}c_{k}\bigg{(}1-\frac{S}{w_{k}}\bigg{)}.$
(55)

Here, we have made use of the two following definitions:

$c_{k}=\bigg{(}2\expectationvalue{\zeta_{k}}\expectationvalue{\frac{1}{\zeta_{k}}}-1\bigg{)}^{-1},$
(56)

where the angle brackets have the same meaning as in Eq. (49), and

$S=\frac{\displaystyle\sum_{k=1}^{N}c_{k}w_{k}}{1+\displaystyle\sum_{k=1}^{N}c_{k}}.$
(57)

###  IV.3 Small Region Limit

For the two possible methods that we have outlined to minimize the weight
after cancellation, we are left with two different possibilities for the value
of $\beta_{k}$:

  * •

Equation (51) stems from replacing $\zeta_{k}$ with
$\expectationvalue{\zeta_{k}}$ in Eq. (48);

  * •

Equation (55) stems from averaging Eq. (48) over the phase space points
$P_{k}^{\prime}$, $\forall$ $k\in\\{1,\dots,N\\}$.

At first glance, these two choices of $\beta_{k}$ look quite different. Upon
closer inspection of the definition of $c_{k}$ in Eq. (56), we notice that
$c_{k}=1$ only if
$1/\expectationvalue{\zeta_{k}}=\expectationvalue{1/\zeta_{k}}$. When this is
the case, it then follows from Eq. (57) that $S=S^{*}$. This then indicates
that the two definitions of $\beta_{k}$ given by Eq. (51) and Eq. (55) are
equivalent only when
$1/\expectationvalue{\zeta_{k}}=\expectationvalue{1/\zeta_{k}}$. In general,
however, Jensen’s inequality implies
$1/\expectationvalue{\zeta_{k}}\leq\expectationvalue{1/\zeta_{k}}$ [26]. If a
cancellation region were defined such that it is small enough that one could
reasonably assume that $\zeta(P_{k}\rightarrow P^{\prime}_{k})$ is nearly
constant within the region, then $c_{k}\approx 1$, leading to the two methods
being equivalent. However, this would likely require such a small region that
it is very unlikely that any other particles would be located within the
region, which will make cancellation very ineffective.

##  V Monte Carlo Implementation

###  V.1 Cancellation with Distributed Memory Simulations

Most production Monte Carlo codes make use of distributed-memory parallel
computing techniques (such as MPI), although the exact algorithm used varies
from code to code [14, 17, 15, 16]. Generally speaking, distributed-memory
parallelization can pose a problem for cancellation, which is by construction
more efficient when there are more particles in each cancellation region. With
distributed-memory parallelization, the fission particles within a given
cancellation region will be distributed amongst several nodes. To ensure the
highest possible efficiency, cancellation must be performed on the entire
fission source. One method to do this is to send all of the fission particles
to the master node between power iteration generations, and then perform
cancellation only on the master node. Another option would be to use a method
inspired by domain decomposition [27], where certain nodes perform
cancellation for certain regions, and fission particles would need to be sent
to the node which corresponds to their cancellation region. This method is
certainly possible, but likely much more difficult to implement in production
Monte Carlo codes. For the proof of concept presented in this paper, we have
chosen to use the former method, sending all fission particles to the master
node for cancellation.

###  V.2 Calculation of $\expectationvalue{\zeta_{k}}$ and
$\expectationvalue{1/\zeta_{k}}$

In order to evaluate $\beta_{k}$ according to Eq. (51), we must have knowledge
of $\expectationvalue{\zeta_{k}}$. For Eq. (55), we additionally need
knowledge of $\expectationvalue{1/\zeta_{k}}$. In general, it is not possible
to analytically calculate either of these quantities for particle $k$, born
with phase space coordinates $P^{\prime}_{k}$ located within cancellation
region $\mathcal{R}$. However, it is possible to estimate both of these
quantities with a Monte Carlo sampling approach. For each particle $k$, we
know the phase space coordinates $P_{k}$ of its previous collision, and we
know the bounds of the outgoing phase space coordinates $P^{\prime}_{k}$ which
define the cancellation region $\mathcal{R}$.

Assume that a set of non-overlapping, hypercuboid cancellation regions are
imposed on top of the problem domain. Then, between each generation of power
iteration, the fission particles (having stored their parent’s previous phase
space coordinates $P_{k}$) may be sorted into the cancellation regions, based
on their phase space coordinates $P^{\prime}_{k}$. Once this is accomplished
for a given cancellation region $\mathcal{R}$, we may iterate over all
particles in $\mathcal{R}$, and estimate their values of
$\expectationvalue{\zeta_{k}}$ and additionally
$\expectationvalue{1/\zeta_{k}}$, depending on which optimisation algorithm is
chosen. The estimates for these quantities may be obtained using

$\expectationvalue{\zeta_{k}}\approx\frac{1}{N_{s}}\sum_{i=1}^{N_{s}}\zeta(P_{k}\rightarrow\tilde{P}_{i})$
(58)

and

$\expectationvalue{\frac{1}{\zeta_{k}}}\approx\frac{1}{N_{s}}\sum_{i=1}^{N_{s}}\frac{1}{\zeta(P_{k}\rightarrow\tilde{P}_{i})},$
(59)

respectively, where $N_{s}$ is the number of samples to be used in the
estimation, and the outgoing phase space coordinates $\tilde{P}_{i}$ are
pseudo-randomly sampled so that $\tilde{P}_{i}\in\mathcal{R}$ and
$\tilde{P}_{i}\sim\mathcal{U}(\mathcal{R})$. This is straightforward to
accomplish with cuboid regions.

With this approach, a better estimate of $\expectationvalue{\zeta_{k}}$ and
$\expectationvalue{1/\zeta_{k}}$ may be obtained by augmenting the number of
samples. As $N_{s}$ is increased, the error on the estimate of the two
expectation values will decrease according to $\mathcal{O}(1/\sqrt{N_{s}})$
[28]. This indicates that a large $N_{s}$ may be required to obtain a suitable
estimate of $\expectationvalue{\zeta_{k}}$ and
$\expectationvalue{1/\zeta_{k}}$. Even more problematic is the fact that
evaluating $\zeta(P_{k}\rightarrow\tilde{P}_{i})$ could be quite costly; this
is especially true in the case of continuous-energy neutron transport
problems, where many evaluations of scattering distributions would be
necessary. In order to reduce $N_{s}$, while still obtaining adequate
estimates for $\expectationvalue{\zeta_{k}}$ and
$\expectationvalue{1/\zeta_{k}}$, we propose the use of a quasi-random
technique, using a Sobol’ sequence [28] to sample the outgoing phase space
coordinates $\tilde{P}_{i}$. This approach generally has a better convergence
rate than using a pseudo-random number generator to sample $\tilde{P}_{i}$, as
it leads to a more uniform exploration of the phase space [28]. At any rate,
we stress that statistical uncertainties on the estimation of
$\expectationvalue{\zeta_{k}}$ and $\expectationvalue{1/\zeta_{k}}$ only
affect the efficiency of the cancellation method, as we have proved that the
method is unbiased for any values of the free parameters $\beta_{k}$.

###  V.3 Heterogeneous Cancellation Regions

In this work, we have proposed two possible approaches to selecting an optimal
value of $\beta$, to optimize the amount of weight which is cancelled. Neither
of these approaches requires the minimum value of the fission emission density
within the cancellation region $\mathcal{R}$. Hence, it is no longer necessary
to restrict the cancellation regions to be cuboids, as we were required to do
in Ref. 22. In light of our proposed sampling methods to estimate
$\expectationvalue{\zeta_{k}}$ and $\expectationvalue{1/\zeta_{k}}$ to obtain
the optimized cancellation parameter $\beta_{k}$, it is evident that a
rejection technique may be applied to isolate different material regions
within a given cuboid cancellation region. If we have a cancellation region
with fuel and water, then all of the fission particles are of course only born
inside the fuel portion, and the fission density everywhere in the water
should be zero. When sampling the random phase space coordinates
$\tilde{P}_{i}$, we must now add the requirement that $\tilde{P}_{i}$ have
spatial coordinates that are located inside of the fuel material.

Other special cases can be handled using this approach. For example,
cancellation can also be performed when there are two non-connected fuel
regions within the same cancellation mesh region. We may also have
cancellation regions which contain multiple different fuel regions. Using
rejection sampling to determine cancellation regions makes it easy to apply
regional cancellation to complex geometries encountered in realistic reactor
physics problems. Of course, the rejection procedure must also be applied to
the sampling of the phase space coordinates of the uniformly distributed
particles.

###  V.4 Monte Carlo Implementation in the open-source code MGMC

For our previous work on regional cancellation, a multi-group Monte Carlo code
called MGMC was used to test cancellation on a well-known reactor physics
benchmark. MGMC is a mini-app which has been developed to facilitate the fast
and easy implementation and testing of new transport algorithms. Being only
$\approx 13\text{ k}$ lines of code, it is much faster to test new ideas in
MGMC than it would be in a large industrial code. General 3D geometries are
supported using a standard constructive solid geometry formalism based on
surfaces, universes, and lattices, familiar to any user of other well-known
Monte Carlo codes [14, 17, 16, 15]. Different mesh tallies are available for
flux or reaction rates, with track-length or collision estimators. MGMC can
solve fixed-source, $k$-eigenvalue, and neutron noise problems, using both
shared and distributed memory parallelism. Shared memory parallelism is
implemented with OpenMP, while the distributed memory parallelism is
implemented using MPI. Different transport methods such as surface-tracking,
delta-tracking, and negative-weighted delta-tracking are also available. All
of the outlined cancellation algorithms have been implemented in MGMC, which
was used to run the simulations presented in the next section. MGMC has been
make publicly available as free software under the CeCILL v2.1 license [29].

##  VI Simulation Results

For our numerical simulations, in this Section we will make use of the
modified C5G7 benchmark which we introduced in our previous work [22]. The
C5G7 is a multi-group neutron transport benchmark which comes from the nuclear
reactor physics community, for the purpose of validating different codes [30].
Our modified version makes use of square profile fuel pins with side lengths
of $0.756\text{\,}\mathrm{cm}$, in lieu of cylindrical pins of radius
$0.54\text{\,}\mathrm{cm}$ as proposed in the original specifications. This
modification allows a regular $170\times 170\times 765$ mesh to be imposed on
top of the geometry over the fuel assemblies to act as cancellation regions,
and guaranteed that each cancellation region contained a unique material. For
continuity, we make use of the same cancellation mesh. For transport, we again
use the negative-weighted delta-tracking variant proposed by Carter et al.
[20], and identical sampling cross sections to the previous study: the
sampling cross section for the first group is 90% of the majorant cross
section, while all other sampling cross sections were taken to be the
majorant. This means that the sampling cross section underestimates the total
cross section in the first energy group for all fuel pins in the problem.
Whenever a virtual collision occurs for a particle in the first energy group,
inside a fuel pin, its weight will then change sign. Once a particle leaves
the first group, it is impossible for the sign to change at a collision
(although signs can possibly change during cancellation). Virtual collisions
lead to the presence of negative weights in the system, and we have shown that
weight cancellation is mandatory for $k$-eigenvalue power iteration problem to
converge when using negative-weighted delta-tracking [22]. All simulations
were initiated with $10^{6}$ particles, and ran for 2500 generations, with the
first 500 being discarded to allow for source convergence.

As we have shown in Ref. 22, the total weight of all the fission particles
between two generations increases without bound if weight cancellation is not
applied. The effect of cancellation is to limit the growth of $W_{\text{tot}}$
to a saturation value; the more efficient cancellation is, the lower the
saturation value will be. Thus, we have chosen to assess the efficiency of
cancellation by using the saturation value of $W_{\text{tot}}$, which is
calculated immediately after applying the cancellation procedure. Note that
$W_{\text{tot}}$ has a lower theoretical limit of $W_{\text{net}}$, which is
kept constant by normalizing all particle weights between generations [22].

###  VI.1 Comparison of Optimization Strategies

To determine which method of choosing $\beta_{k}$ leads to the most efficient
cancellation of positive and negative weights, the optimization techniques
described in Sec. IV.1 and Sec. IV.2 were compared against the original
implementation using the minimum value of the fission density within the
region. The values of $W_{\text{tot}}$ are plotted against the number of
generations in Figure 2. For comparison, curves corresponding to no
cancellation and approximate cancellation have also been presented.
Approximate cancellation imposes a mesh on top of the geometry, and sorts
fission particle into this mesh. The average weight of all particles in each
mesh element can then be calculated and assigned to the particles [31, 32].
This method is quite efficient, but is not exact, and imposes a bias on the
fission source and on the eigenvalue (though the bias can be made arbitrarily
small by using a sufficiently fine mesh).

Figure 2: Behavior of $W_{\text{tot}}$ as a function of generation, for
different cancellation methods, including no cancellation and approximate
cancellation.

If no cancellation technique is used, the total weight increases
exponentially, without bound. This phenomenon is expected when using negative-
weighted delta-tracking with $k$-eigenvalue power iteration, as described
previously [22]. When taking $\beta_{k}=\min_{\mathcal{R}}(\zeta_{k})$, an
asymptotic value of $W_{\text{tot}}\approx 4.9\cdot 10^{6}$ was seen. The most
efficient method of determining $\beta_{k}$ was that obtained from optimizing
$\expectationvalue{\Gamma_{2}}$ in Eq. (55), resulting in
$W_{\text{tot}}\approx 2.7\cdot 10^{6}$, almost half the amount of total
weight obtained with the minimum strategy. Calculating $\beta_{k}$ from Eq.
(51) for the case of replacing $\zeta_{k}$ with $\expectationvalue{\zeta_{k}}$
is less efficient than using the minimum value of $\zeta_{k}$ within the
cancellation region, resulting in $W_{\text{tot}}\approx 8.0\cdot 10^{6}$. It
is not known why this approximation does not perform as well as using the
minimum of $\zeta_{k}$, and this intriguing question calls for future
investigations.

Approximate cancellation yielded the lowest total weight (and therefore the
highest cancellation efficiency), with $W_{\text{tot}}\approx 1.4\cdot
10^{6}$, but is not an exact approach. Currently, we do not know of any way to
estimate, or to put a limit on the bias imposed by this method without running
several realizations, each with a different mesh size.

###  VI.2 Strategies for evaluating the average fission emission densities

The analysis in Sec. VI.1 shows that the optimal choice for determining
$\beta_{k}$ for the C5G7 benchmark is Eq. (55), from the optimization of
$\expectationvalue{\Gamma_{2}}$. We now consider the optimal strategy for
estimating the requisite values of $\expectationvalue{\zeta_{k}}$ and
$\expectationvalue{1/\zeta_{k}}$ for each particle. Figure 3 depicts the
behavior of $W_{\text{tot}}$ where points $\tilde{P}_{i}$ (from Eqs. (58) and
(59)) are sampled with either a pseudo-random number generator (PRNG) or a
Sobol’ sequence.

Figure 3: Behaviour of $W_{\text{tot}}$ as a function of generation, for
different cancellation parameters. For all curves, $\beta_{k}$ is determined
from the optimization of $\expectationvalue{\Gamma_{2}}$, Eq. (55). The
estimated values of $\expectationvalue{\zeta_{k}}$ and
$\expectationvalue{1/\zeta_{k}}$ are determined with a varying number of
points, using either a pseudo-random number generator (PRNG) or a Sobol’
sequence.

First of all, the spread among the different strategies for estimating
$\expectationvalue{\zeta_{k}}$ and $\expectationvalue{1/\zeta_{k}}$ is much
smaller than the spread among the different minimization strategies of Fig. 2.
It is observed that in general, when $N_{s}<50$, using Sobol’ points leads to
more efficient weight cancellation. This effect is most apparent for
$N_{s}=3$, where the Sobol’ sequence leads to approximately 9.2% less total
weight being transported, compared to the PRNG estimation strategy. The
increased efficiency observed in the Sobol’ points diminishes however with
increasing $N_{s}$. Sobol’ estimation gives a 2.8% improvement for $N_{s}=5$,
1.8% improvement for $N_{s}=10$, and only a 0.5% improvement for $N_{s}=50$.
This would indicate that the estimated values for
$\expectationvalue{\zeta_{k}}$ and $\expectationvalue{1/\zeta_{k}}$ start to
become independent of the evaluation strategy at around $N_{s}=50$.

In addition to achieving more weight cancellation, the Sobol’ points also have
the added benefit of being slightly easier to compute, as the quasi-random
numbers used the compute the points can be tabulated in advance, and written
in the code. All that is then needed is a table lookup to get a Sobol’ value,
whereas several mathematical operations must be performed to calculate each
value generated from a PRNG. However, a drawback with the use of Sobol’ points
is that one does not necessarily know in advance how many points will be
needed, when considering heterogeneous cancellation regions.

###  VI.3 Demonstration of Heterogeneous Cancellation Regions on the C5G7
Benchmark

We also tested the rejection-based sampling technique described in Section
V.3, for performing regional cancellation in cuboid regions which contain
multiple materials. Instead of using our modified version of the C5G7
benchmark, we have opted to use the original version with cylindrical fuel
pins [30], in combination with the the same $170\times 170\times 765$ mesh as
used in our previous simulations. A reference calculation was performed using
standard delta-tracking and obtained a multiplication factor of
$k_{\text{eff}}=1.18383\pm 0.00003$, which is in agreement with the reference
solution for the 3D version of the benchmark [30].

Figure 4: Histogram distribution of the Student $t$-variable, comparing the
flux computed from delta-tracking to the flux computed with negative weighted
delta-tracking. A normal distribution is plotted on top of the histogram for a
reference. Figure 5: Plots of the positive, negative, net, and total weights
as a function of generation, for the original C5G7 benchmark. Exact regional
cancellation with heterogeneous regions was used to perform the simulation.

When running the same simulation with negative-weighted delta-tracking, using
the same sampling cross sections as before, and exact cancellation, a
multiplication factor of $k_{\text{eff}}=1.18382\pm 0.00009$ was obtained,
which is in agreement with the delta-tracking value. A comparison of the two
estimations of the flux were also made, looking at the the Student
$t$-variable, which is defined as

$t_{i}=\frac{\varphi_{i,A}-\varphi_{i,B}}{\sqrt{\sigma^{2}_{i,A}+\sigma^{2}_{i,B}}}\text{,}$
(60)

where $\varphi_{i,x}$ is the average value of the flux in the $i$-th bin for
calculation $x$, and $\sigma_{i,x}$ is its standard error. For independent,
normally distributed variables with a large number of degrees of freedom, the
distribution of the $t$-variable should approach a normal distribution. A plot
of the empirical distribution of the $t$-variable is provided in Figure 4. We
excluded from the comparison all the bins where the flux was estimated to be
zero, or where the relative standard error was greater than 20%. This was done
in an effort to ensure that each bin was approximately normally distributed,
for the Student $t$-variable distribution assumptions to be reasonable. It is
clear from Figure 4 that the two flux estimates are in good agreement. The
$t$-variable distribution is not perfectly normal, which is to be expected, as
there are correlations between the scores in different flux bins. In general,
this is a very strong indication that our cancellation method has not imposed
any bias on the fission source, and that the method is still exact when
applied to heterogeneous cancellation regions.

The behavior of the total weight is shown in Fig. 5. While $W_{\text{tot}}$ is
larger than in the case of square fuel pins with homogeneous cancellation
region presented in Fig. 2, the behavior is in general similar. A large
$W_{\text{tot}}$ indicates that there are more negative particles, which will
increase the variance in scores. This is indeed the case, as our previous work
obtained an uncertainty for $k_{\text{eff}}$ of only $5\times 10^{-5}$ (albeit
for a slightly different problem) [22], while an uncertainty of $9\times
10^{-5}$ was obtained for this problem.

##  VII Conclusions

This work has leveraged the integral form of the Boltzmann transport equation
to provide a more in-depth mathematical analysis of the exact regional weight
cancellation technique, considerably expanding on previous works on the
subject [23, 22]. Not only has this formal approach given a much better
understanding as to the mechanics of regional cancellation in simplified
isotropic multi-group problems, but it has illuminated the non-trivial path to
performing exact cancellation in more complex problems, where scattering is
anisotropic, and the fission spectrum may depend on the incident energy of a
particle. The analysis highlights the fact that the implementation of exact
regional cancellation is rather straightforward in a multi-group Monte Carlo
code, but will be more difficult in a continuous-energy code, as one will need
access to conditional scattering distributions in the laboratory frame, which
are not always available for all reactions. The implementation of exact
regional cancellation in a continuous-energy code is thus a subject which will
require further research.

Additionally, a strategy to determine an optimal value of the cancellation
parameter $\beta$ for each particle undergoing cancellation has been
conceived. For each particle $k$ within a cancellation region, its optimal
cancellation parameter $\beta_{k}$ can be computed if both the average value
of the fission emission density in the region and the average value of the
inverse of the fission emission density in the region are known. Our previous
implementation required that cancellation regions be homogeneous and cubical,
which restricted its applicability to simple problems. Thanks to the
improvements in the optimization technique proposed in this work, both
requirements have been relaxed.

On a modified version of the C5G7 benchmark, our technique to optimize weight
cancellation was demonstrated to reduce the total weight in the simulation by
approximately 45%, when compared to using the minimum value of the fission
emission density in the region for $\beta$, as previously suggested in the
literature. In order to estimate the average fission emission density and the
average inverse of the fission emission density, a sampling approach has been
proposed, where the averages are estimated using the values of the fission
emission density at pseudo-random or quasi-random points. It was demonstrated
that the quasi-random Sobol’ sequence requires slightly fewer points than the
pseudo-random sequence to reach the asymptotic limit of the optimized
cancellation algorithm. As a comparison, the use of 3 Sobol’ points had very
similar performance to the use of 5 pseudo-random points for the modified
version of the C5G7 benchmark examined here. However, such results are likely
to be highly problem-dependent, and more systems should be analyzed to
ascertain what sort of performance improvements could be expected in general.
We also tested the use of heterogeneous cancellation regions on the original
C5G7 benchmark, with cylindrical fuel pins. No bias was observed in the
resulting fundamental eigenvalue, or flux tally.

##  Appendix A Optimization of $\Gamma_{2}^{*}$

We remind the reader of the definition of $\Gamma_{2}^{*}$:

$\Gamma^{*}_{2}=\sum_{k=1}^{N}\bigg{(}\frac{\expectationvalue{\zeta_{k}}-\beta_{k}}{\expectationvalue{\zeta_{k}}}w_{k}\bigg{)}^{2}+\bigg{(}\sum_{k=1}^{N}\frac{\beta_{k}}{\expectationvalue{\zeta_{k}}}w_{k}\bigg{)}^{2}.$
(61)

We optimize $\Gamma_{2}^{*}$ simultaneously for all particles by
differentiating with respect to $\beta_{j}$, and setting the partial
derivative equal to zero:

$\frac{\partial\Gamma_{2}^{*}}{\partial\beta_{j}}=-2\frac{\expectationvalue{\zeta_{j}}-\beta_{j}}{\expectationvalue{\zeta_{j}}^{2}}w^{2}_{j}+2\frac{w_{j}}{\expectationvalue{\zeta_{j}}}\sum_{k=1}^{N}\frac{\beta_{k}w_{k}}{\expectationvalue{\zeta_{k}}}=0.$
(62)

This may be simplified to

$-w_{j}+\frac{\beta_{j}w_{j}}{\expectationvalue{\zeta_{j}}}+\sum_{k=1}^{N}\frac{\beta_{k}w_{k}}{\expectationvalue{\zeta_{k}}}=0.$
(63)

On the left-hand-side, the second term matches the argument of the sum in the
third term. Summing over index $j$, we see that

$-W+\sum_{j=1}^{N}\frac{\beta_{j}w_{j}}{\expectationvalue{\zeta_{j}}}+N\sum_{k=1}^{N}\frac{\beta_{k}w_{k}}{\expectationvalue{\zeta_{k}}}=0.$
(64)

Here, we used the definition provided in Eq. (53). This allows us to isolate
the sum

$\sum_{k=1}^{N}\frac{\beta_{k}w_{k}}{\expectationvalue{\zeta_{k}}}=S^{*}=\frac{W}{N+1}.$
(65)

Applying this substitution to Eq. (63) while also using Eqs. (52) and (53), we
find that the optimized value of $\beta_{j}$ is

$\beta_{j}=\expectationvalue{\zeta_{j}}\left(1-\frac{W}{(N+1)w_{j}}\right)=\expectationvalue{\zeta_{j}}\left(1-\frac{S^{*}}{w_{j}}\right).$
(66)

##  Appendix B Optimization of $\expectationvalue{\Gamma_{2}}$

Substituting Eq. (48) into Eq. (54), and partially expanding the squared
terms, we see that

$\displaystyle\expectationvalue{\Gamma_{2}}$
$\displaystyle=\frac{\displaystyle\int_{\mathcal{R}}\bigg{[}\sum_{k=1}^{N}\bigg{(}1-\frac{2\beta_{k}}{\zeta_{k}}+\frac{\beta^{2}_{k}}{\zeta^{2}_{k}}\bigg{)}w^{2}_{k}+\sum_{k=1}^{N}\sum_{l=1}^{N}\frac{\beta_{k}}{\zeta_{k}}\frac{\beta_{l}}{\zeta_{l}}w_{k}w_{l}\bigg{]}\prod_{m=1}^{N}\zeta_{m}\differential
P^{\prime}_{m}}{\displaystyle\int_{\mathcal{R}}\prod_{n=1}^{N}\zeta_{n}\differential
P^{\prime}_{n}}$ (67)
$\displaystyle=\frac{\displaystyle\int_{\mathcal{R}}\bigg{[}\sum_{k=1}^{N}\bigg{(}1-\frac{2\beta_{k}}{\zeta_{k}}+\frac{\beta^{2}_{k}}{\zeta^{2}_{k}}\bigg{)}w^{2}_{k}+\sum_{k=1}^{N}\sum_{\begin{subarray}{c}l=1\\\
l\neq
k\end{subarray}}^{N}\frac{\beta_{k}}{\zeta_{k}}\frac{\beta_{l}}{\zeta_{l}}w_{k}w_{l}+\sum_{k=1}^{N}\frac{\beta_{k}^{2}}{\zeta_{k}^{2}}w_{k}^{2}\bigg{]}\prod_{m=1}^{N}\zeta_{m}\differential
P^{\prime}_{m}}{\displaystyle\prod_{n=1}^{N}\expectationvalue{\zeta_{n}}}$
(68)
$\displaystyle=\sum_{k=1}^{N}\bigg{(}1-\frac{2\beta_{k}}{\expectationvalue{\zeta_{k}}}+\frac{\beta^{2}_{k}}{\expectationvalue{\zeta_{k}}}\expectationvalue{\frac{1}{\zeta_{k}}}\bigg{)}w^{2}_{k}+\sum_{k=1}^{N}\sum_{\begin{subarray}{c}l=1\\\
l\neq
k\end{subarray}}^{N}\frac{\beta_{k}}{\expectationvalue{\zeta_{k}}}\frac{\beta_{l}}{\expectationvalue{\zeta_{l}}}w_{k}w_{l}+\sum_{k=1}^{N}\frac{\beta_{k}^{2}}{\expectationvalue{\zeta_{k}}}\expectationvalue{\frac{1}{\zeta_{k}}}w_{k}^{2}$
(69)
$\displaystyle=\sum_{k=1}^{N}\bigg{(}w^{2}_{k}-\frac{2\beta_{k}w^{2}_{k}}{\expectationvalue{\zeta_{k}}}\bigg{)}+\sum_{k=1}^{N}\frac{\beta^{2}_{k}w^{2}_{k}}{\expectationvalue{\zeta_{k}}}\bigg{(}2\expectationvalue{\frac{1}{\zeta_{k}}}-\frac{1}{\expectationvalue{\zeta_{k}}}\bigg{)}+\sum_{k=1}^{N}\sum_{l=1}^{N}\frac{\beta_{k}}{\expectationvalue{\zeta_{k}}}\frac{\beta_{l}}{\expectationvalue{\zeta_{l}}}w_{k}w_{l}.$
(70)

It is convenient to use the constant $c_{k}$, defined by Eq. (56), which may
be substituted into Eq. (70) to produce

$\expectationvalue{\Gamma_{2}}=\sum_{k=1}^{N}\bigg{(}w^{2}_{k}-\frac{2\beta_{k}w^{2}_{k}}{\expectationvalue{\zeta_{k}}}\bigg{)}+\sum_{k=1}^{N}\frac{\beta^{2}_{k}w^{2}_{k}}{c_{k}\expectationvalue{\zeta_{k}}^{2}}+\sum_{k=1}^{N}\sum_{l=1}^{N}\frac{\beta_{k}}{\expectationvalue{\zeta_{k}}}\frac{\beta_{l}}{\expectationvalue{\zeta_{l}}}w_{k}w_{l}.$
(71)

Now that all of the integrals have been simplified, we are left with
$\expectationvalue{\Gamma_{2}}$ as a function of $\beta_{k}$,
$\expectationvalue{\zeta_{k}}$, and $\expectationvalue{\frac{1}{\zeta_{k}}}$
$\forall k=1,\dots,N$. We now optimize $\expectationvalue{\Gamma_{2}}$ with
respect to the cancellation parameter $\beta_{j}$ by solving for

$\frac{\partial\expectationvalue{\Gamma_{2}}}{\partial\beta_{j}}=0.$ (72)

From Eq. (71), one may then proceed by solving

$\frac{\partial\expectationvalue{\Gamma_{2}}}{\partial\beta_{j}}=-\frac{2w_{j}^{2}}{\expectationvalue{\zeta_{j}}}+\frac{2\beta_{j}w_{j}^{2}}{c_{j}\expectationvalue{\zeta_{j}}^{2}}+\frac{2w_{j}}{\expectationvalue{\zeta_{j}}}\sum_{k=1}^{N}\frac{\beta_{k}w_{k}}{\expectationvalue{\zeta_{k}}}=0.$
(73)

Upon a division by $2w_{j}/\expectationvalue{\zeta_{j}}$ on both sides, we are
left with

$-w_{j}+\frac{\beta_{j}w_{j}}{c_{j}\expectationvalue{\zeta_{j}}}+\sum_{k=1}^{N}\frac{\beta_{k}w_{k}}{\expectationvalue{\zeta_{k}}}=0.$
(74)

It is possible to isolate the sum in the third term on the left-hand-side by
multiplying by $c_{j}$, and then summing over $j$:

$-\sum_{j=1}^{N}c_{j}w_{j}+\sum_{j=1}^{N}\frac{\beta_{j}w_{j}}{\expectationvalue{\zeta_{k}}}+\sum_{j=1}^{N}c_{j}\sum_{k=1}^{N}\frac{\beta_{k}w_{k}}{\expectationvalue{\zeta_{k}}}=0.$
(75)

We will now define

$S=\sum_{k=1}^{N}\frac{\beta_{k}w_{k}}{\expectationvalue{\zeta_{k}}},$ (76)

and substitute Eq. (76) into Eq. (75), allowing one to solve for $S$. Doing
so, one may obtain the result provided by Eq. (57). Now that the summation
term, $S$, can be computed without knowledge of $\beta_{j}$, we may substitute
Eq. (57) and Eq. (76) into Eq. (74), and solve for $\beta_{j}$, producing

$\beta_{j}=\expectationvalue{\zeta_{j}}c_{j}\bigg{(}1-\frac{S}{w_{j}}\bigg{)}.$
(77)

## References

  * Schneider  _et al._ [2016] D. Schneider, F. Dolci, F. Gabriel, J. Palau, M. Guillo, B. Pothet, P. Archier, K. Ammar, F. Auffret, R. Baron, _et al._ , in  _PHYSOR 2016_ (2016).
  * Collins  _et al._ [2016] B. Collins, S. Stimpson, B. W. Kelley, M. T. Young, B. Kochunas, A. Graham, E. W. Larsen, T. Downar, and A. Godfrey, J. Comput. Phys. 326, 612–628 (2016).
  * Rhodes  _et al._ [2006] J. Rhodes, K. Smith, and D. Lee, in  _PHYSOR 2006_ , Vol. 144 (2006).
  * Lux and Koblinger [1991] I. Lux and L. Koblinger,  _Monte Carlo Particle Transport Methods: Neutron and Photon Calculations_ (CRC Press, 1991).
  * Romano  _et al._ [2020] P. K. Romano, S. P. Hamilton, R. O. Rahaman, A. Novak, E. Merzari, S. M. Harper, P. C. Shriwise, and T. M. Evans, Nucl. Sci. Eng. 195, 1–21 (2020).
  * Mancusi  _et al._ [2022] D. Mancusi, M. Faucher, and A. Zoia, Eur. Phys. J. Plus 137, 127 (2022).
  * Yamamoto [2013] T. Yamamoto, Ann. Nucl. Energy 58, 72–79 (2013).
  * Booth [2003] T. E. Booth, Nucl. Sci. Eng. 143, 291–300 (2003).
  * Yamamoto [2012] T. Yamamoto, Ann. Nucl. Energy 47, 14–20 (2012).
  * Rouchon  _et al._ [2017] A. Rouchon, A. Zoia, and R. Sanchez, Ann. Nucl. Energy 102, 465–475 (2017).
  * Belanger  _et al._ [2020] H. Belanger, D. Mancusi, and A. Zoia, Eur. Phys. J. Plus 135, 877 (2020).
  * Arnow  _et al._ [1982] D. M. Arnow, M. H. Kalos, M. A. Lee, and K. E. Schmidt, J. Chem. Phys. 77, 5562–5572 (1982).
  * Assaraf  _et al._ [2007] R. Assaraf, M. Caffarel, and A. Khelif, J. of Phys. A: Math. Theor. 40, 1181 (2007).
  * Brun  _et al._ [2015] E. Brun, F. Damian, C. Diop, E. Dumonteil, F. Hugot, C. Jouanne, Y. Lee, F. Malvagi, A. Mazzolo, O. Petit, J. Trama, T. Visonneau, and A. Zoia, Ann. Nucl. Energy 82, 151–160 (2015).
  * Romano  _et al._ [2015] P. K. Romano, N. E. Horelik, B. R. Herman, A. G. Nelson, B. Forget, and K. Smith, Ann. Nucl. Energy 82, 90–97 (2015).
  * Leppänen  _et al._ [2015] J. Leppänen, M. Pusa, T. Viitanen, V. Valtavirta, and T. Kaltiaisenaho, Ann. Nucl. Energy 82, 142–150 (2015).
  * Goorley  _et al._ [2012] T. Goorley, M. James, T. Booth, F. Brown, J. Bull, L. J. Cox, J. Durkee, J. Elson, M. Fensin, R. A. Forster, J. Hendricks, H. G. Hughes, R. Johns, B. Kiedrowski, R. Martz, S. Mashnik, G. McKinney, D. Pelowitz, R. Prael, J. Sweezy, L. Waters, T. Wilcox, and T. Zukaitis, Nucl. Technol. 180, 298–315 (2012).
  * Woodcock  _et al._ [1965] E. R. Woodcock, T. Murphy, P. J. Hemmings, and T. C. Longworth,  _Techniques used in the GEM code for Monte Carlo neutronics calculations in reactors and other systems of complex geometry_ , Tech. Rep. (Argonne National Laboratory, 1965) ANL-7050.
  * Leppänen [2017] J. Leppänen, Ann. Nucl. Energy 105, 161–167 (2017).
  * Carter  _et al._ [1972] L. L. Carter, E. D. Cashwell, and W. M. Taylor, Nucl. Sci. Eng. 48, 403–411 (1972).
  * Legrady  _et al._ [2017] D. Legrady, B. Molnar, M. Klausz, and T. Major, Ann. Nucl. Energy 102, 116–123 (2017).
  * Belanger  _et al._ [2021a] H. Belanger, D. Mancusi, and A. Zoia, Phys. Rev. E 104, 015306 (2021a).
  * Booth and Gubernatis [2010] T. E. Booth and J. E. Gubernatis, Nucl. Sci. Eng. 165, 283–291 (2010).
  * Spanier and Gelbard [2008] J. Spanier and E. M. Gelbard,  _Monte Carlo Principles and Neutron Transport Problems_ (Dover Publications, 2008).
  * Coleman [1968] W. A. Coleman, Nucl. Sci. Eng. 32, 76–81 (1968).
  * Feller [1968] W. Feller,  _An Introduction to Probability Theory and its Applications_ , Vol. 1 (Wiley, New York, NY, 1968) p. 229.
  * Brunner and Brantley [2009] T. A. Brunner and P. S. Brantley, J. Comput. Phys. 228, 3882–3890 (2009).
  * Niederreiter [1978] H. Niederreiter, Bull. Am. Math. Soc. 84, 957–1041 (1978).
  * Belanger [2022] H. Belanger, 10.5281/zenodo.6546715 (2022), https://github.com/HunterBelanger/mgmc.
  * Lewis  _et al._ [2001] E. E. Lewis, M. A. Smith, N. Tsoulfanidis, G. Palmiotti, T. A. Taiwo, and R. N. Blomquist,  _Benchmark specification for Deterministic 2-D/3-D MOX fuel assembly transport calculations without spatial homogenisation (C5G7 MOX)_ , Tech. Rep. JT00105087 (NEA/NSC, 2001).
  * Zhang  _et al._ [2016] P. Zhang, H. Lee, and D. Lee, J. Comput. Phys. 305, 387–402 (2016).
  * Belanger  _et al._ [2021b] H. Belanger, D. Mancusi, and A. Zoia, in  _M &C 2021_ (2021).