# Exact weight cancellation in Monte Carlo eigenvalue transport problems

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

Random walks are frequently used as a model for very diverse physical
phenomena. The Monte Carlo method is a versatile tool for the study of the
properties of systems modelled as random walks. Often, each walker is
associated with a statistical weight, used in the estimation of observable
quantities. Weights are typically assumed to be positive; nonetheless, some
applications require the use of positive and negative weights or complex
weights, and often pose particular challenges with convergence. In this paper,
we examine such a case from the field of nuclear reactor physics, where the
negative particle weights prevent the power iteration algorithm from
converging on the sought fundamental eigenstate of the Boltzmann transport
equation. We demonstrate how the use of weight cancellation allows convergence
on the physical eigenstate. To this end, we develop a novel method to perform
weight cancellation in an exact manner, in three spatial dimensions. The
viability of this algorithm is then demonstrated on a reactor physics problem.

##  I Introduction

Many physical processes can be be represented by the random movement of
particles, or “walkers”, through phase space. Such phenomena include radiation
transport, propagation of active molecules in living bodies, or the spread of
epidemics [1, 2, 3, 4]. Quite often, it is desirable to compute the properties
of such systems with the aid of Monte Carlo simulations. This is certainly the
case in neutron transport, where Monte Carlo techniques are used to sample the
random walk process and thus solve the Boltzmann neutron transport equation
[5]. In this context, the Monte Carlo method may be preferred to (faster)
deterministic methods because it requires very few (if any) approximations:
Monte Carlo simulations are typically used to obtain reference solutions to
which deterministic solutions may be compared.

The random walkers have an associated statistical “weight”, which is typically
positive and real. This weight is used to estimate the observable quantities
which are sought in the simulation. Certain forms of transport problems,
however, require the use of positive and negative weights, or even complex
weights. Some examples include diffusion quantum Monte Carlo [6], the solution
of neutron noise equations in the frequency domain [7], the determination of
the second harmonic of the Boltzmann eigenvalue equation [8], and the
determination of transmittance in graphics rendering [9]. Monte Carlo
simulations that use negative or complex weights are notoriously difficult, as
they often do not converge to the desired solution [6, 8], or have a very high
variance in the observed quantities [7]. Most literature recognizes that
applying weight cancellation, where walkers carrying positive and negative
weights may annihilate with one another, is highly beneficial for the solution
of these problems [10, 8, 11, 7, 6].

Several methods for performing weight cancellation have been proposed. One
method, stemming from neutron transport, achieves weight cancellation by
pairing walkers individually (as opposed to pooling weights of walkers within
a defined region) [8]. While in theory this algorithm is exact (i.e. it does
not introduce any approximation), it is difficult to implement, and has a
quadratic computational complexity in the number of random walkers being
simulated. Another linear, exact regional cancellation algorithm has been
proposed [11], but is only valid in 1D geometries. One alternative approach
consists of using an approximate technique of averaging the weight of all
walkers within the same region [7]. This is easily applied to 3D, and has
linear complexity, but does not provide an exact solution. Finally, it must be
mentioned that not all random walker problems are posed in a manner which
makes cancellation difficult. In the full configuration interaction quantum
Monte Carlo technique, for example, cancellation is straightforward, as
walkers explore a discrete state space, as opposed to a continuous state
space; when two walkers of different signs land in the same state, they may
immediately annihilate [12].

In this paper, we review some of the existing weight cancellation methods and
propose a novel exact 3D regional cancellation algorithm, that has linear
complexity in the number of walkers, and that can be easily generalized to any
number of dimensions. The development of this method has been motivated by
certain neutron transport problems that require the use of negative walker
weights in conjunction with the power iteration technique (for obtaining the
dominant eigenstate), and that fail to converge to the correct solution
without weight cancellation; a simplified example of such a system is
presented in Section II. Subsequently, in Section III, we will develop a
mathematical model to explain why the power iteration fails in the presence of
negative weights, and why weight cancellation can resolve this problem. We
will demonstrate how the existing 1D regional cancellation algorithm allows
successful convergence of our 1D power iteration problem in Section IV, and we
will develop our novel, 3D version of the cancellation algorithm in Section V,
testing it on a reactor physics benchmark in Section VI. Finally, we will
present our conclusions Section VII.

##  II Power Iteration Debacle

In neutron transport, a main quantity of interest is the fundamental eigenpair
of the Boltzmann $\displaystyle k$-eigenvalue equation for neutrons, which we
write here in its one-speed form, as a function of position
$\displaystyle\bm{r}$ and direction $\displaystyle\hat{\bm{\Omega}}$:

$\displaystyle\displaystyle\bm{\hat{\Omega}}\cdot\nabla\varphi_{k}(\bm{r},\hat{\bm{\Omega}})+\Sigma_{t}(\bm{r})\varphi_{k}(\bm{r},\hat{\bm{\Omega}})=$
$\displaystyle\displaystyle\int_{4\pi}\Sigma_{s}(\bm{r},\hat{\bm{\Omega}}^{\prime}\rightarrow\hat{\bm{\Omega}})\varphi_{k}(\bm{r},\hat{\bm{\Omega}}^{\prime})\differential\hat{\bm{\Omega}}^{\prime}+$
$\displaystyle\displaystyle\frac{\nu(\bm{r})\Sigma_{f}(\bm{r})}{4\pi
k}\int_{4\pi}\varphi_{k}(\bm{r},\hat{\bm{\Omega}}^{\prime})\differential\hat{\bm{\Omega}}^{\prime}.$
(1)

Here $\displaystyle\varphi_{k}$ is the eigenfunction, $\displaystyle k$ is the
eigenvalue, $\displaystyle\Sigma_{t}$ is the total cross section,
$\displaystyle\Sigma_{s}$ is the scattering cross section,
$\displaystyle\Sigma_{f}$ is the fission cross section, and $\displaystyle\nu$
is the average number of new particles produced per fission. While it does not
explicitly appear in Eq. (1), there is an implicitly defined capture cross
section $\displaystyle\Sigma_{c}=\Sigma_{t}-\Sigma_{s}-\Sigma_{f}$, which
results in the death of neutrons. In general, Eq. (1) admits several
eigenvalue-eigenfunction pairs; the fundamental (largest) $\displaystyle
k$-eigenvalue is the multiplication factor, which is typically written as
$\displaystyle k_{\text{eff}}$. The corresponding eigenfunction
$\displaystyle\varphi_{0}$ is known as the angular neutron flux, and it
represents the average number of neutrons crossing a unit surface area per
unit time [5].

In order to obtain $\displaystyle k_{\text{eff}}$ and
$\displaystyle\varphi_{0}$, a numerical technique known as power iteration is
often employed [5, 13]. Although Eq. (1) is not written in such a form, we may
assume that it represents an eigenvalue problem of the type $\displaystyle
Lv_{0}=\lambda_{0}v_{0}$, with $\displaystyle v_{0}$ being the fundamental
eigenstate, and $\displaystyle\lambda_{0}$ being the fundamental eigenvalue.
In relation to Monte Carlo transport problems, the operator $\displaystyle L$
can be interpreted as the propagation of particles through a system from one
fission event to another. With power iteration, one may obtain the fundamental
state $\displaystyle v_{0}$ from any state $\displaystyle b$ such that the
inner product $\displaystyle\langle b,v_{0}\rangle\not=0$, by repeated
application of the operator $\displaystyle L$:

$\lim_{n\rightarrow\infty}\frac{L^{n}b}{\absolutevalue{L^{n}b}}=v_{0}.$ (2)

Therefore, starting from almost any initial state, it is possible to converge
to the fundamental mode [13].

We wish to use Monte Carlo to yield a solution to Eq. (1), in the sense that
the average density of walkers in phase space should satisfy Eq. (1), and
represent the fundamental eigenmode. In the context of Monte Carlo, the idea
of the power iteration technique needs to be adapted as follows. Starting from
an arbitrary set of walkers, we sample random walks for all walkers until
their death. Along the walk, a walker may be randomly killed if the magnitude
of its weight becomes too low, or split if the magnitude becomes too large
[14]. During the random walks of this first generation, new particles will be
born from fission. These fission particles are stored in a bank and are
attributed to the second generation. Once the first generation has finished,
the banked particles undergo population control, where those with small
weights may randomly be killed, those with large weights may be split into
multiple particles, and the net weight of all particles is normalised. This is
done to keep the net weight of all the particles at the beginning of a
generation constant, and to keep the number of particles in the simulation
constant on average [14]. This culled and normalized particle bank is
subsequently used as the source for the second generation. These particles
then undergo the random walk, producing fission particles which will belong to
the third generation. This application of the random walk mechanics on a
generation of particles may be continued indefinitely. After a number of a
generations, the positions (and directions) of the fission particles will
settle on an equilibrium distribution, representing the converged fission
source for the problem [14]. The combined random walks of all particles from
this converged fission source represent the dominant eigenfunction of Eq. (1).
With the fission source converged, observable quantities (such as the angular
flux and the multiplication factor) may be estimated in each generation. With
an estimation of the angular flux and the multiplication factor obtained by
each generation, an average may be calculated, although the estimation of the
uncertainty is not trivial [15].

Within a generation, the random walk process is sampled as follows. A particle
begins with an initial position and direction, which for the first generation
may be sampled from a somewhat arbitrary distribution. The first task is to
sample the distance the particle will fly before having a collision. In most
applications of Monte Carlo simulation to reactor physics problems, the total
cross section is usually assumed to be piece-wise constant within each
macroscopic geometric region composing the modeled system. To sample a flight
distance in a material region with a spatially constant cross section, a
random variable $\displaystyle\xi\sim\mathcal{U}[0,1)$ is drawn, and the
cumulative distribution function for the flight distance must be inverted to
obtain the distance to collision $\displaystyle d_{c}$111It is for this reason
that the underlying random walk is called “exponential flight” [16].:

$\xi=\int_{0}^{d_{c}}\Sigma_{t}\exp(-\Sigma_{t}s)\differential s\Rightarrow
d_{c}=-\frac{\ln(1-\xi)}{\Sigma_{t}}.$ (3)

The sampled distance is only valid in the given material region. If the
distance $\displaystyle d_{m}$ to the next material boundary is less than
$\displaystyle d_{c}$, the particle is only moved by $\displaystyle d_{m}$,
the cross section is updated, and a new distance to collision is sampled. Once
at a collision site, the collision mechanics may be simulated; the particle
can be captured with probability $\displaystyle\Sigma_{c}/\Sigma_{t}$,
whereupon the history is terminated; alternatively, scattering events may
occur with probability $\displaystyle\Sigma_{s}/\Sigma_{t}$, or the particle
may produce new fission particles with probability
$\displaystyle\Sigma_{f}/\Sigma_{t}$, which will be part of the next
generation [14]. The multiplication factor may be estimated in several
different manners: here we use the collision estimator. When a particle of
weight $\displaystyle w$ has a collision, we add the expected number of
fission neutrons produced per collision at that location, namely
$\displaystyle w\nu(\bm{r})\Sigma_{f}(\bm{r})/\Sigma_{t}(\bm{r})$, to an
accumulator which we shall refer to as $\displaystyle K$. Once the generation
has finished, the estimate for the multiplication factor for the generation
$\displaystyle G$ is

$k_{\text{eff}}^{(G)}=\frac{K}{W},$ (4)

$\displaystyle G$ being the generation number, and $\displaystyle W$ being the
sum of the weights of all the particles at the beginning of the generation
(which is a constant, because of the normalization procedure described above).
The accumulator $\displaystyle K$ must of course be reset to zero before each
new generation begins.

A more realistic description of the system would consist of relaxing the
hypothesis of piece-wise constant cross sections, especially in view of multi-
physics problems [17] where the cross sections for neutron transport depend on
complex space-dependent physical feedback mechanisms such as temperature and
material density fields. Recently, we examined different methods of sampling
the flight distance for the case of spatially varying cross sections [18].
Such methodologies are highly desirable for the next generation of Monte Carlo
transport codes, as they could allow for a better representation of the system
being simulated, reduce memory requirements, and pair well with the reactor
physics community’s goals of taking into account a multi-physics approach
[19]. In principle, sampling the distance to collision for spatially dependent
cross sections would require inverting the following equation for
$\displaystyle d_{c}$:

$\xi=\int_{0}^{d_{c}}\Sigma_{t}(s\hat{\bm{\Omega}}+\bm{r}_{0})\exp(\int_{0}^{s}\Sigma_{t}(u\hat{\bm{\Omega}}+\bm{r}_{0})\differential
u)\differential s.$ (5)

In practice, this is quite difficult to accomplish. Our previous work
considered several algorithms which could sample flight distances from
spatially continuous cross sections, without needing to directly invert Eq.
(5). In particular, we examined the traditional delta tracking [20, 21], and
the negative weighted delta tracking methods [22, 23].

In order to illustrate how these algorithms work, we revisit a simple one-
dimensional transport problem that was considered in our previous paper. In
the so-called ‘rod model’, particles may only move in the forward or backward
direction, along a line segment. For our purposes, this line segment has a
finite length, and it is possible for particles to leak out of either end of
the line. For our application we will set the boundaries of the segment at
$\displaystyle x=0$ and $\displaystyle x=2$. As an example of a space-
dependent cross section, we will use the broad-Gaussian cross section from our
earlier paper [18]. This cross section profile was chosen as it is not
monotonic, and better represents the idea that particles may see both an
increase and decrease in the cross section along their ray of flight. The
chosen $\displaystyle\Sigma_{t}(x)$ has the form

$\Sigma_{t}(x)=\sqrt{\frac{2}{\pi}}e^{-(x-1.23)^{2}}+0.1\text{ }\forall
x\in[0,2]\text{.}$ (6)

The scattering, absorption and fission probabilities are chosen to be
spatially constant, having values of $\displaystyle\Sigma_{s}/\Sigma_{t}=0.7$,
$\displaystyle\Sigma_{c}/\Sigma_{t}=0.1$, and
$\displaystyle\Sigma_{f}/\Sigma_{t}=0.2$. The average number of neutrons born
per fission is $\displaystyle\nu=2.5$; both scattering and fission are
isotropic (equal probability of emission forward or backward).

1 Sample uniform random variable $\displaystyle\xi_{1}$;

2 $\displaystyle d:=-\ln(\xi_{1})/\Sigma_{\text{smp}}$;

3 $\displaystyle\bm{r}_{1}:=d\hat{\bm{\Omega}}+\bm{r}_{0}$;

4 Sample random variable $\displaystyle\xi_{2}$;

5 if  _$\displaystyle\xi_{2} <q$_ then

6 $\displaystyle\displaystyle
w:=w\frac{\Sigma_{t}(\bm{r}_{1})}{q\,\Sigma_{\text{smp}}}$;

7 Perform real collision;

8

9else

10 $\displaystyle\displaystyle
w:=w\frac{1-\frac{\Sigma_{t}(\bm{r}_{1})}{\Sigma_{\text{smp}}}}{1-q}$;

11 Virtual collision; goto line 1;

12

13 end if

Algorithm 1 Negative weighted delta tracking

Transport is conducted using the negative weighted delta tracking (NWDT)
method [22, 23]. The basic procedures for this method are presented in Alg. 1.
It relies on the concept of real and virtual collisions. NWDT requires two
parameters: a sampling cross section $\displaystyle\Sigma_{\text{smp}}$, which
is used to sample the distance to a tentative collision site, and a
probability $\displaystyle q$, which is used to determine whether a collision
is real or virtual. At a real collision, the mechanics of a collision are used
to change the direction of the particle accordingly. A virtual collision does
not simulate the collision mechanics; instead, the direction of the particle
is left unchanged, and the distance to a new tentative collision site is
sampled, where the process begins again. Both
$\displaystyle\Sigma_{\text{smp}}$ and $\displaystyle q$ are allowed to be
functions of the position. According to the NWDT algorithm, the particle
weight changes sign when a virtual collision is sampled, and
$\displaystyle\Sigma_{t}>\Sigma_{\text{smp}}$. Here, we have chosen
$\displaystyle\Sigma_{\text{smp}}=0.85\Sigma_{t}(x=1.23)$, which leads to
particles changing sign within the region $\displaystyle 0.8\lesssim x\lesssim
1.66$. As for $\displaystyle q$, we have used $\displaystyle
q=\Sigma_{t}/(\Sigma_{t}-\absolutevalue{\Sigma_{\text{smp}}-\Sigma_{t}})$,
which corresponds to the strategy proposed by Carter, Cashwell and Taylor [24,
23]. This choice of $\displaystyle q$ was motivated by our previous work,
where the method demonstrated very reasonable performance compared to the
standard of delta tracking [18].

Delta tracking is a special case of NWDT, where the sampling cross section
$\displaystyle\Sigma_{\text{smp}}$ is taken to be a strict majorant cross
section $\displaystyle\Sigma_{\text{maj}}$, with
$\displaystyle\Sigma_{\text{maj}}\geq\Sigma_{t}$ and the acceptance
probability $\displaystyle q$ is taken equal to
$\displaystyle\Sigma_{t}/\Sigma_{\text{maj}}$ [20, 21]222Delta tracking is
also sometimes referred to as self-scattering in the electron transport
community [25], as thinning [26], or as the null-collisions method [27].. With
these choices, particle weights never change sign; however, insisting on the
use of delta tracking can raise a number of practical problems, because it
could be very challenging to obtain a majorant when using spatially continuous
cross sections in a more realistic context. If one were to conduct a
simulation where both the isotopic density and the temperature of materials
varied spatially, there is, in general, no way to exactly determine the
majorant cross section for the system. It is certainly possible to probe the
phase space of the problem, testing the cross section at each point. However,
if the selected sampling cross section is not in fact a majorant, it will be
impossible to handle collisions for which
$\displaystyle\Sigma_{t}>\Sigma_{\text{smp}}$; even if such an event does not
occur, the results will be slightly biased. Conversely, if
$\displaystyle\Sigma_{\text{smp}}$ is chosen by applying a large safety margin
to the maximum known cross section (in an effort to ensure underestimation of
the majorant does not occur), the method becomes quite inefficient, as many
unnecessary virtual collisions will occur [21]. NWDT tolerates
underestimations in the majorant cross section, and also allows one to avoid
gross overestimations of the majorant which reduce efficiency.

With NWDT, not only may there be negatively weighted particles in the system,
but the weights may also change in magnitude. These two traits can lead to an
increase in variance and simulation time. To mitigate these effects, roulette
is used on particles with $\displaystyle\absolutevalue{w}<0.6$, with the
survival weight being $\displaystyle w=\pm 1$ (ensuring that the particle
keeps its initial sign). Particles are also split if
$\displaystyle\absolutevalue{w}\geq 2$.

In our previous work, we neglected fission. In view of testing NWDT in the
framework of power iteration for eigenvalue problems, fission has been added
to our model. The number $\displaystyle n$ of new fission neutrons generated
at any collision is taken to be

$n=\left\lfloor\absolutevalue{w}\frac{\nu\Sigma_{f}}{\Sigma_{t}}\cdot\frac{1}{k_{\text{eff}}^{(G-1)}}+\xi\right\rfloor,$
(7)

where $\displaystyle\xi\sim\mathcal{U}[0,1)$ is a uniform random variable, and
$\displaystyle k_{\text{eff}}^{(G-1)}$ is the estimated value of
$\displaystyle k_{\text{eff}}$ for the previous generation of particles.
Without dividing by $\displaystyle k_{\text{eff}}^{(G-1)}$, the number of
particles in the simulation would either increase exponentially if the system
is super-critical ($\displaystyle k_{\text{eff}}>1$), or decrease
exponentially if the system is sub-critical ($\displaystyle
k_{\text{eff}}<1$). This algorithm closely follows standard methods applied in
Monte Carlo codes [28, 29], with slight modifications to accommodate fission
with negative particles. All new fission particles are born with a weight
$\displaystyle w=\pm 1$, keeping the sign of the weight of the particle
inducing fission. Between particle generations, all particle weights are
multiplied by a normalization coefficient to ensure that the net weight of the
system (sum of all particle weights) is always a constant value.

Figure 1: The number of positive and negative particles per generation, when
using the population control scheme provided in Eq. (7). Figure 2: The weight
of positive and negative particles per generation, when using the population
control scheme provided in Eq. (7).

Figure 1 shows the number of positive and negative particles per generation,
when starting with $\displaystyle 10^{6}$ particles uniformly distributed
within the rod. Also shown is the net number of particles $\displaystyle
N_{\text{net}}=N_{+}-N_{-}$, and the total number of particles $\displaystyle
N_{\text{tot}}=N_{+}+N_{-}$. Similarly, Figure 2 shows the total positive
weight $\displaystyle W_{+}$ (the sum of the weights of all positive
particles), the total negative weight $\displaystyle W_{-}$ (the sum of the
magnitude of the weights of all negative particles), together with the net
weight $\displaystyle W_{\text{net}}=W_{+}-W_{-}$ and the total weight
$\displaystyle W_{\text{tot}}=W_{+}+W_{-}$. The number of particles is
somewhat less meaningful than the weight, as one can use methods such as
weight combing [30] to adjust how many particles are followed in the
simulation, but this says nothing about the magnitude of their weight. Weight
is the more natural quantity to examine, as its behavior is not controlled by
the number of particles. It is instructive, however, to consider the interplay
between weight and number of particles. While Eq. (7) keeps the net number of
particles stable as expected, the total number of particles is left to
increase exponentially. Eventually, this increase in the number of particles
causes each generation to take a longer amount of time to process, and the
memory requirements increase as well. At some point, the computer memory is
overwhelmed, and the simulation fails. In an attempt to remedy this problem,
we modified Eq. (7) by replacing $\displaystyle k_{\text{eff}}^{(G-1)}$; using
the value of $\displaystyle k_{\text{eff}}$ for the previous generation keeps
the net number of particles constant on average, so a new quantity
$\displaystyle k_{\text{tot}}$ was defined, representing the the increase in
the total weight $\displaystyle W_{\text{tot}}$ of the simulation. The
quantity $\displaystyle k_{\text{tot}}$ may be estimated in a similar manner
to the multiplication factor, accumulating the collision estimator
$\displaystyle\absolutevalue{w}\nu(\bm{r})\Sigma_{f}(\bm{r})/\Sigma_{t}(\bm{r})$.
Using $\displaystyle k_{\text{tot}}^{(G-1)}$ in Eq. (7), we obtain

$n=\left\lfloor\absolutevalue{w}\frac{\nu\Sigma_{f}}{\Sigma_{t}}\cdot\frac{1}{k_{\text{tot}}^{(G-1)}}+\xi\right\rfloor.$
(8)

Figure 3: The number of positive and negative particles per generation, when
using the population control scheme provided in Eq. (8). Figure 4: The weight
of positive and negative particles per generation, when using the population
control scheme provided in Eq. (8).

The effects on the particle populations from using this normalization
technique are presented in Figure 3 and Figure 4. While the total number of
particles indeed remains constant on average, the net number of particles now
decreases exponentially. The net weight of course remains constant because we
are normalizing it at the end of the generation. With fewer net particles to
maintain the entire initial weight of the system, the magnitude of the weight
of the particles increases drastically. This is indicated by the total weight
of the system continuing to increase exponentially in Figure 4. This leads to
near equal quantities of positive and negative particles, each with very large
weight magnitudes, causing large fluctuations in the flux estimators, and
therefore a larger variance in scores.

Other population control mechanisms were tested as well, hoping that they
might stabilize the system. One such method was particle combing [30],
conducted between each generation of particles. While combing did reduce the
total number of particles going into each generation, the weight magnitude of
each particle was quite large. Once they entered the simulation, the particles
were split, causing an abrupt increase in memory consumption. After several
generations, the simulation is killed by the operating system due to the large
memory consumption.

Although the rod model example is admittedly simple, similar results have been
obtained when using the NWDT algorithm in combination with power iteration on
more realistic reactor physics problems in 3D, and with energy-dependent cross
sections (see Sec. VI). It is clear that the power iteration with NWDT cannot
be used to estimate the equilibrium distribution of walkers, and therefore the
solution to Eq. (1). It is therefore pertinent to better understand why the
power iteration is failing, and how the problem might be addressed.

##  III Necessity of Weight Cancellation

###  III.1 The failure of power iteration

The Boltzmann equation presented in Eq. (1) is not capable of describing the
interaction between the positive and negative particle populations in a
simulation, only their combined average angular flux. It is fairly
straightforward, however, to propose a formulation of the transport equation
which is capable of representing both positively and negatively weighted
particles. Considering Alg. 1, outlining the NWDT method, we may write a set
of coupled transport equations, one for the angular flux
$\displaystyle\varphi_{+}$ of the positive particles, and a second for the
angular flux $\displaystyle\varphi_{-}$ of the negative particles, resulting
in the following coupled equations:

$\displaystyle\displaystyle\hat{\bm{\Omega}}\cdot\nabla\varphi_{\pm}+\Sigma_{\text{smp}}\varphi_{\pm}=\mathcal{S}\varphi_{\pm}+\frac{1}{k}\mathcal{F}\varphi_{\pm}+$
$\displaystyle\displaystyle\Delta(\Sigma_{\text{smp}}-\Sigma_{t})\,\varphi_{\pm}+\Delta(\Sigma_{t}-\Sigma_{\text{smp}})\,\varphi_{\mp}.$
(9)

Here, $\displaystyle\mathcal{S}$ is defined as

$\displaystyle\displaystyle\mathcal{S}\varphi=\int_{4\pi}\Sigma_{s}(\bm{r},\hat{\bm{\Omega}}^{\prime}\rightarrow\hat{\bm{\Omega}})\varphi(\bm{r},\hat{\bm{\Omega}}^{\prime})\differential\hat{\bm{\Omega}}^{\prime}$
(10)

and $\displaystyle\mathcal{F}$ is defined as

$\displaystyle\displaystyle\mathcal{F}\varphi=\frac{\nu(\bm{r})\Sigma_{f}(\bm{r})}{4\pi}\int_{4\pi}\varphi(\bm{r},\hat{\bm{\Omega}}^{\prime})\differential\hat{\bm{\Omega}}^{\prime}.$
(11)

We have also made use of the function

$\Delta(x)=\begin{cases}x&x>0\\\ 0&x\leq 0\end{cases}.$ (12)

A rigorous derivation of this statement is given in the Appendix. There is a
subtlety in this formulation: instead of having one species of particles which
may have a weight that is positive or negative, we now have two species
(positive and negative), and both are represented by particles having a
strictly positive weight. In this framework, our unknown is now the set of
eigenfunctions

$\displaystyle\displaystyle\bm{\zeta}=\begin{bmatrix}\varphi_{+}\\\
\varphi_{-}\end{bmatrix},$ (13)

which is the solution of the generalized eigenvalue problem

$\bm{A}\bm{\zeta}=\frac{1}{k}\bm{F}\bm{\zeta},$ (14)

operator $\displaystyle\bm{A}$ being

$\bm{A}=\begin{bmatrix}\mathcal{A}_{11}&\mathcal{A}_{12}\\\
\mathcal{A}_{12}&\mathcal{A}_{11}\end{bmatrix}\text{,}$ (15)

with

$\begin{gathered}\mathcal{A}_{11}=\hat{\bm{\Omega}}\cdot\nabla+\Sigma_{\text{smp}}-\mathcal{S}-\Delta(\Sigma_{\text{smp}}-\Sigma_{t})\text{,}\\\
\mathcal{A}_{12}=-\Delta(\Sigma_{t}-\Sigma_{\text{smp}}),\end{gathered}$ (16)

and $\displaystyle\bm{F}$ being defined as

$\bm{F}=\begin{bmatrix}\mathcal{F}&&0\\\ 0&&\mathcal{F}\end{bmatrix}.$ (17)

If we use $\displaystyle\mathcal{V}$ to denote the vector space of the
physical flux, then $\displaystyle\bm{\zeta}$ is an element of the expanded
vector space $\displaystyle\mathcal{V}\times\mathcal{V}$.

Given the previous definitions, the physical flux is interpreted to be

$\varphi=\varphi_{+}-\varphi_{-}$ (18)

and may be retrieved through the application of
$\displaystyle\bm{D}:\mathcal{V}\times\mathcal{V}\rightarrow\mathcal{V}$,
defined as

$\bm{D}=\begin{bmatrix}I&&-I\end{bmatrix}$ (19)

and whose action is

$\bm{D}\begin{bmatrix}f\\\ g\end{bmatrix}=f-g.$ (20)

The eigenvalue equation presented in Eq. (14) has several interesting
properties, which may be elegantly outlined by the introduction of the parity
operator, namely

$\bm{P}=\begin{bmatrix}0&&I\\\ I&&0\end{bmatrix}.$ (21)

Eq. (14) is invariant under the action of the parity operator: that is, if
$\displaystyle\bar{\bm{\zeta}}$ is a solution of Eq. (14) with eigenvalue
$\displaystyle\bar{k}$, then so is $\displaystyle\bm{P}\bar{\bm{\zeta}}$. It
is easily verifiable that $\displaystyle\bm{P}\bm{A}\bm{P}=\bm{A}$ and
$\displaystyle\bm{P}\bm{F}\bm{P}=\bm{F}$, and therefore

$\bm{A}(\bm{P}\bar{\bm{\zeta}})=\frac{1}{\bar{k}}\bm{F}(\bm{P}\bar{\bm{\zeta}}).$
(22)

Thus, if the eigenvalue $\displaystyle\bar{k}$ is non-degenerate, we must have
$\displaystyle\bm{P}\bar{\bm{\zeta}}=\pm\bar{\bm{\zeta}}$. This partitions the
eigenstates into two sets, depending on the sign of their eigenvalue with
respect to $\displaystyle\bm{P}$. The odd eigenstates have the form

$\bm{\zeta}_{o}=\begin{bmatrix}\varphi\\\ -\varphi\end{bmatrix}.$ (23)

By inserting this ansatz into Eq. (14), we can verify that
$\displaystyle\varphi$ must solve the physical eigenvalue equation

$\bm{\hat{\Omega}}\cdot\nabla\varphi+\Sigma_{t}\varphi=\mathcal{S}\varphi+\frac{1}{k}\mathcal{F}\varphi,$
(24)

which is just Eq. (1). We denote the eigenfunctions and eigenvalues of this
equation as $\displaystyle\varphi_{i}$ and $\displaystyle k_{\varphi,i}$,
respectively, with $\displaystyle i\in\\{0,1,\ldots\\}$. The even eigenstates
have the form

$\bm{\zeta}_{e}=\begin{bmatrix}\eta\\\ \eta\end{bmatrix},$ (25)

where $\displaystyle\eta$ must solve the modified Boltzmann equation

$\bm{\hat{\Omega}}\cdot\nabla\eta+\Sigma_{t,\eta}\eta=\mathcal{S}\eta+\frac{1}{k}\mathcal{F}\eta,$
(26)

which is verified by substitution into Eq. (14). Here we have defined

$\Sigma_{t,\eta}=\Sigma_{\text{smp}}-\absolutevalue{\Sigma_{\text{smp}}-\Sigma_{t}}.$
(27)

We denote the eigenfunctions and eigenvalues of Eq. (26) as
$\displaystyle\eta_{i}$ and $\displaystyle k_{\eta,i}$, respectively, with
$\displaystyle i\in\\{0,1,\ldots\\}$.

Since we always have $\displaystyle\Sigma_{t,\eta}\leq\Sigma_{t}$, the
equation for $\displaystyle\eta$ has the same scattering and fission terms as
Eq. (1), but capture has been decreased, namely

$\displaystyle\displaystyle\Sigma_{c,\eta}=\Sigma_{t,\eta}-\Sigma_{s}-\Sigma_{f}\leq\Sigma_{t}-\Sigma_{s}-\Sigma_{f}.$
(28)

Thus, on physical grounds, we correspondingly expect a larger dominant
eigenvalue:

$k_{\eta,0}\geq k_{\varphi,0}=k_{\text{eff}}.$ (29)

Therefore, the dominant eigenvalue of Eq. (14) is not $\displaystyle
k_{\text{eff}}$, but the nonphysical eigenvalue $\displaystyle k_{\eta,0}$.
Applying the power iteration method to Eq. (14) will result in convergence
towards the latter eigenstate. Additionally, the equilibrium distribution is
even under the exchange of positive and negative particles, i.e. it contains
the same amount of positive and negative particles, and thus zero net
particles. Converging to a state which has zero net particles is incompatible
with traditional population control mechanisms (such as combing, or
normalizing by $\displaystyle k_{\text{eff}}$), which are designed to keep the
net weight constant, leading to the divergence in particle populations which
was observed in Sec. II.

It is worth stressing the similarity between our analysis and the study by
Spencer et al. [12] on the origin of the sign problem in full configuration
interaction quantum Monte Carlo without weight cancellation.

###  III.2 Modeling weight cancellation

Previous investigations have shown that weight cancellation can be very
effective in dealing with Monte Carlo problems using particles of both
positive and negative weights [31, 8, 6, 12]. In view of these considerations,
in the following we will formally address the effect of cancellation on Eq.
(14). Weight cancellation can be modeled in the following manner: we start
with some function of the phase space variables
$\displaystyle(\bm{r},\hat{\bm{\Omega}})$, i.e. an element of
$\displaystyle\mathcal{V}$. This function may be embedded in
$\displaystyle\mathcal{V}\times\mathcal{V}$ by applying a suitable mapping
$\displaystyle\bm{E}:\mathcal{V}\rightarrow\mathcal{V}\times\mathcal{V}$.
There is some latitude in the definition of $\displaystyle\bm{E}$; one
possible definition is

$\bm{E}_{1}f=\begin{bmatrix}f\\\ 0\end{bmatrix},$ (30)

but equally valid choices could be

$\displaystyle\displaystyle\bm{E}_{0}f=$
$\displaystyle\displaystyle\begin{bmatrix}0\\\ -f\end{bmatrix}$ (31)
$\displaystyle\displaystyle\bm{E}_{1/2}f=$
$\displaystyle\displaystyle\begin{bmatrix}f/2\\\ -f/2\end{bmatrix}$ (32)

or the non-linear variant

$\displaystyle\displaystyle\tilde{\bm{E}}f=$
$\displaystyle\displaystyle\begin{bmatrix}\max(f,0)\\\
-\min(f,0)\end{bmatrix}.$ (34)

The only property that we require of $\displaystyle\bm{E}$ is that it should
be right-inverse to $\displaystyle\bm{D}$, viz.

$\bm{D}\bm{E}=\bm{I}.$ (35)

This property expresses the fact that lifting a function from
$\displaystyle\mathcal{V}$ into $\displaystyle\mathcal{V}\times\mathcal{V}$
(the action of $\displaystyle\bm{E}$) followed by collapsing back into the
space $\displaystyle\mathcal{V}$ of physical fluxes (the action of
$\displaystyle\bm{D}$) should not change the function we started with. It is
easily verifiable that the previously proposed definitions of
$\displaystyle\bm{E}$ satisfy this property. In matrix notation, the operators
$\displaystyle\bm{E}_{1}$, $\displaystyle\bm{E}_{0}$, and
$\displaystyle\bm{E}_{1/2}$ are given by the general formula

$\displaystyle\displaystyle\bm{E}_{z}=\begin{bmatrix}zI\\\
(z-1)I\end{bmatrix}.$ (36)

It is not possible to express $\displaystyle\tilde{\bm{E}}$ as a matrix, as it
is non-linear.

Let us now consider the operator product with the opposite ordering,
$\displaystyle\bm{C}=\bm{E}\bm{D}$. This operator,
$\displaystyle\bm{C}:\mathcal{V}\times\mathcal{V}\rightarrow\mathcal{V}\times\mathcal{V}$
is a projector, by virtue of Eq. (35):

$\bm{C}^{2}=\bm{E}\bm{D}\bm{E}\bm{D}=\bm{E}(\bm{D}\bm{E})\bm{D}=\bm{E}\bm{D}=\bm{C}.$
(37)

By construction, the null space of $\displaystyle\bm{C}$ coincides with the
null space of $\displaystyle\bm{D}$, which is the space of even (unphysical)
vectors:

$\bm{C}\begin{bmatrix}\eta\\\ \eta\end{bmatrix}=0.$ (38)

Therefore, the operator $\displaystyle\bm{C}$ can be regarded as a model for
perfect cancellation. We have shown that the even vectors represent the
solutions of the nonphysical Boltzmann equation in Eq. (26). This is supported
by the explicit form of $\displaystyle\bm{C}$ which follows from the previous
definitions for $\displaystyle\bm{E}$:

$\displaystyle\displaystyle\bm{C}_{z}\begin{bmatrix}f_{+}\\\
f_{-}\end{bmatrix}=$
$\displaystyle\displaystyle\bm{E}_{z}\bm{D}\begin{bmatrix}f_{+}\\\
f_{-}\end{bmatrix}=$
$\displaystyle\displaystyle\begin{bmatrix}z(f_{+}-f_{-})\\\
(1-z)(f_{+}-f_{-})\end{bmatrix}$ (39)
$\displaystyle\displaystyle\tilde{\bm{C}}\begin{bmatrix}f_{+}\\\
f_{-}\end{bmatrix}=$
$\displaystyle\displaystyle\tilde{\bm{E}}\bm{D}\begin{bmatrix}f_{+}\\\
f_{-}\end{bmatrix}=$
$\displaystyle\displaystyle\begin{bmatrix}\max(f_{+}-f_{-},0)\\\
-\min(f_{+}-f_{-},0)\end{bmatrix}.$ (40)

Note that only $\displaystyle\tilde{\bm{C}}$ can guarantee that both vector
components are non-negative, but it requires that both $\displaystyle\bm{C}$
and $\displaystyle\bm{E}$ be non-linear. In matrix notation,
$\displaystyle\bm{C}_{z}$ reads

$\bm{C}_{z}=\begin{bmatrix}zI&-zI\\\ (z-1)I&-(z-1)I\end{bmatrix}.$ (41)

Going back to our eigenvalue equation, Eq. (14) may be modified to include the
application of cancellation on the fission source:

$\bm{A}\bm{\zeta}=\frac{1}{k}\bm{C}\bm{F}\bm{\zeta}.$ (42)

If $\displaystyle\bm{\zeta}$ is even, $\displaystyle\bm{F}\bm{\zeta}$ must
also be even as $\displaystyle\commutator{\bm{F}}{\bm{P}}=0$. Since
$\displaystyle\bm{C}$ maps even vectors to $\displaystyle 0$, the cancellation
operator in Eq. (42) causes the nonphysical eigenmodes corresponding to Eq.
(26) to vanish.

At this point, it remains to be shown that cancellation does not perturb the
physical eigenvalues of Eq. (14), that were initially associated with the odd
eigenstates. To accomplish this, we start with Eq. (42), and apply
$\displaystyle\bm{D}$ from the left on both sides, as this is the operator
which performs the mapping
$\displaystyle\mathcal{V}\times\mathcal{V}\rightarrow\mathcal{V}$:

$\bm{D}\bm{A}\bm{\zeta}=\frac{1}{k}\bm{D}\bm{C}\bm{F}\bm{\zeta}.$ (43)

From Eq. (35), we may substitute $\displaystyle\bm{D}\bm{C}=\bm{D}$. Using
Eqs. (13), (15), and (18), this then simplifies to

$(\mathcal{A}_{11}-\mathcal{A}_{12})\varphi=\frac{1}{k}\mathcal{F}\varphi.$
(44)

Eq. (12) indicates that $\displaystyle\Delta(x)-\Delta(-x)=x$, which, when
combined with Eq. (16), yields

$\big{[}\hat{\bm{\Omega}}\cdot\nabla+\Sigma_{t}\big{]}\varphi=\mathcal{S}\varphi+\frac{1}{k}\mathcal{F}\varphi.$
(45)

This is the physical Boltzmann transport equation. Therefore, cancellation
leaves the physical eigenstates unchanged, so long as Eq. (35) holds true,
regardless of the form of $\displaystyle\bm{E}$.

In summary, the introduction of cancellation conserves all the (physical and
nonphysical) eigenstates of Eq. (9); however, the nonphysical part of the
spectrum is relegated to zero, while the physical eigenvalues are unperturbed
and dominate the power iteration.

###  III.3 Deterministic proof of concept

In order to better understand the effects of cancellation, and how it
positively affects the convergence of the power iteration, it is fruitful to
visualize the spectrum of the associated equations. For this purpose, we have
written a deterministic solver that computes the full spectrum of eigenvalues
and eigenfunctions (using the LAPACK library [32]) for a discretized version
of the one-dimensional system introduced in Sec. II. A deterministic solver
was necessary, since Monte Carlo can only be used to estimate the fundamental
mode by power iteration. This solver uses a finite-difference method, with
cross section values being taken at the midpoint of each spatial bin; both the
positive and negative particle fluxes are explicitly considered.

Before considering cancellation, it is helpful to observe the behaviour of the
eigenvalue spectrum of Eq. (14) as a function of
$\displaystyle\Sigma_{\text{smp}}/\Sigma_{\text{maj}}$. This quantity roughly
measures of the amount of negative weight that is produced by NWDT; the
smaller the ratio is, the more negative weight is introduced. The spectrum is
the union of the sets of eigenvalues $\displaystyle k_{\varphi,i}$ and
$\displaystyle k_{\eta,i}$ and is presented in Figure 5. Naturally, the
spectrum associated with the physical Boltzmann equation does not change with
the ratio of $\displaystyle\Sigma_{\text{smp}}/\Sigma_{\text{maj}}$. The
eigenvalues of the nonphysical Boltzmann equation diverge from the physical
values as this ratio decreases from unity. The dominant nonphysical eigenvalue
$\displaystyle k_{\eta,0}$ continually increases above $\displaystyle
k_{\varphi,0}$, while the other depicted eigenvalues decrease from their
physical counterparts. In Eq. (29), we heuristically demonstrated that
$\displaystyle k_{\eta,0}\geq k_{\text{eff}}$, and this is indeed the case in
Figure 5. Our heuristics however, are not able to make any remarks as to the
behavior of the higher eigenvalues, and while for this system $\displaystyle
k_{\eta,i}$ for $\displaystyle i>0$ always appears to decrease when
$\displaystyle\Sigma_{smp}/\Sigma_{maj}$ decreases, we have observed other
systems which do not exhibit this behavior.

Figure 5: Behavior of the eigenvalue spectrum as a function of the ratio of
the sampling cross section to the majorant cross section. Eigenvalues
associated with the physical Boltzmann equation have a solid line, while those
associated with the nonphysical Boltzmann equation have a dotted line.

When the cancellation operator $\displaystyle\bm{C_{1}}$ is added in
accordance with Eqs. (41) and (42), we find that cancellation indeed
suppresses the even eigenstates, leaving only the odd ones associated with the
physical solutions, so long as
$\displaystyle\Sigma_{\text{smp}}/\Sigma_{\text{maj}}>0$. This is the case for
a cancellation operator which is “perfect”, in that it is able to conduct
cancellation in a manner which always completely neutralizes the negative
particle population. While this is of course desired, and possible to
implement in a deterministic solver, we will later show in Sec. IV that 100%
cancellation efficiency is not necessarily achievable in a Monte Carlo
approach. To mimic this fact, we introduce an imperfect cancellation operator
$\displaystyle\bm{C}^{(\alpha)}$:

$\displaystyle\displaystyle\bm{C}^{(\alpha)}=\alpha\bm{C}_{1}+(1-\alpha)\bm{I}.$
(46)

Here $\displaystyle\alpha$ represents the cancellation efficiency: when
$\displaystyle\alpha=1$ there is perfect cancellation, and when
$\displaystyle\alpha=0$ there is no cancellation. This is a highly idealized
approach to model imperfect cancellation, and is not necessarily a faithful
model of the effect of cancellation in a Monte Carlo setting.

Figure 6: Behaviour of the eigenvalue spectrum as a function of cancellation
efficiency $\displaystyle\alpha$, for the case of
$\displaystyle\Sigma_{\text{smp}}/\Sigma_{\text{maj}}=0.7$. Eigenvalues
associated with the physical Boltzmann equation have a solid line, while those
associated with the nonphysical Boltzmann equation have a dotted line.

Figure 6 presents the effects of cancellation with varying levels of
efficiency, for the case of
$\displaystyle\Sigma_{\text{smp}}/\Sigma_{\text{maj}}=0.7$. For values of
$\displaystyle\alpha$ between 1 and approximately 0.21, the cancellation
operator remains efficient enough that the physical eigenstate is the dominant
one. A critical point is reached near $\displaystyle\alpha\approx 0.21$, where
the eigenvalues of the physical and nonphysical systems are equal. Further
reducing the efficiency of cancellation leads to the nonphysical eigenstate
being the dominant one. This example indicates that for our system, there is a
_minimum amount of weight cancellation_ required in order for power iteration
to converge on the fundamental physical eigenstate. This is likely true for
all systems which can be described by coupled transport equations for positive
and negative particles.

Even in settings where enough cancellation is present to make the
$\displaystyle k_{\varphi,0}$ mode dominant, cancellation might still be
unable to reduce $\displaystyle k_{\eta,0}$ to a level that allows for rapid
statistical convergence towards the fundamental mode. To clarify this, we
introduce the dominance ratio, which is defined as the ratio of the second-
largest eigenvalue to the largest one. The dominance ratio is equal to
$\displaystyle k_{\eta,0}/k_{\varphi,0}$ for values of $\displaystyle\alpha$
where $\displaystyle k_{\varphi,0}>k_{\eta,0}>k_{\varphi,1}$. If instead
$\displaystyle\alpha$ is such that $\displaystyle k_{\eta,0}<k_{\varphi,1}$,
then the dominance ratio is $\displaystyle k_{\varphi,1}/k_{\varphi,0}$, which
is the dominance ratio of the physical system. The dominance ratio provides an
indication of the rate of convergence of the power iteration. When it is very
close to unity, more iterations are required to converge on the dominant
eigenvalue. Should $\displaystyle\alpha$ be sufficiently large to ensure
$\displaystyle k_{\eta,0}<k_{\varphi,1}$, the convergence rate of the problem
will no longer be bound by the efficiency of cancellation, but by the physical
properties of the system being examined.

This model for cancellation in a deterministic context provides valuable
insight as to the behavior of this coupled system of positive and negative
particles, and to the possible behavior of implementing cancellation of
particle weights in a Monte Carlo context. It is for this reason that we do
not go beyond the provided surface-level analysis of the effects of
cancellation efficiency or choice of $\displaystyle\Sigma_{\text{smp}}$ on the
convergence of the deterministic model. We will now continue by discussing
techniques of weight cancellation in Monte Carlo simulations.

##  IV An Exact Regional Cancellation Scheme for 1D Problems

As mentioned in the introduction, several cancellation strategies have been
proposed in the past; one of these, devised by Booth and Gubernatis, is exact
in 1D geometries [11]. This method works by partitioning all fissile domains
of the problem domain into regions. Between each generation, weight
cancellation amongst all new fission particles born in the region occurs. An
important feature of this algorithm is that it has linear computational
complexity with the number of particles partaking in cancellation, making it a
good candidate for inclusion in general-purpose Monte Carlo transport codes.
We have therefore chosen to focus on this algorithm. We will now provide a
brief overview of how Booth and Gubernatis’ method works, as it is essential
for understanding our proposed 3D algorithm, which will be developed in Sec.
V.

Booth and Gubernatis make use of a quantity referred to as the fission density
[11]. In an effort to reproduce and expand upon their work, we have chosen the
following definition for the expected fission density: considering a particle
starting a flight at position $\displaystyle x_{0}$ and traveling in direction
$\displaystyle\mu=\pm 1$, its expected fission density $\displaystyle
f(x|x_{0},\mu)$ is the expected value of the number of fission events per unit
length around $\displaystyle x$. It can be written as the product of the
probability density of flying from $\displaystyle x_{0}$ to $\displaystyle x$
and having a collision at position $\displaystyle x$, and the probability of
that collision being in the fission reaction channel. The exact form of
$\displaystyle f(x|x_{0},\mu)$ depends on the transport methodology being
employed. For NWDT, this formula can be deduced by examination of Alg. 1. The
probability density of flying from $\displaystyle x_{0}$ to $\displaystyle x$
and having a real collision is

$P_{c}(x|x_{0},\mu)=\begin{cases}q\frac{\Sigma_{t}(x)}{q\Sigma_{\text{smp}}}\Sigma_{\text{smp}}e^{-\Sigma_{\text{smp}}\absolutevalue{x-x_{0}}}&\frac{x-x_{0}}{\absolutevalue{x-x_{0}}}=\mu\\\
\\\ 0&\frac{x-x_{0}}{\absolutevalue{x-x_{0}}}\not=\mu\end{cases}.$ (47)

The
$\displaystyle\Sigma_{\text{smp}}e^{-\Sigma_{\text{smp}}\absolutevalue{x-x_{0}}}$
portion is the probability density of flying from $\displaystyle x_{0}$ to
$\displaystyle x$, and having either a real or virtual collision. The factor
$\displaystyle q$ is the probability of the collision being real, while
$\displaystyle\Sigma_{t}(x)/(q\Sigma_{\text{smp}})$ is the weight correction
factor for real collisions. This must of course be combined with the fission
probability

$P_{f}(x)=\frac{\Sigma_{f}(x)}{\Sigma_{t}(x)}.$ (48)

As $\displaystyle f(x|x_{0},\mu)=P_{f}(x)P_{c}(x|x_{0},\mu)$, for the case of
NWDT we may write

$f(x|x_{0},\mu)=\begin{cases}\Sigma_{f}(x)e^{-\Sigma_{\text{smp}}(x_{0})\absolutevalue{x-x_{0}}}&\frac{x-x_{0}}{\absolutevalue{x-x_{0}}}=\mu\\\
\\\ 0&\frac{x-x_{0}}{\absolutevalue{x-x_{0}}}\not=\mu\end{cases},$ (49)

noting that the fission density is zero for all positions which cannot be
reached by the particle during the considered flight.

Figure 7: Presented is a neutron (the colored circle), beginning a flight at
$\displaystyle x_{0}$, and flying in the $\displaystyle\mu=1$ direction
(direction of flight is indicated by the arrow). For this flight, the expected
fission density, $\displaystyle f(x|x_{0},\mu)$ has been plotted as a blue
line. The limits of the region are denoted at $\displaystyle x_{R_{-}}$ and
$\displaystyle x_{R_{+}}$. The value of $\displaystyle\beta$ displayed here is
the minimum value of $\displaystyle f(x|x_{0},\mu)$ within the region. Note
that $\displaystyle f(x|x_{0},\mu)=0$ for positions that cannot be reached by
the particle.

With the expected fission density for a particle having been defined, let us
consider a fission particle $\displaystyle p$ belonging to the fission source
and ready to start a random walk, at position $\displaystyle x_{p}$, located
in an arbitrary region $\displaystyle R$ (with bounds $\displaystyle
x_{R_{-}}$ and $\displaystyle x_{R_{+}}$). The parent of this particle (which
was a member of the previous fission generation) was previously at position
$\displaystyle x_{0}$ before flying to $\displaystyle x_{p}$, and producing
the current particle of interest. Figure 7 illustrates the expected fission
density of the parent particle, as well as the bounds of region $\displaystyle
R$, and $\displaystyle x_{0}$. The depicted value of $\displaystyle\beta$ in
the plot is the minimum value of the fission density within the region, for
the flight of the parent particle:

$\beta=\min_{x\in R}f(x|x_{0},\mu).$ (50)

The idea of Booth and Gubernatis’ method is to estimate the fission density
associated with the flight as the sum of a uniform component over the region
$\displaystyle R$ and a pointwise component concentrated at the actual fission
site. For a fission particle born in $\displaystyle R$, which was induced by a
particle beginning a flight at $\displaystyle x_{0}$ and traveling in the
direction $\displaystyle\mu$, the bottom shaded portion in Fig. 7 represents
the fraction of the fission density which is uniformly distributed in
$\displaystyle R$, while the top portion represents the non-uniform fraction
of the fission density, which depends on the position $\displaystyle x_{p}$
where the parent lands and induces a fission. This indicates that, for the
particle $\displaystyle p$, the fraction of the fission density concentrated
at $\displaystyle x_{p}$ is
$\displaystyle(f(x_{p}|x_{0},\mu)-\beta)/f(x_{p}|x_{0},\mu)$, and the fraction
uniformly distributed through region $\displaystyle R$ is
$\displaystyle\beta/f(x_{p}|x_{0},\mu)$. As such, we can set the weight of
$\displaystyle p$ to be

$w_{p}=w\frac{f(x_{p}|x_{0},\mu)-\beta}{f(x_{p}|x_{0},\mu)}$ (51)

and at the same time create a uniform weight portion $\displaystyle w_{u}$ of
the particle, representing the portion of the particle which is evenly
distributed throughout $\displaystyle R$, namely

$w_{u}=w\frac{\beta}{f(x_{p}|x_{0},\mu)}.$ (52)

It is important to note that $\displaystyle w_{p}+w_{u}=w$: the net weight in
the system has not been modified, we have simply displaced a portion of
$\displaystyle w$, distributing it uniformly through region $\displaystyle R$.

Dividing particles into a pointwise and uniform portion does not in and of
itself accomplish any weight cancellation. When there are many fission
particles within the same region, however, all of their uniform weight
portions may be combined into a single weight which represents the uniformly
distributed portion of the fission source for the region. This quantity shall
be denoted as $\displaystyle U_{R}$, and is the sum of all uniform weight
portions for all particles born within the region, each coming with a sign.
Negative particles will yield negative contributions to $\displaystyle U_{R}$,
while positive particles will yield positive contributions, leading to a
cancellation. Once all particles have contributed their uniform portion to
$\displaystyle U_{R}$, $\displaystyle n$ new particles are sampled uniformly
within the region, where

$n=\left\lceil\absolutevalue{U_{R}}\right\rceil.$ (53)

Each of these uniformly sampled particles within $\displaystyle R$ then has a
weight of $\displaystyle U_{R}/n$.

These newly sampled, uniformly distributed particles belong to the same
generation as particle $\displaystyle p$, and behave exactly like
traditionally generated particles from this point on. Under this scheme, a
particle which was initially positive will produce a positive uniform portion,
and keep a positive weight; conversely, negative particles will remain
negative and produce negative uniform portions. In the case of $\displaystyle
x_{0}$ being in the same region as $\displaystyle x_{p}$, the uniform weight
portion $\displaystyle w_{u}$ must be zero, as the minimum of the expected
fission density “behind” the starting point of the flight is zero. Therefore,
particles with $\displaystyle x_{0}$ and $\displaystyle x_{p}$ in the same
region effectively do not partake in the cancellation process. These two
properties are ensured by $\displaystyle\beta$ being the minimum value of
$\displaystyle f$ in the region. However, it is possible to relax the
requirement that $\displaystyle\beta$ should be the minimum of the expected
fission density of the region. Booth and Gubernatis show in their original
work that any value of $\displaystyle\beta$ may be used, while still producing
an unbiased result [11]. In the event of
$\displaystyle\beta>f(x_{p}|x_{0},\mu)$, the sign of the particle’s pointwise
weight portion will change.

With this regional method, the efficiency of cancellation can never be 100%.
Only the uniform weight portions contribute to cancellation. For negative
weighted fission particles, the pointwise weight portions do not cancel,
always leaving some residual negative weight. Given this information, one
might think that using a value of $\displaystyle\beta$ which is larger than
the minimum fission density might be beneficial, as it would increase the
magnitude of the weight which goes into the uniform portion. This is not
necessarily the case, however, as it could lead to the pointwise portions of
positive particles becoming negative, via Eq. (51), possibly defeating the
purpose of cancellation. Sadly, this method of cancellation can not be modeled
as a linear operator, and therefore we are unable to put it in the context of
of the framework which was put forth in Section III. This makes regional
cancellation conceptually similar to (though not the same as) the
$\displaystyle\tilde{\bm{C}}$ operator mentioned in Eq. (40). It is
inconsequential if our cancellation is linear or not; so long as it
corresponds to a choice of $\displaystyle\bm{E}$ that preserves Eq. (35), the
method will be exact and unbiased.

###  IV.1 Results of the 1D methodology

We added the exact 1D cancellation method of Booth and Gubernatis to the power
iteration problem outlined in Sec. II, performing the cancellation between
generations, and before the normalization of the system weight occurs. Thirty
evenly spaced cancellation regions were used to partition the rod, and
$\displaystyle\beta$ was always taken to be the minimum of the expected
fission density within the region. With cancellation, the simulation was
stable, and was able to complete without issue. An eigenvalue of
$\displaystyle k_{\text{eff}}=0.33577\pm 0.00005$ was obtained after 120
generations (20 inactive generations). This is in excellent agreement with the
value obtained for the system using delta tracking, which was $\displaystyle
k_{\text{eff}}=0.33573\pm 0.00005$ after the same number of generations. The
weights $\displaystyle W_{+}$, $\displaystyle W_{-}$, $\displaystyle
W_{\text{net}}$, and $\displaystyle W_{\text{tot}}$ are shown in Figure 8.

Figure 8: The positive weight ($\displaystyle W_{+}$), negative weight
($\displaystyle W_{-}$), net weight ($\displaystyle W_{\text{net}}$), and
total weight ($\displaystyle W_{\text{tot}}$) in the 1D rod system, when using
30 cancellation regions and $\displaystyle 10^{6}$ initial particles.

We see that the negative weight quickly rises to an equilibrium level, near
5000, and fluctuates then about that value for the duration of the simulation.
The positive weight must naturally increase by the same amount to keep the net
weight of the system constant. This leads to an increase in the total
transported weight of approximately $\displaystyle 10^{4}$ (the net weight
still being $\displaystyle 10^{6}$). The evolution of the number of positive
and negative particles is presented in Figure 9. The number of particles in
memory stabilizes near almost $\displaystyle 2\cdot 10^{6}$, twice as many
particles as were initially used.

Figure 9: The number of positive ($\displaystyle N_{+}$), negative
($\displaystyle N_{-}$), net ($\displaystyle N_{\text{net}}$), and total
($\displaystyle N_{\text{tot}}$) particles in the 1D rod system, when using 30
cancellation regions and $\displaystyle 10^{6}$ initial particles.

The cancellation algorithm itself has linear computational complexity with
respect to the number of fission particles at any given generation. However,
it is less obvious whether the total number of fission particles present once
the simulation has settled to equilibrium is itself a linear function of the
net weight of the system. It is evident that the total amount of weight
transported must increase until equilibrium is reached, but is the increase
linear in the net weight of the system? To examine this, we performed several
runs with different initial values of $\displaystyle W_{\text{net}}$ and we
looked at the average value of $\displaystyle W_{-}/W_{\text{net}}$, as a
function of $\displaystyle W_{\text{net}}$. The total weight can be deduced by
looking at only the behavior of the negative weight, as the increase in the
positive weight will mirror the negative weight. The result of this study is
presented in Fig. 10.

Figure 10: Negative weight fraction ($\displaystyle W_{-}/W_{\text{net}}$) as
a function of net weight ($\displaystyle W_{\text{net}}$) in the 1D rod
system, when using cancellation with 30 cancellation regions.

When too little net weight is injected in the system, a higher percentage of
negative weight will be present in the system at equilibrium. There appears to
be a critical point (in this case, near $\displaystyle W_{\text{net}}=400$)
above which the fraction of negative weight no longer decreases by adding more
net weight. This might at first seem counter-intuitive, but positive particles
are always being converted to negative particles during the random walk
(because of the NWDT algorithm), and the cancellation process is not 100%
efficient. Above this critical point, no matter how much positive weight is
added to the system at the beginning of the simulation, the equilibrium amount
of negative weight will on average be a set fraction of the net weight. In
other words, the total equilibrium weight is a linear function of the starting
weight only for sufficiently large values of the latter.

Another consideration is the effect of the number of cancellation regions on
the amount of negative weight. When there are fewer regions, each one must
become larger; this makes cancellation less efficient, as $\displaystyle\beta$
will decrease if one is taking $\displaystyle\beta$ to be the minimum of
$\displaystyle f$ over the region. Conversely, making regions too small will
result in too few particles which can partake in cancellation, also reducing
the efficiency. To examine this, we have plotted the negative weight fraction
for the number of cancellation regions in Figure 11.

Figure 11: Negative weight fraction ($\displaystyle W_{-}/W_{\text{net}}$) as
a function of the number of cancellation regions, for different net weights.

Indeed, from the figure, there is a range for which the fraction of negative
weight decreases (i.e. the efficiency of cancellation increases) by adding
more regions. The number of regions where the minimum occurs depends on the
net weight of the system. Here, this occurs near $\displaystyle
6\text{\times}{10}^{2}$ regions for $\displaystyle
W_{\text{net}}=$\displaystyle{10}^{4}$$, $\displaystyle
4\text{\times}{10}^{3}$ for $\displaystyle
W_{\text{net}}=$\displaystyle{10}^{5}$$ and $\displaystyle
2\text{\times}{10}^{4}$ for $\displaystyle
W_{\text{net}}=$\displaystyle{10}^{6}$$, suggesting that the optimal number of
regions does not quite scale linearly with the net initial weight. The three
tested net weights all had very similar behavior, and the efficiency improved
at the same rate when adding cancellation regions, until the minimum value was
reached. There is also a large range over which the negative weight fraction
is less than 1% of the net weight. For all net weights, this range starts at
approximately 20 regions, and goes up to $\displaystyle
2\text{\times}{10}^{3}$ regions for $\displaystyle
W_{\text{net}}=$\displaystyle{10}^{4}$$, or $\displaystyle
3\text{\times}{10}^{5}$ regions for $\displaystyle
W_{\text{net}}=$\displaystyle{10}^{6}$$.

This section has outlined the methodology behind Booth and Gubernatis’
regional cancellation. The potential of the method was successfully
demonstrated, having canceled enough negative weight to allow for the proper
convergence of the power iteration algorithm when using NWDT. In light of this
success, we move on to extending the method to work with more realistic
problems in higher spatial dimensions.

##  V Exact 3D, Multi-Group, Regional Cancellation Scheme

Unfortunately, it is not possible to directly implement the 1D regional
cancellation scheme of Booth and Gubernatis in higher spatial dimensions. One-
dimensional systems have the special property that lines coincide with
volumes, ensuring the lines of flight of all particles traversing the region
overlap, as well as their uniform weight portions. This fact allows us to add
the uniform weight portions of all particles in the region, effectively
leading to a cancellation of negative and positive weight. For higher
dimensions, this is not the case. The uniform portion of the particle weights
may only be distributed along their rays of flight within the cancellation
region, but the flight rays will never overlap completely. It is then no
longer possible to combine the uniform weight portions as before by simply
taking their sum.

The problem with the method of calculating the fission density in Eq. (49) is
that it is only valid along the ray of flight of the parent particle, after
the flight direction has been sampled. To extend the regional cancellation to
higher dimensions, it is required to consider the probability of the parent
particle scattering into the solid angle which is subtended by the
cancellation region. This is done by examining
$\displaystyle\hat{\bm{\Omega}}^{\prime}$, the direction of the parent
particle at position $\displaystyle\bm{r}_{0}$ before the direction is changed
by the scattering kernel. After the parent has been modified by the process of
scattering, it then has a direction $\displaystyle\hat{\bm{\Omega}}$, which
must intersect the cancellation region (given that we are attempting to
perform cancellation for a particle in the region which was induced by the
parent). The scattering cosine for the interaction which the parent underwent
at $\displaystyle\bm{r}_{0}$ is then
$\displaystyle\mu=\hat{\bm{\Omega}}^{\prime}\cdot\hat{\bm{\Omega}}$. For
neutron transport, the azimuthal direction of a scatter is almost always
isotropic, and only the scattering cosine is anisotropic (should there be any
anisotropy) [5]. The probability density of scattering in direction
$\displaystyle\hat{\bm{\Omega}}$ from direction
$\displaystyle\hat{\bm{\Omega}}^{\prime}$ is therefore

$P({\hat{\bm{\Omega}}}|{\hat{\bm{\Omega}}}^{\prime})=\frac{P(\mu)}{2\pi},$
(54)

with $\displaystyle\int_{-1}^{1}P(\mu)d\mu=1$. Modifying Eq. (49) to use the
flight kernel for 3D, and applying this factor to consider the scattering
angle, we arrive at the expected fission density

$\displaystyle\displaystyle
f(\bm{r}|\bm{r}_{0},\hat{\bm{\Omega}}^{\prime})=\frac{P(\mu)\Sigma_{f}(\bm{r})}{2\pi\absolutevalue{\bm{r}-\bm{r}_{0}}^{2}}e^{-\Sigma_{\text{smp}}(\bm{r}_{0})\absolutevalue{\bm{r}-\bm{r}_{0}}},$
(55)

where

$\displaystyle\displaystyle\mu=\frac{\bm{r}-\bm{r}_{0}}{\absolutevalue{\bm{r}-\bm{r}_{0}}}\cdot\hat{\bm{\Omega}}^{\prime}.$
(56)

Taking a closer look at Eq. (55), one may recognize it as being what Lux and
Koblinger refer to as the next-event estimator, for the fission rate at
$\displaystyle\bm{r}$ given a collision at $\displaystyle\bm{r}_{0}$ [14, Sec.
6.IV.A]. Note that Eq. (55) represents the expected fission density over the
scattering and the free flight following it; this should be contrasted with
Eq. (49), which represents the expected fission density over the following
free flight _only_.

Figure 12: Depiction of the 3D regional cancellation process. To find the
uniform weight fraction of the fission particle at $\displaystyle\bm{r}$,
$\displaystyle\beta$ must first be determined. This is done by calculating the
expected fission density at all eight corners of the cancellation region
(though only two are depicted by the dashed arrows) from the parent particles
previous position and direction ($\displaystyle\bm{r}_{0}$ and
$\displaystyle\hat{\bm{\Omega}}^{\prime}$), and taking the minimum value.

While there is only one natural shape for a cancellation region in 1D, there
are an infinite number of possible shapes in 3D which one could use. In this
work, we will consider right rectangular prisms for our cancellation region.
This is due to the fact that it is very simple to sample positions within a
right rectangular prism uniformly, which is required in the cancellation
processes. We will choose $\displaystyle\beta$ to be the minimum value of the
expected fission density within the cancellation region. If
$\displaystyle\Sigma_{f}(\bm{r})$ is spatially varying within the cancellation
region and scattering is anisotropic, it is difficult to determine the true
minimum; when $\displaystyle\Sigma_{f}(\bm{r})$ is homogeneous within the
region and scattering is isotropic, however, one may determine the true
minimum of the expected fission density for the flight from
$\displaystyle\bm{r}_{0}$ to $\displaystyle\bm{r}$ by evaluating
$\displaystyle f(\bm{r}|\bm{r}_{0},\hat{\bm{\Omega}}^{\prime})$ at all eight
corners of the prism, as the minimum must occur at one of those eight points.
Figure 12 provides a depiction of the cancellation process. Using the minimum
as the value of $\displaystyle\beta$ for the flight, it is possible to
calculate a pointwise weight portion and a uniform weight portion of the
fission particle, as was done before in Eq. (51) and Eq. (52). Once the
uniform portions of all fission particles in the region have been collected,
one may take the net uniform weight of the region, and sample new particles to
add to the fission bank uniformly within the prism.

Until now, for simplicity, we have only considered single-speed transport,
with all cross sections being energy-independent. With no energy dependence,
one does not need to consider how to sample the energy of the uniform
particles. In all realistic applications in reactor physics, however, all
cross sections are continuous functions of the neutron energy, and the
emission spectrum for the energy of fission neutrons depends on the collided
nuclide and on the incoming energy of the neutron which induced the fission
[5]. Under this assumption, it would be impossible to collect uniform portions
of all the fission particles together. Without this action, no cancellation
occurs. This fact poses a difficulty in applying weight cancellation to
continuous energy systems.

An approximation which is often used in the reactor physics community is the
multi-group treatment of the energy variable in the phase space, by using
cross sections which are piece-wise constant in energy. These energy intervals
in which the cross section is constant are referred to as groups. We denote
the $\displaystyle g$-th energy group as $\displaystyle E_{g}$, with
$\displaystyle g\in\\{1,2,..,N_{g}\\}$. By convention, the group corresponding
with the highest energy is $\displaystyle E_{1}$, while the group with the
lowest energy is $\displaystyle E_{N_{g}}$ [5]. The multi-group $\displaystyle
k$-eigenvalue Boltzmann equation reads

$\displaystyle\displaystyle\hat{\bm{\Omega}}\cdot\nabla\varphi(\bm{r},\hat{\bm{\Omega}},E_{g})+\Sigma_{t}(\bm{r},E_{g})\varphi(\bm{r},\hat{\bm{\Omega}},E_{g})=\sum_{g^{\prime}=1}^{N_{g}}\int_{4\pi}\Sigma_{s}(\bm{r},\hat{\bm{\Omega}}^{\prime}\rightarrow\hat{\bm{\Omega}},E_{g^{\prime}}\rightarrow
E_{g})\varphi(\bm{r},\hat{\bm{\Omega}}^{\prime},E_{g^{\prime}})\differential\hat{\bm{\Omega}}^{\prime}+$
$\displaystyle\displaystyle\frac{1}{4\pi
k}\sum_{g^{\prime}=1}^{N_{g}}\chi(E_{g^{\prime}}\rightarrow
E_{g})\nu(\bm{r},E_{g^{\prime}})\Sigma_{f}(\bm{r},E_{g^{\prime}})\int_{4\pi}\varphi(\bm{r},\hat{\bm{\Omega}}^{\prime},E_{g^{\prime}})\differential\hat{\bm{\Omega}}^{\prime}$
(57)

The structure of Eq. (57) is formally that of a system of particles with
$\displaystyle N_{g}$ species, coupled with each other by means of scattering
or fission.

Here, the fission spectrum $\displaystyle\chi(E_{g^{\prime}}\rightarrow
E_{g})$ is written to show a dependence of the fission particle energy on the
parent particle’s incident energy. This dependence on the incoming particle
energy is relatively weak and is often ignored in the multi-group
approximation [5]. In this case, we may simply write
$\displaystyle\chi(E_{g})$, moving it outside the sum in Eq. (57). Under this
approximation, we may collect the uniform portions of all fission neutrons in
the region together; regardless of the energy of their parent, their energy
spectra are all the same. However, one must change the definition of the
expected fission density to take energy into account. Equation (55) is then
modified to include the energy group $\displaystyle E_{g}$ of the particle
that induces fission:

$\displaystyle\displaystyle
f(\bm{r}|\bm{r}_{0},\hat{\bm{\Omega}}^{\prime},E_{g})=\frac{P(\mu,E_{g})\Sigma_{f}(\bm{r},E_{g})}{2\pi\absolutevalue{\bm{r}-\bm{r}_{0}}^{2}}$
$\displaystyle\displaystyle\times
e^{-\Sigma_{\text{smp}}(\bm{r}_{0},E_{g})\absolutevalue{\bm{r}-\bm{r}_{0}}}.$
(58)

A similar problem can occur for the direction of the fission particles as
well. In continuous energy transport, both the direction and energy of the
fission particles may be a function of the energy and direction of the
incoming neutron. While the nuclear data representations allow for this, it is
a very marginal occurrence333One exception to this is the evaluation for
[232]Th which has anisotropic distributions for prompt fission neutrons in the
laboratory frame, in both ENDF/B-VIII.0 and JEFF-3.3. in major nuclear-data
evaluation libraries [33, 34]. Therefore, we do not need to make any special
considerations for the angular distribution of fission particles in order to
apply regional cancellation to multi-group problems.

We have seen that weight cancellation takes place at fission; however, one may
wonder what is special about fission, and whether one could perform
cancellation at scattering events, instead. At fission, the energy and angular
distributions of the secondary fission particles are independent of the
properties of the particle that induces the fission event. Fission represents
the natural cancellation event because the three-dimensional distribution of
the expected fission density (Eq. (58)) encodes all the six-dimensional
distribution of the fission emission density in phase space.

##  VI Implementation and Results

Figure 13: Geometric configuration for our modified version of the C5G7
benchmark, with boundary conditions. Each color represents a material,
associated with a unique set of cross sections.

To test our exact, 3D multi-group regional cancellation, we used a modified
version of the C5G7 international benchmark, which we have depicted in Figure
13 [35]. The C5G7 represents a small 1/8th nuclear reactor core, with four
fuel assemblies and 7 energy groups, that is customarily used to assess and
compare deterministic transport codes. In its original specifications, the
fuel pins within the assemblies are cylindrical, with a radius of
$\displaystyle 0.54\text{\,}\mathrm{cm}$ [35]. The use of cylindrical fuel
pins however makes it difficult to use arbitrarily small cancellation regions
while also being able to directly calculate the minimum value of Eq. (58), to
be used as $\displaystyle\beta$. In order to facilitate cancellation, we have
modified the benchmark so that the fuel pins are square in the $\displaystyle
x$-$\displaystyle y$ plane, with side lengths of $\displaystyle
0.756\text{\,}\mathrm{cm}$. This choice allows the fuel cells to be easily cut
into an integral number of rectangular prisms to be used as cancellation
regions, ensuring that no material other than fuel is present in the region.
An important difference from the 1D rod system is that, in the C5G7 benchmark,
the cross sections are piece-wise spatially constant, being homogeneous within
a given material cell. Also, while we derived the 3D cancellation formulas for
the general case of anisotropic scattering, all scattering in the C5G7
benchmark is assumed to be isotropic, simplifying Eq. (58) to

$f(\bm{r}|\bm{r}_{0},E_{g})=\frac{\Sigma_{f}(\bm{r},E_{g})}{4\pi\absolutevalue{\bm{r}-\bm{r}_{0}}^{2}}e^{-\Sigma_{\text{smp}}(\bm{r}_{0},E_{g})\absolutevalue{\bm{r}-\bm{r}_{0}}}.$
(59)

Under this assumption, we have no need to keep track of the parent particle’s
initial direction $\displaystyle\hat{\bm{\Omega}}^{\prime}$, as there is equal
probability of scattering in any direction.

The C5G7 benchmark makes use of 3 reflective boundary conditions, as shown in
Figure 13. One caveat of the regional cancellation method is that Eq. (58) is
akin to a next-event point reaction rate estimator. Due to this, a trivial
implementation of the method can not be used with reflective boundary
conditions [28]. Use of this method with reflective boundary conditions is
only possible if special book-keeping is done to ensure a proper calculation
of the expected fission density. If a particle begins at
$\displaystyle\bm{r}_{0}$ and encounters a reflection before having a real
collision at $\displaystyle\bm{r}$, then the distances which one would obtain
by simply computing $\displaystyle\absolutevalue{\bm{r}-\bm{r}_{0}}$ will not
represent the actual distance traveled. It is also not enough to simply store
the distance traveled by the particle, as one should ideally be able to
calculate the value of $\displaystyle f$ for any given $\displaystyle\bm{r}$
in the region, in order to calculate $\displaystyle\beta$. To ensure a proper
calculation of $\displaystyle\absolutevalue{\bm{r}-\bm{r}_{0}}$, the initial
position $\displaystyle\bm{r}_{0}$ of the parent particle is stored: upon each
reflection, this position is transformed to the mirrored location on the other
side of the plane of reflection. When doing this, the point will generally lie
outside of the defined geometry for the system, but this should not pose a
problem, as we only need to be able to calculate the flight distance between
$\displaystyle\bm{r}_{0}$ and other points in the cancellation region. With
delta tracking and NWDT, $\displaystyle\Sigma_{\text{maj}}(E_{g})$ and
$\displaystyle\Sigma_{\text{smp}}(E_{g})$ do not depend on the position within
the geometry, but only on the energy group. We therefore do not need to worry
about trying to find what material this fictitious point would be located in,
or the distance to the next surface.

A problem of this type is admittedly a simple case, where NWDT is not strictly
necessary. One is able to trivially determine the majorant cross section so
that delta tracking could be used, avoiding the problem of negative weights
all together. As mentioned in Sec. II, however, in the context of multi-
physics problems which are represented with spatially continuous cross
sections, determining the majorant exactly is likely to be difficult, or
impossible. It would certainly be interesting to test this algorithm on a more
complex problem, where delta tracking would not be possible (or would be very
inefficient). It would be difficult to verify that the presented cancellation
method is working with such a problem, however, as it could not easily be
solved with existing methods. This is why we have chosen to examine this
simpler problem.

To test and evaluate cancellation and NWDT, a multi-group Monte Carlo code,
called MGMC, was written to solve $\displaystyle k$-eigenvalue power iteration
problems. MGMC supports general geometries using traditional surface-based
descriptions of volumes. Either delta tracking [20] or the variant of negative
weighted delta tracking developed by Carter, Cashwell, and Taylor [24, 22] may
be selected for transport. Scalar flux and the fission reaction rate may be
scored over a rectilinear mesh, using collision estimators444Scalar flux is
defined as
$\displaystyle\int_{4\pi}\varphi(\bm{r},\hat{\bm{\Omega}})\differential\hat{\bm{\Omega}}\text{
}$[5].. Monte Carlo estimates are saved as binary Numpy files [36] for easy
analysis and plotting with Python. Cancellation regions are defined by a
rectilinear mesh imposed on top of the problem, and can be used with both
delta tracking and negative weighted delta tracking. Shared memory parallelism
is implemented with OpenMP. Geometry, material properties, scores,
cancellation, and simulation settings are all controlled with a YAML input
file. MGMC is written in C++17, and has been made available under the
CeCILL-v2.1 license [37]. One will also find the necessary input files to
replicate our results there.

All simulations were started with $\displaystyle 10^{6}$ particles uniformly
distributed across the four fuel assemblies, all in the first energy group,
and ran for 2200 generations, with the first 200 generations being discarded
to allow for fission source convergence. Delta tracking was first run (without
weight cancellation) to obtain a reference multiplication factor for the
system, which was found to be $\displaystyle k_{\text{eff}}=1.21912\pm
0.00002$. When using the Carter, Cashwell, and Taylor variant of negative
weighted delta tracking, the majorant cross section was used for the sampling
cross section in all the energy groups except for the first; for the first
energy group, $\displaystyle 0.9\Sigma_{\text{maj}}(E_{1})$ was used as the
sampling cross section. With this choice of
$\displaystyle\Sigma_{\text{smp}}$, the total cross section is underestimated
in all of the fuel pins in the first energy group. Any virtual collisions
which occur in these regions of the phase space will therefore result in the
particle weight changing sign.

Running this simulation with NWDT, and without the use of any weight
cancellation, causes the particle populations to diverge and the simulation to
fail. The exponential increase in total weight is depicted in Figure 14. The
behavior is nearly identical to the 1D case presented in Section II.
Asymptotically, the negative weight will increase at the same rate as the
positive weight, with the difference between the two remaining constant. The
exponential increase in $\displaystyle W_{\text{tot}}$ also leads to an
exponential increase in the total number of particles in the simulation,
overwhelming computer memory, exactly as in the 1D case. This supports the
point made by Eq (29); weight cancellation will always be necessary when
attempting to perform power iteration simulations using NWDT.

Figure 14: The positive, negative, net, and total weights in the modified C5G7
benchmark, using negative weighted delta tracking and no weight cancellation.

In order to implement regional cancellation, all fuel pins were divided into
cubical cancellation regions with side lengths of $\displaystyle
0.252\text{\,}\mathrm{cm}$. This mesh was chosen by trial and error, as we
have found no practical way to know in advance for any system the minimum
number of cancellation regions required to stabilize the simulation. Cubical
cancellation regions also seemed to be more efficient than regions which had
aspect ratios much larger than one, as compact region shapes lead to a higher
minimum value of $\displaystyle f$ in the region, and therefore increases
$\displaystyle\beta$.

Our 3D multi-group exact regional cancellation algorithm was able to stabilize
the particle populations when using NWDT in conjunction with the previously
outlined parameters, and an eigenvalue of $\displaystyle
k_{\text{eff}}=1.21915\pm 0.00005$ was obtained. This is in very good
agreement with the value of $\displaystyle k_{\text{eff}}$ obtained from
traditional delta tracking with no weight cancellation, as the two estimates
differ by less than one standard deviation. The thermal (7th energy group)
scalar flux which was tallied is shown in Figure 15. This was compared with
the flux tally from the delta tracking simulation by performing a Student
$\displaystyle t$-test in each mesh cell [38]. The plot of the absolute value
of the Student $\displaystyle t$-statistic is provided in Figure 16.

Figure 15: Scalar flux of the 7th energy group in the center axial slice of
the core, obtained using negative weighted delta tracking and 3D regional
cancellation. Figure 16: Absolute value of the Student $\displaystyle
t$-statistic for the difference between the delta tracking and NWDT flux
estimates. The portion displayed is for the 7th energy group, at center axial
slice. Figure 17: Histogram of the Student $\displaystyle t$-statistic for the
difference between the flux estimates obtained with delta tracking and NWDT.
The curve plotted on top is the expected normal distribution.

It can be seen there that the two flux estimates are also in good agreement
with one another, with no apparent spatial dependence in their difference.

These two maps only show a small portion of the examined phase space. For a
more thorough comparison, a histogram of the distribution of the Student
$\displaystyle t$-statistic for the NWDT flux and the delta tracking flux is
given in Figure 17 [38]. Large portions of the flux tally mesh have average
values of zero, due to the large water reflectors. These elements were removed
before producing Figure 17. Elements which had a relative error greater than
20% (in either the NWDT or the delta tracking score) were also removed, as the
Student $\displaystyle t$-test is only applicable on normally distributed
variables; therefore, we need to ensure that sufficient statistics are
accumulated, so that the central limit theorem may apply. The 20% cut-off is
admittedly somewhat arbitrary, but it is likely that bins with such a high
relative error are not normally distributed. It can be seen that the histogram
has very good agreement with the theoretical distribution, which we assume to
be normal given the large number of degrees of freedom for the comparison.

Finally, Figure 18 shows the behaviour of the weights in the system for the
first 100 generations in the simulation. As with the 1D case, there is
initially no negative weight. This increases rapidly in the first 10-20
generations, before beginning to level out. By the time 80 generations have
passed, the weights have reached their equilibrium values. This is much longer
than the 1D case, which only required 2–3 generations before the weights
reached equilibrium. The longer time to convergence is attributed to the C5G7
benchmark having a higher dominance ratio than the 1D case.

Figure 18: The positive, negative, net, and total weights in the modified C5G7
benchmark, using negative weighted delta tracking with exact 3D regional
weight cancellation.

##  VII Conclusions

Over the course of this work, we have presented a previously undocumented
population control problem which arises when attempting to run $\displaystyle
k$-eigenvalue power iteration Monte Carlo simulations with both negative and
positive weights in negative weighted delta tracking [22]. Modeling the
transport process through a set of coupled Boltzmann transport equations for
negative and positive particles, we were able to show that the power iteration
technique applied to negative weighted delta tracking will always fail to
converge, as the sought physical fundamental eigenvalue is not the dominant
eigenvalue of the system. Instead, a fictitious eigenvalue, associated with a
system with a lower amount of absorption, is now dominant. We have formally
developed a cancellation operator, demonstrating theoretically, and with a
deterministic model, how particle weight cancellation can suppress the
fictitious eigenstate and restore convergence of power iteration to the
physical fundamental eigenvalue.

To demonstrate weight cancellation in a Monte Carlo context, we have
implemented the exact 1D region cancellation algorithm of Booth and Gubernatis
[11] in a 1D rod model. Weight cancellation did stabilize the particle
populations through the fission generations, allowing the simulation to
finish. As the method outlined by Booth and Gubernatis is only valid in 1D
single-speed problems, we have developed an exact 3D regional cancellation
method. We tested our 3D algorithm on a modified version of the C5G7 reactor-
physics benchmark, using negative weighted delta tracking. Cancellation
stabilized the particle populations, and also resulted in estimates for the
multiplication factor and flux which were in agreement with the reference
results obtained through delta tracking. Our exact 3D multi-group regional
cancellation could potentially be useful for other applications in neutron
transport as well. One such case is the use of Monte Carlo methods to obtain
higher harmonics of the $\displaystyle k$-eigenvalue equation [8]. Our
cancellation algorithm could also potentially prove useful in transport
problems involving complex particle weights, such as the search for critical
buckling or the solution of the neutron noise equations [39, 7].

Several questions remain to be settled in regard to the methodology of 3D
cancellation. First, in this work we have required that a cancellation region
consist of only one material. Being able to have multiple materials in a
cancellation region would make it easier to perform cancellation on different
geometric forms (cylindrical pins). Second, extension of this algorithm to
continuous energy will require the application of new techniques, because in
continuous-energy transport the fission spectrum depends on the incoming
neutron energy. Finally, the efficient choice of $\displaystyle\beta$ should
be evaluated. While we have chosen to use the minimum value of the expected
fission density in the region for $\displaystyle\beta$, Booth and Gubernatis
made clear that any value of $\displaystyle\beta$ results in an unbiased
estimate of the fission density. However, the variance of the weights of the
fission particles clearly depends on $\displaystyle\beta$. It is desirable to
characterize which choices of $\displaystyle\beta$ result in more efficient
cancellation. All of these issues will require further investigation to
improve the algorithm, and to probe its possible extension to continuous-
energy problems.

*

##  Appendix A Derivation of coupled Boltzmann transport equations

We provide here the derivation of the coupled system of Boltzmann equations,
where the populations of positive and negative particles are treated
separately and are assumed to be governed by the rules of negative weighted
delta tracking (Alg. 1) [22]. The quantities that we wish to describe are the
angular fluxes of positive and negative particles. These quantities are not
physical observables, but they are useful to characterize the behaviour of the
random walk. Consider any one of the usual Monte Carlo estimators for the
(physical) angular flux; the unbiasedness condition requires any such
estimator to be proportional to the particle weight $\displaystyle w$.
Replacing the particle weight $\displaystyle w$ in the estimator with

$w_{+}=\begin{cases}w&w>0\\\ 0&w\leq 0\end{cases}$ (60)

results in a modified estimator; we define the angular flux of positive
particles $\displaystyle\varphi_{+}$ to be the expected value of the modified
estimator. Likewise, replacing the particle weight with

$w_{-}=\begin{cases}0&w>0\\\ -w&w\leq 0\end{cases}$ (61)

results in a modified estimator, whose expected value is defined to be the
angular flux of negative particles, $\displaystyle\varphi_{-}$.

A few properties are worth stressing. First, if all particle weights are
positive, then $\displaystyle\varphi_{+}=\varphi$ (the physical angular flux)
and $\displaystyle\varphi_{-}=0$. Second, since $\displaystyle w_{+}-w_{-}=w$,
then $\displaystyle\varphi_{+}-\varphi_{-}=\varphi$. Finally, since
$\displaystyle w_{+}\geq 0$ and $\displaystyle w_{-}\geq 0$, then both
$\displaystyle\varphi_{+}$ and $\displaystyle\varphi_{-}$ are always non-
negative.

With these postulations, we may commence our derivation of the angular flux of
positive particles. We shall obey the rules for negative weighted delta
tracking [22], presented in Alg. 1. The Boltzmann equation for neutron
transport is canonically written in the form of a balance equation, where
losses = gains. These losses and gains refer to the change in the neutron flux
at the phase space point $\displaystyle(\bm{r},\hat{\bm{\Omega}})$. This can
be referenced in Eq. (1), where the left hand side (LHS) of the equation
represents losses, while the right hand side (RHS) represents gains. More
detail as to the reasoning and derivation behind particular terms can be found
in Nuclear Reactor Theory, by Bell and Glasstone [5].

Starting with the LHS, we first consider losses due to particle streaming,
$\displaystyle L_{S}$. This term is identical to that found in Eq. (1), only
replacing the physical angular flux with the angular flux of positive
particles:

$\displaystyle\displaystyle
L_{S}=\hat{\bm{\Omega}}\cdot\nabla\varphi_{+}(\bm{r},\hat{\bm{\Omega}}).$ (62)

Next, we consider losses due to collisions, $\displaystyle L_{C}$. In the case
of NWDT, it is possible to be removed from the phase space point by having a
collision at $\displaystyle\bm{r}$ (real or virtual), which occurs with a
cross section $\displaystyle\Sigma_{\text{smp}}$. We therefore consider the
losses due to collisions with the term

$\displaystyle\displaystyle
L_{C}=\Sigma_{\text{smp}}(\bm{r})\varphi_{+}(\bm{r},\hat{\bm{\Omega}}).$ (63)

In the event of a real collision, the particle direction will change, so there
is indeed a removal from the phase space point. If a virtual collision occurs
and the positive particle changes sign to be a negative particle, this is also
a loss. If a virtual collision does not result in a positive particle becoming
negative, then it is not actually a loss. This case will be treated on the RHS
by adding a gain term for virtual collisions which do not result in a sign
flip. The only possible sources for loss is then $\displaystyle L_{S}$ and
$\displaystyle L_{C}$, and we then have a LHS of

$L_{S}+L_{C}=\hat{\bm{\Omega}}\cdot\nabla\varphi_{+}(\bm{r},\hat{\bm{\Omega}})+\Sigma_{\text{smp}}(\bm{r})\varphi_{+}(\bm{r},\hat{\bm{\Omega}}).$
(64)

Next, we may consider gains at the phase space point
$\displaystyle(\bm{r},\hat{\bm{\Omega}})$. Gains may only come from particles
which have had either a real or virtual collision, allowing them to enter the
phase space point in question. First, we shall consider the gains from real
collisions, $\displaystyle G_{R}$. A real collision occurs with probability
$\displaystyle q(\bm{r})$, and when this happens, a weight modification also
occurs, which is a multiplication of the factor
$\displaystyle\Sigma_{t}(\bm{r})/(q(\bm{r})\Sigma_{\text{smp}}(\bm{r}))$:

$\displaystyle\displaystyle
G_{R}=q(\bm{r})\frac{\Sigma_{t}(\bm{r})}{q(\bm{r})\Sigma_{\text{smp}}(\bm{r})}R,$
(65)

here $\displaystyle R$ is the sum of the gains of all real collision channels.
We have two possible reaction channels which contribute to gains at
$\displaystyle(\bm{r},\hat{\bm{\Omega}})$: particles having a collision at
$\displaystyle\bm{r}$ and scattering from direction
$\displaystyle\hat{\bm{\Omega}}^{\prime}$ into direction
$\displaystyle\hat{\bm{\Omega}}$ ($\displaystyle R_{S}$), and particles having
a collision at $\displaystyle\bm{r}$, inducing a fission particle which is
born traveling in direction $\displaystyle\hat{\bm{\Omega}}$ ($\displaystyle
R_{F}$). Treating the scattering term, we must sum over all possible incoming
directions, resulting in

$R_{S}=\frac{\Sigma_{\text{smp}}(\bm{r})}{\Sigma_{t}(\bm{r})}\int_{4\pi}\Sigma_{s}(\bm{r},\hat{\bm{\Omega}}^{\prime}\rightarrow\hat{\bm{\Omega}})\varphi_{+}(\bm{r},\hat{\bm{\Omega}}^{\prime})\differential\hat{\bm{\Omega}}^{\prime}.$
(66)

We must multiply by the ratio $\displaystyle\Sigma_{\text{smp}}/\Sigma_{t}$ to
account for the fact that we are conditioning on the collision being real. We
now perform a similar operation for the fission channel, with the assumption
that the distribution of the direction of fission particles is isotropic:

$\displaystyle\displaystyle
R_{F}=\frac{\Sigma_{\text{smp}}(\bm{r})}{\Sigma_{t}(\bm{r})}\frac{\nu(\bm{r})}{4\pi
k}\int_{4\pi}\Sigma_{f}(\bm{r})\varphi_{+}(\bm{r},\hat{\bm{\Omega}}^{\prime})\differential\hat{\bm{\Omega}}^{\prime}.$
(67)

These being the only two real collision channels ($\displaystyle
R=R_{S}+R_{F}$), we may combine the definitions of $\displaystyle R_{S}$,
$\displaystyle R_{F}$, and $\displaystyle G_{R}$ to obtain the result
presented in Eq. (68), where we have employed our previous definitions for
$\displaystyle\mathcal{S}$ and $\displaystyle\mathcal{F}$ from Eq. (10) and
Eq. (11).

$\displaystyle\displaystyle G_{R}=$ $\displaystyle\displaystyle
q(\bm{r})\frac{\Sigma_{t}(\bm{r})}{q(\bm{r})\Sigma_{\text{smp}}(\bm{r})}\bigg{[}\frac{\Sigma_{\text{smp}}(\bm{r})}{\Sigma_{t}(\bm{r})}\int_{4\pi}\Sigma_{s}(\bm{r},\hat{\bm{\Omega}}^{\prime}\rightarrow\hat{\bm{\Omega}})\varphi_{+}(\bm{r},\hat{\bm{\Omega}}^{\prime})\differential\hat{\bm{\Omega}}^{\prime}+\frac{\Sigma_{\text{smp}}(\bm{r})}{\Sigma_{t}(\bm{r})}\frac{\nu(\bm{r})}{4\pi
k}\int_{4\pi}\Sigma_{f}(\bm{r})\varphi_{+}(\bm{r},\hat{\bm{\Omega}}^{\prime})\differential\hat{\bm{\Omega}}^{\prime}\bigg{]}$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\int_{4\pi}\Sigma_{s}(\bm{r},\hat{\bm{\Omega}}^{\prime}\rightarrow\hat{\bm{\Omega}})\varphi_{+}(\bm{r},\hat{\bm{\Omega}}^{\prime})\differential\hat{\bm{\Omega}}^{\prime}+\frac{\nu(\bm{r})}{4\pi
k}\int_{4\pi}\Sigma_{f}(\bm{r})\varphi_{+}(\bm{r},\hat{\bm{\Omega}}^{\prime})\differential\hat{\bm{\Omega}}^{\prime}$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathcal{S}\varphi_{+}+\frac{1}{k}\mathcal{F}\varphi_{+}$
(68)

All that remains are gains from virtual collisions $\displaystyle G_{V}$. The
probability of a virtual collision is $\displaystyle 1-q(\bm{r})$, and is
accompanied by a weight modification of

$\displaystyle\displaystyle\frac{\absolutevalue{1-\frac{\Sigma_{t}(\bm{r})}{\Sigma_{\text{smp}}(\bm{r})}}}{1-q(\bm{r})}.$
(69)

We take the absolute value here, as our particle weights may never become
negative. Instead, a change in sign is modeled by the transfer of a particle
from the positive population (corresponding to $\displaystyle\varphi_{+}$) to
the negative population (corresponding to $\displaystyle\varphi_{-}$). Our
virtual gains are

$\displaystyle\displaystyle
G_{V}=(1-q(\bm{r}))\frac{\absolutevalue{1-\frac{\Sigma_{t}(\bm{r})}{\Sigma_{\text{smp}}(\bm{r})}}}{1-q(\bm{r})}V,$
(70)

$\displaystyle V$ being the sum of gains due to virtual collisions. The first
channel for virtual collisions is the previously mentioned case of a positive
particle having a virtual collision, and remaining positive ($\displaystyle
V_{+}$). This only occurs when
$\displaystyle\Sigma_{\text{smp}}(\bm{r})\geq\Sigma_{t}(\bm{r})$, and can be
modeled with the Heaviside function

$\displaystyle\displaystyle\Theta(x)=\begin{cases}1&x\geq 0\\\
0&x<0\end{cases}.$ (71)

Our gains from virtual collisions of positive particles is then

$\displaystyle\displaystyle
V_{+}=\Sigma_{\text{smp}}(\bm{r})\Theta\big{(}\Sigma_{\text{smp}}(\bm{r})-\Sigma_{t}(\bm{r})\big{)}\varphi_{+}(\bm{r},\hat{\bm{\Omega}}).$
(72)

In the event that
$\displaystyle\Sigma_{\text{smp}}(\bm{r})<\Sigma_{t}(\bm{r})$, negative
particles will flip sign, joining the positive particles ($\displaystyle
V_{-}$). This is modeled in a similar manner, simply flipping the argument in
the Heaviside function, and replacing the positive angular flux with the
negative angular flux:

$\displaystyle\displaystyle
V_{-}=\Sigma_{\text{smp}}(\bm{r})\Theta\big{(}\Sigma_{t}(\bm{r})-\Sigma_{\text{smp}}(\bm{r})\big{)}\varphi_{-}(\bm{r},\hat{\bm{\Omega}}).$
(73)

With these terms defined, we may combine our definitions of $\displaystyle
G_{V}$ and $\displaystyle V=V_{+}+V_{-}$ to obtain the gains from virtual
collisions, presented in Eq. (74), having used the function from Eq. (12).

All loss and gain terms have now been defined, leaving us with the final form
of the equation, describing the angular flux of the positive particles in Eq.
(75). By symmetry, the equation for negative particles must have the exact
same form, and is given in Eq. (76). Regardless of whether a particle is
negative or positive, it is transported in the same manner, and by the same
rules. These two equations are of course the form presented in Eq. (9), and
demonstrate the coupling relationship between positive and negative particles.

$\displaystyle\displaystyle
G_{V}=(1-q(\bm{r}))\frac{\absolutevalue{1-\frac{\Sigma_{t}(\bm{r})}{\Sigma_{\text{smp}}(\bm{r})}}}{1-q(\bm{r})}\Sigma_{\text{smp}}(\bm{r})\bigg{[}\Theta\big{(}\Sigma_{\text{smp}}(\bm{r})-\Sigma_{t}(\bm{r})\big{)}\varphi_{+}(\bm{r},\hat{\bm{\Omega}})+\Theta\big{(}\Sigma_{t}(\bm{r})-\Sigma_{\text{smp}}(\bm{r})\big{)}\varphi_{-}(\bm{r},\hat{\bm{\Omega}})\bigg{]}$
$\displaystyle\displaystyle=\Theta\big{(}\Sigma_{\text{smp}}(\bm{r})-\Sigma_{t}(\bm{r})\big{)}\big{[}\Sigma_{\text{smp}}(\bm{r})-\Sigma_{t}(\bm{r})\big{]}\varphi_{+}(\bm{r},\hat{\bm{\Omega}})+\Theta\big{(}\Sigma_{t}(\bm{r})-\Sigma_{\text{smp}}(\bm{r})\big{)}\big{[}\Sigma_{t}(\bm{r})-\Sigma_{\text{smp}}(\bm{r})\big{]}\varphi_{-}(\bm{r},\hat{\bm{\Omega}})$
$\displaystyle\displaystyle=\Delta\big{(}\Sigma_{\text{smp}}(\bm{r})-\Sigma_{t}(\bm{r})\big{)}\varphi_{+}(\bm{r},\hat{\bm{\Omega}})+\Delta\big{(}\Sigma_{t}(\bm{r})-\Sigma_{\text{smp}}(\bm{r})\big{)}\varphi_{-}(\bm{r},\hat{\bm{\Omega}})$
(74)
$\displaystyle\displaystyle\hat{\bm{\Omega}}\cdot\nabla\varphi_{+}+\Sigma_{\text{smp}}\varphi_{+}=\mathcal{S}\varphi_{+}+\frac{1}{k}\mathcal{F}\varphi_{+}+\Delta\big{(}\Sigma_{\text{smp}}-\Sigma_{t}\big{)}\varphi_{+}+\Delta\big{(}\Sigma_{t}-\Sigma_{\text{smp}}\big{)}\varphi_{-}$
(75)
$\displaystyle\displaystyle\hat{\bm{\Omega}}\cdot\nabla\varphi_{-}+\Sigma_{\text{smp}}\varphi_{-}=\mathcal{S}\varphi_{-}+\frac{1}{k}\mathcal{F}\varphi_{-}+\Delta\big{(}\Sigma_{\text{smp}}-\Sigma_{t}\big{)}\varphi_{-}+\Delta\big{(}\Sigma_{t}-\Sigma_{\text{smp}}\big{)}\varphi_{+}\text{.}$
(76)

## References

  * Berg [1993] H. C. Berg,  _Random Walks in Biology_ (Princeton University Press, 1993).
  * Bartlett [1960] M. Bartlett,  _Stochastic Population Models in Ecology and Epidemiology_, Methuen’s monographs on applied probability and statistics (Methuen, 1960).
  * Williams [1974] M. M. R. Williams,  _Random Processes in Nuclear Reactors_ , 1st ed. (Pergamon Press Oxford, New York, 1974).
  * Modest [2013] M. F. Modest,  _Radiative Heat Transfer_ (Academic press, 2013).
  * Bell and Glasstone [1970] G. I. Bell and S. Glasstone,  _Nuclear Reactor Theory_ (Van Nostrand Reinhold Company, 1970).
  * Arnow  _et al._ [1982] D. M. Arnow, M. H. Kalos, M. A. Lee, and K. E. Schmidt, J. Chem. Phys. 77, 5562 (1982).
  * Yamamoto [2013] T. Yamamoto, Ann. Nucl. Energy 58, 72–79 (2013).
  * Booth [2003] T. E. Booth, Nucl. Sci. Eng. 143, 291 (2003).
  * Szirmay‐Kalos  _et al._ [2017] L. Szirmay‐Kalos, I. Georgiev, M. Magdics, B. Molnár, and D. Légrády, Computer Graphics Forum 36, 9–19 (2017).
  * Booth and Gubernatis [2009] T. E. Booth and J. E. Gubernatis, Phys. Rev. E 80, 046704 (2009).
  * Booth and Gubernatis [2010] T. E. Booth and J. E. Gubernatis, Nucl. Sci. Eng. 165, 283 (2010).
  * Spencer  _et al._ [2012] J. S. Spencer, N. S. Blunt, and W. M. Foulkes, J. Chem. Phys. 136, 054110 (2012).
  * Holmes [2016] M. H. Holmes,  _Introduction to Scientific Computing and Data Analysis_ , Texts in Computational Science and Engineering, Vol. 13 (Springer, 2016).
  * Lux and Koblinger [1991] I. Lux and L. Koblinger,  _Monte Carlo Particle Transport Methods: Neutron and Photon Calculations_ (CRC Press, 1991).
  * Nowak  _et al._ [2016] M. Nowak, J. Miao, E. Dumonteil, B. Forget, A. Onillon, K. S. Smith, and A. Zoia, Ann. Nucl. Energy 94, 856–868 (2016).
  * Zoia  _et al._ [2011] A. Zoia, E. Dumonteil, and A. Mazzolo, Phys. Rev. E 83, 041137 (2011).
  * Demazière [2013] C. Demazière, Int. J. Nucl. Energy Sci. Technol. 7, 10.1504/IJNEST.2013.054368 (2013).
  * Belanger  _et al._ [2020] H. Belanger, D. Mancusi, and A. Zoia, Eur. Phys. J. Plus 135, 877 (2020).
  * Demazière [2019] C. Demazière,  _Modelling of Nuclear Reactor Multi-physics: From Local Balance Equations to Macroscopic Models in Neutronics and Thermal-hydraulics_ (Academic Press, 2019).
  * Woodcock  _et al._ [1965] E. Woodcock, T. Murphy, P. Hemmings, and S. 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 [2010] J. Leppänen, Ann. Nucl. Energy 37, 715 (2010).
  * Legrady  _et al._ [2017] D. Legrady, B. Molnar, M. Klausz, and T. Major, Ann. Nucl. Energy 102, 116 (2017).
  * Molnar  _et al._ [2018] B. Molnar, G. Tolnai, and D. Legrady, Nucl. Sci. Eng. 190, 1 (2018).
  * Carter  _et al._ [1972] L. L. Carter, E. D. Cashwell, and W. M. Taylor, Nucl. Sci. Eng. 48, 403 (1972).
  * Jacoboni and Reggiani [1983] C. Jacoboni and L. Reggiani, Rev. Mod. Phys. 55, 645–705 (1983).
  * Lemaire  _et al._ [2018] V. Lemaire, M. Thieullen, and N. Thomas, J. Sci. Comput. 75, 1776 (2018).
  * Galtier  _et al._ [2013] M. Galtier, S. Blanco, C. Caliot, C. Coustet, J. Dauchet, M. E. Hafi, V. Eymet, R. Fournier, J. Gautrais, A. Khuong, and et al., J. Quant. Spectrosc. Radiat. Transfer 125, 57–68 (2013).
  * X-5 Monte Carlo Team [2003] X-5 Monte Carlo Team,  _MCNP — A General Monte Carlo N-Particle Transport Code, Version 5 Volume I: Overview and Theory_ (2003), LA-UR-03-1987.
  * 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 (2015).
  * Booth [1996] T. E. Booth,  _A weight (charge) conserving importance-weighted comb for Monte Carlo_, Tech. Rep. (Los Alamos National Laboratory, 1996) LA-UR 96-0051.
  * Rouchon  _et al._ [2017] A. Rouchon, A. Zoia, and R. Sanchez, Ann. Nucl. Energy 102, 465–475 (2017).
  * Demmel [1989] J. Demmel, in  _IEEE Control Systems Society Workshop on Computer-Aided Control System Design_ (1989) p. 1–7.
  * Brown  _et al._ [2018] D. Brown, M. Chadwick, R. Capote, A. Kahler, A. Trkov, M. Herman, A. Sonzogni, Y. Danon, A. Carlson, M. Dunn, and et al., Nucl. Data Sheets 148, 1 (2018).
  * Plompen  _et al._ [2020] A. J. M. Plompen, O. Cabellos, C. D. S. Jean, M. Fleming, A. Algora, M. Angelone, P. Archier, E. Bauge, O. Bersillon, A. Blokhin, and et al., Eur. Phys. J. A 56, 181 (2020).
  * 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. (NEA/NSC, 2001).
  * Harris  _et al._ [2020] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, and et al., Nature 585, 357–362 (2020).
  * Belanger [2021] H. Belanger, 10.5281/zenodo.4636567 (2021), https://github.com/HunterBelanger/mgmc.
  * Rice [1995] J. Rice,  _Mathematical Statistics and Data Analysis_ (Duxbury Press, Belmont, California, 1995).
  * Yamamoto [2012] T. Yamamoto, Ann. Nucl. Energy 47, 14–20 (2012).