# Variance Reduction and Noise Source Sampling Techniques for Monte Carlo
Simulations of Neutron Noise Induced by Mechanical Vibrations

Hunter Belanger<EMAIL_ADDRESS>Davide Mancusi<EMAIL_ADDRESS>Amélie Rouchon<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

Neutron noise in nuclear power reactors refers to the small fluctuations
around the average neutron flux at steady state resulting from time-dependent
perturbations inside the core. The neutron noise equations in the frequency
domain can be solved using Monte Carlo simulation codes, which are capable of
obtaining reference solutions involving almost no approximations, but are
hindered by severe issues affecting the statistical convergence: the
simultaneous presence of positive and negative particles, which is required by
the nature of the complex noise equations, leads to catastrophically large
variance in the tallies. In this work, we consider the important case of
neutron noise problems induced by mechanical vibrations. First, we derive a
new exact sampling strategy for the noise source. Then, building upon our
previous findings in other contexts, we show that weight cancellation methods
can be highly beneficial in dealing with the presence of negative weights,
enabling extremely large gains in the figure of merit. We successfully
demonstrate our results on a benchmark configuration consisting of a fuel
assembly with a vibrating pin and we discuss possible pathways for further
improvements.

##  I Introduction

Neutron noise in nuclear power reactors is defined as the small fluctuations
occurring around the average neutron flux at steady state due to perturbations
typically induced by fluid-structure interactions or moderator density
fluctuations inside the core [1, 2, 3]. These perturbations may affect
different reactor components, at the scale of single or multiple fuel rods,
entire assemblies, or the full core vessel, depending on the specific physical
origin of the disturbances. Generally speaking, neutron noise is an unwanted
phenomenon, which in some extreme cases might lead to a significant decrease
in reactor power in order to ensure safe operation, or even to the shut-down
of the reactor [4]. Notwithstanding, such disturbances carry an information
content that can be usefully exploited to improve reactor diagnostic through
the application of inverse problem techniques, for instance by locating
anomalous control rod or fuel assembly vibrations, or monitoring moderator
speed and void fraction [5, 6, 7, 8, 9, 10].

The neutron noise equations are established by considering the effect of
temperature and density changes and/or material displacements on the operators
occurring in the Boltzmann equation and in the precursor equations: the noise
is considered as a small time-dependent perturbation with respect to the
stationary flux, and the time-dependent operators are similarly decomposed
into a stationary part and a residual time-dependent perturbation [3]. Then,
assuming that the products of perturbed quantities can be neglected (the so-
called ‘orthodox’ linearization [11]), a Fourier transform finally leads to a
complex-valued equation for the neutron noise in the frequency domain [3].

A number of analytical and semi-analytical results concerning the solutions of
the noise equations have been obtained using diffusion theory, simplified few-
group transport equations, and Green’s function methods [12, 13, 14, 15, 16,
17]: although these approaches shed light on the behaviour of the neutron
noise field in a variety of physical configurations, their application to
real-world systems is somewhat limited. Recently, a renewed interest in
neutron noise stimulated by industrial problems related to core diagnostics
has fostered the development of novel numerical codes capable of solving the
noise equations in the time or frequency domain using state-of-the-art methods
[18, 19, 20, 21, 22, 23, 24]. In this respect, a remarkable contribution has
been provided by the EU H2020 CORTEX project (2017-2021) [25], during which
several deterministic and Monte Carlo solvers for the noise equation have been
conceived [26, 27, 28] and validated against benchmark problems [29, 30, 31]
and experimental data stemming from dedicated measurement campaigns in
research reactors [32, 33].

Since the number of available experimental datasets for the qualification of
faster, albeit approximate, deterministic solvers for the noise equations is
still rather small, one would like to rely on Monte Carlo simulation as a
gold-standard tool to compute reference solutions involving almost no
approximations, similarly to what is customarily done for the regular
Boltzmann equation in reactor physics and radiation shielding problems [30,
31]. Monte Carlo sampling methods devoted to the noise equations have been
devised and successfully tested in several applications, but numerical
investigations have shown that their statistical convergence may be extremely
poor in some cases: this is basically due to the fact that solving the noise
equations requires sampling particles carrying complex statistical weights
with positive and negative real and imaginary components [20, 21]. When
positive and negative contributions are tallied in order to estimate the
neutron noise, the variance might become overwhelmingly large and a huge
number of sampled particles is mandatory to achieve a reasonable statistical
uncertainty, which makes such simulations unwieldy even for relatively simple
configurations [30, 31].

The possibility of using Monte Carlo simulations as a routine tool for
industrial calculations in the field of neutron noise therefore crucially
depends on the availability of ad hoc variance reduction and population
control techniques capable of handling the simultaneous presence of positive
and negative particles. Recent insights have shown that weight cancellation
methods might be the key to successfully dealing with such issues [20, 34,
35]. Preliminary investigations involving cross section oscillation problems
in a fuel assembly suggest that the introduction of weight cancellation
actually leads to astounding improvements in the figure of merit of the noise
simulations by Monte Carlo codes, at the expense of extensive restructuring of
the sampling routines [36]. The benefits induced by the vastly reduced
statistical uncertainty establish Monte Carlo simulations as a practical tool
to obtain reference solutions and outweighs the efforts involved in code
rewriting. In this work, we extend and considerably generalize our previous
findings to the case of neutron noise induced by mechanical vibrations, which
poses distinct challenges in both sampling the noise source particles and
handling the variance of the resulting noise field.

This work is organized as follows. In Sec. II we recall the derivation of the
neutron noise equations in the frequency domain, and we show that they can be
interpreted as a complex-valued fixed-source transport problem, starting from
a complex-valued source. Subsequently, in Sec. III we review existing sampling
methods for the noise source, and introduce a new exact sampling strategy for
the case of mechanical vibrations. After briefly recalling in Sec. IV the
Monte Carlo methods for solving the transport of noise particles, in Sec. V we
discuss the application of weight cancellation to reduce the variance of noise
problems. The effectiveness of these methods is illustrated in Sec. VI on a
benchmark configuration consisting of a fuel assembly where a pin is subject
to periodic mechanical vibrations. Conclusions are finally drawn in Sec. VII.

##  II The Neutron Noise Equations

We start by considering a multiplying system which is initially critical. To
allow for possible model bias due to technological uncertainties and/or
nuclear data, we assume that the fission production is normalized by the
fundamental eigenvalue $k_{\text{eff}}$, which enforces the stationary state.
Now, if a perturbation is introduced into the system, so that the macroscopic
quantities such as cross sections, yields, and scattering laws become time-
dependent, then the resulting neutron flux
$\varphi(\bm{r},\bm{\hat{\Omega}},E,t)$ will also be time-dependent, and will
be therefore described by the time-dependent Boltzmann equation, coupled to
the evolution equations for the delayed neutron precursors. After injecting
the precursor equations into the Boltzmann equation, we concisely denote the
time-dependent transport equation as
$\mathcal{B}_{k}(t)\varphi(\bm{r},\bm{\hat{\Omega}},E,t)=0$, where

$\mathcal{B}_{k}(t)=\frac{1}{v}\frac{\partial}{\partial
t}+\bm{\hat{\Omega}}\cdot\gradient+\Sigma_{t}(\bm{r},E,t)-\iint\nu_{s}(\bm{r},E^{\prime},t)\Sigma_{s}(\bm{r},E^{\prime},t)f_{s}(E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}},t)\differential
E^{\prime}\differential\bm{\hat{\Omega}}^{\prime}-\\\
\frac{1}{k_{\text{eff}}}\iint\nu_{f,p}(\bm{r},E^{\prime},t)\Sigma_{f}(\bm{r},E^{\prime},t)f_{f,p}(E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}},t)\differential
E^{\prime}\differential\bm{\hat{\Omega}}^{\prime}-\\\
\frac{1}{k_{\text{eff}}}\sum_{j}\iiint\lambda_{j}e^{-\lambda_{j}(t-t^{\prime})}\nu_{f,d}^{j}(\bm{r},E^{\prime},t^{\prime})\Sigma_{f}(\bm{r},E^{\prime},t^{\prime})f_{f,d}^{j}(E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}},t^{\prime})\differential
E^{\prime}\differential\bm{\hat{\Omega}}^{\prime}\differential
t^{\prime}\text{.}$ (1)

Due to the effect of the perturbation, it is assumed that all of the time-
dependent terms in Eq. (1) can be be decomposed into a stationary part and an
additional small, periodic, time-dependent contribution111In this work, we do
not consider the possibility of perturbations in the precursor decay
constants.. For example, the total macroscopic cross section is decomposed as

$\Sigma_{t}(\bm{r},E,t)=\Sigma_{t}(\bm{r},E)+\delta\Sigma_{t}(\bm{r},E,t)\text{.}$
(2)

This can similarly be done for reaction yields, scattering laws and fission
spectra. As it would become slightly cumbersome to write each integral
explicitly, we use the following shorthand for the decomposition of a generic
term $\alpha$:

$\nu_{\alpha}(\bm{r},E,t)\Sigma_{\alpha}(\bm{r},E,t)f_{\alpha}(E\to
E^{\prime},\bm{\hat{\Omega}}\to\bm{\hat{\Omega}}^{\prime},t)=\\\
\nu_{\alpha}(\bm{r},E)\Sigma_{\alpha}(\bm{r},E)f_{\alpha}(E\to
E^{\prime},\bm{\hat{\Omega}}\to\bm{\hat{\Omega}}^{\prime})\\\
+\delta[\nu_{\alpha}\Sigma_{\alpha}f_{\alpha}](\bm{r},E\to
E^{\prime},\bm{\hat{\Omega}}\to\bm{\hat{\Omega}}^{\prime},t)\text{,}$ (3)

where all of the time dependence has been placed in the single
$\delta[\nu_{\alpha}\Sigma_{\alpha}f_{\alpha}]$ term. A thorough discussion on
this decomposition will be provided in Sec. III. In this formalism, we may
break the operator $\mathcal{B}_{k}(t)$ into two different components

$\mathcal{B}_{k}(t)=\mathcal{B}(t)+\delta\mathcal{B}(t)\text{,}$ (4)

with the operator

$\mathcal{B}(t)=\frac{1}{v}\frac{\partial}{\partial
t}+\bm{\hat{\Omega}}\cdot\gradient+\Sigma_{t}(\bm{r},E)-\iint\nu_{s}(\bm{r},E^{\prime})\Sigma_{s}(\bm{r},E^{\prime})f_{s}(E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}})\differential
E^{\prime}\differential\bm{\hat{\Omega}}^{\prime}-\\\
\frac{1}{k_{\text{eff}}}\iint\nu_{f,p}(\bm{r},E^{\prime})\Sigma_{f}(\bm{r},E^{\prime})f_{f,p}(E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}})\differential
E^{\prime}\differential\bm{\hat{\Omega}}^{\prime}-\\\
\frac{1}{k_{\text{eff}}}\sum_{j}\iiint\lambda_{j}e^{-\lambda_{j}(t-t^{\prime})}\nu_{f,d}^{j}(\bm{r},E^{\prime})\Sigma_{f}(\bm{r},E^{\prime})f_{f,d}^{j}(E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}})\differential
E^{\prime}\differential\bm{\hat{\Omega}}^{\prime}\differential
t^{\prime}\text{,}$ (5)

and the perturbation operator

$\delta\mathcal{B}(t)=\delta\Sigma_{t}(\bm{r},E,t)-\iint\delta[\nu_{s}\Sigma_{s}f_{s}](\bm{r},E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}},t)\differential
E^{\prime}\differential\bm{\hat{\Omega}}^{\prime}-\\\
\frac{1}{k_{\text{eff}}}\iint\delta[\nu_{f,p}\Sigma_{f,p}f_{f,p}](\bm{r},E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}},t)\differential
E^{\prime}\differential\bm{\hat{\Omega}}^{\prime}-\\\
\frac{1}{k_{\text{eff}}}\sum_{j}\iiint\lambda_{j}e^{-\lambda_{j}(t-t^{\prime})}\delta[\nu_{f,d}^{j}\Sigma_{f,d}f_{f,d}^{j}](\bm{r},E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}},t^{\prime})\differential
E^{\prime}\differential\bm{\hat{\Omega}}^{\prime}\differential
t^{\prime}\text{.}$ (6)

Consequently, we postulate that the time-dependent flux which solves Eq. (1)
can be similarly decomposed as

$\varphi(\bm{r},\bm{\hat{\Omega}},E,t)=\varphi_{c}(\bm{r},\bm{\hat{\Omega}},E)+\delta\varphi(\bm{r},\bm{\hat{\Omega}},E,t)\text{,}$
(7)

where $\varphi_{c}$ is the fundamental eigenmode of the Boltzmann equation
associated to the fundamental eigenvalue $k_{\text{eff}}$, before the time-
dependent disturbance was introduced, and the flux perturbation
$\delta\varphi$ is the neutron noise. From these definitions, we see that

$\displaystyle\mathcal{B}_{k}$
$\displaystyle(t)\varphi(\bm{r},\bm{\hat{\Omega}},E,t)=0$
$\displaystyle=\bigg{[}\mathcal{B}(t)+\delta\mathcal{B}(t)\bigg{]}\bigg{[}\varphi_{c}(\bm{r},\bm{\hat{\Omega}},E)+\delta\varphi(\bm{r},\bm{\hat{\Omega}},E,t)\bigg{]}$
$\displaystyle=\mathcal{B}(t)\delta\varphi(\bm{r},\bm{\hat{\Omega}},E,t)+\delta\mathcal{B}(t)\varphi_{c}(\bm{r},\bm{\hat{\Omega}},E)$
$\displaystyle+\mathcal{B}(t)\varphi_{c}(\bm{r},\bm{\hat{\Omega}},E)+\delta\mathcal{B}(t)\delta\varphi(\bm{r},\bm{\hat{\Omega}},E,t)\text{.}$
(8)

The term $\mathcal{B}(t)\varphi_{c}(\bm{r},\bm{\hat{\Omega}},E)$ vanishes,
since $\varphi_{c}$ is the critical (stationary) flux. Next, we apply the
‘orthodox’ linearization, where we assume that for small perturbations the
terms involving products of perturbed quantities can be neglected, namely,
$\delta\mathcal{B}\delta\varphi\approx 0$. The validity of this assumption has
been discussed by several authors [11, 37, 12], and has been recently
revisited for the case of mechanical vibrations [16, 17]. Upon rearranging the
two remaining terms, we obtain a Boltzmann-like transport equation for the
neutron noise which resembles a fixed-source problem:

$\mathcal{B}(t)\delta\varphi(\bm{r},\bm{\hat{\Omega}},E,t)=-\delta\mathcal{B}(t)\varphi_{c}(\bm{r},\bm{\hat{\Omega}},E)\text{,}$
(9)

the term $-\delta\mathcal{B}\varphi_{c}$ representing the ‘noise source’.

Since we are typically interested in the periodic solution to the noise
equation, we apply then a Fourier transform

$g(\omega)={\mathcal{F}}\left[g(t)\right](\omega)=\int_{-\infty}^{+\infty}e^{-i\omega
t}g(t)dt$ (10)

to Eq. (9), $\omega$ denoting the angular frequency, which yields the
linearized noise equation in the frequency domain:

$\mathcal{B}(\omega)\delta\varphi(\bm{r},\bm{\hat{\Omega}},E,\omega)=-\delta\mathcal{B}(\omega)\varphi_{c}(\bm{r},\bm{\hat{\Omega}},E)\text{.}$
(11)

The Fourier-transformed operators appearing in Eq. (11) are complex-valued and
read

$\mathcal{B}(\omega)=i\frac{\omega}{v}+\bm{\hat{\Omega}}\cdot\gradient+\Sigma_{t}(\bm{r},E)-\iint\nu_{s}(\bm{r},E^{\prime})\Sigma_{s}(\bm{r},E^{\prime})f_{s}(E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}})\differential
E^{\prime}\differential\bm{\hat{\Omega}}^{\prime}-\\\
\frac{1}{k_{\text{eff}}}\iint\nu_{f,p}(\bm{r},E^{\prime})\Sigma_{f}(\bm{r},E^{\prime})f_{f,p}(E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}})\differential
E^{\prime}\differential\bm{\hat{\Omega}}^{\prime}-\\\
\frac{1}{k_{\text{eff}}}\sum_{j}\frac{\lambda_{j}}{\lambda_{j}+i\omega}\iint\nu_{f,d}^{j}(\bm{r},E^{\prime})\Sigma_{f}(\bm{r},E^{\prime})f_{f,d}^{j}(E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}})\differential
E^{\prime}\differential\bm{\hat{\Omega}}^{\prime}\text{,}$ (12)

and

$\delta\mathcal{B}(\omega)=\delta\Sigma_{t}(\bm{r},E,\omega)-\iint\delta[\nu_{s}\Sigma_{s}f_{s}](\bm{r},E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}},\omega)\differential
E^{\prime}\differential\bm{\hat{\Omega}}^{\prime}-\\\
\frac{1}{k_{\text{eff}}}\iint\delta[\nu_{f,p}\Sigma_{f,p}f_{f,p}](\bm{r},E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}},\omega)\differential
E^{\prime}\differential\bm{\hat{\Omega}}^{\prime}-\\\
\frac{1}{k_{\text{eff}}}\sum_{j}\frac{\lambda_{j}}{\lambda_{j}+i\omega}\iint\delta[\nu_{f,d}^{j}\Sigma_{f,d}f_{f,d}^{j}](\bm{r},E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}},\omega)\differential
E^{\prime}\differential\bm{\hat{\Omega}}^{\prime}\text{.}$ (13)

For symmetry, it is sometimes convenient to write the
$\delta\Sigma_{t}(\bm{r},E,\omega)$ term in Eq. (13) in the same form as the
other terms:

$\delta\Sigma_{t}(\bm{r},E,\omega)=\\\
\iint\delta\Sigma_{t}(\bm{r},E^{\prime},\omega)\delta(E^{\prime}-E)\delta(\bm{\hat{\Omega}}^{\prime}-\bm{\hat{\Omega}})\differential
E^{\prime}\differential\bm{\hat{\Omega}}^{\prime}\text{,}$ (14)

which shows that $\delta\Sigma_{t}(\bm{r},E,\omega)$ is associated to a unit
yield and a ‘delta-copy’ spectrum
$\delta(E^{\prime}-E)\delta(\bm{\hat{\Omega}}^{\prime}-\bm{\hat{\Omega}})$.
The terms $\delta\Sigma_{t}(\bm{r},E,\omega)$ and
$\delta[\nu_{\alpha}\Sigma_{\alpha}f_{\alpha}](\bm{r},E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}},\omega)$ appearing in the
noise source involve the Fourier-transformed perturbation of the transport
operators, and as such their precise functional form is determined by the
specific type of perturbation of the noise model (e.g., mechanical vibrations
or oscillations). Regardless of the nature of the solver (deterministic or
Monte Carlo), these quantities must thus be prepared in advance, in the source
code, before solving the noise equation (11). In the next section we will
discuss the exact form of the noise source for two different types of
perturbations, and the Monte Carlo strategies that we have implemented to
sample them.

##  III Sampling the Neutron Noise Source

Sampling the neutron noise source requires the knowledge of the Fourier-
transformed perturbation operator $-\delta\mathcal{B}(\omega)$ and the
fundamental eigenmode $\varphi_{c}(\bm{r},\bm{\hat{\Omega}},E)$. The
fundamental eigenmode can be obtained by running the well-known power
iteration algorithm, which leads to a collection of neutrons whose flux is
precisely $\varphi_{c}(\bm{r},\bm{\hat{\Omega}},E)$, after a sufficiently
large numbers of inactive cycles. When convergence has been attained,
$-\delta\mathcal{B}(\omega)\varphi_{c}(\bm{r},\bm{\hat{\Omega}},E)$ can be
estimated by running one or more additional power iteration cycle, where the
noise source particles are sampled from the histories of the neutrons, and are
placed in a dedicated bank [28]. This sampling strategy is inspired from the
algorithm used to implement Monte Carlo kinetic simulations [38]. The explicit
form of $-\delta\mathcal{B}(\omega)$ depends as mentioned on the type of noise
source to be examined.

###  III.1 Cross Section Oscillations

First, we shall consider the simple case of cross section oscillations, which,
although somewhat artificial, is nonetheless important in that it will play a
key role in addressing the more involved case of mechanical vibrations,
covered in Sec. III.2. Cross section oscillations are relatively
straightforward, since they do not involve any spatial movement of the
materials. Suppose there is a spatial region $\mathcal{P}$ where the
macroscopic cross sections exhibit a periodic time dependence. We write the
cross sections as [30]

$\Sigma_{\alpha}(\bm{r},E,t)=\Sigma_{\alpha}(\bm{r},E)\left[1+\varepsilon_{\alpha}\sin(\omega_{0}t)\mathbbm{1}_{\mathcal{P}}(\bm{r})\right]\text{,}$
(15)

where $\varepsilon_{\alpha}$ is the amplitude of the perturbation, and

$\mathbbm{1}_{\mathcal{P}}(\bm{r})=\begin{cases}1&\text{for
}\bm{r}\in\mathcal{P}\\\ 0&\text{for
}\bm{r}\not\in\mathcal{P}\end{cases}\text{.}$ (16)

Upon taking the Fourier transform, the cross section perturbation in the
frequency domain is then

$\delta\Sigma_{\alpha}(\bm{r},E,\omega)=\\\
-i\pi\varepsilon_{\alpha}\Sigma_{\alpha}(\bm{r},E)\mathbbm{1}_{\mathcal{P}}(\bm{r})\left[\delta(\omega-\omega_{0})+\delta(\omega+\omega_{0})\right]\\\
=\Sigma_{\alpha}(\bm{r},E)\mathbbm{1}_{\mathcal{P}}(\bm{r})h_{\alpha}(\omega,\omega_{0},\varepsilon_{\alpha})\text{.}$
(17)

We assume that no other material properties display any time dependence. Such
a problem is admittedly academic, but provides an excellent benchmark for
comparison of codes and verification of algorithms. The noise source terms in
Eq. (13) can be thus written as

$\iint\nu_{\alpha}(\bm{r},E^{\prime})\delta\Sigma_{\alpha}(\bm{r},E^{\prime},\omega)\times\\\
f_{\alpha}(\bm{r},E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}})\varphi_{c}(\bm{r},\bm{\hat{\Omega}}^{\prime},E^{\prime})\differential
E^{\prime}\differential\bm{\hat{\Omega}}^{\prime}\text{.}$ (18)

By multiplying and dividing by $\Sigma_{\alpha}(\bm{r},E^{\prime})$ in the
integrand, and rearranging, we then have

$\iint\underbrace{\frac{\delta\Sigma_{\alpha}(\bm{r},E^{\prime},\omega)}{\Sigma_{\alpha}(\bm{r},E^{\prime})}}_{\text{complex
importance factor}}\times\\\
\underbrace{\nu_{\alpha}(\bm{r},E^{\prime})\Sigma_{\alpha}(\bm{r},E^{\prime})f_{\alpha}(\bm{r},E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}})}_{\text{standard production
rate}}\times\\\
\varphi_{c}(\bm{r},\bm{\hat{\Omega}}^{\prime},E^{\prime})\differential
E^{\prime}\differential\bm{\hat{\Omega}}^{\prime}\text{.}$ (19)

Eq. (19) can be handled by Monte Carlo codes with importance sampling, where
we use the standard sampling methods for the production rate, and multiply the
weight of the sampled noise particles corresponding to reaction $\alpha$ by
the complex factor
$\delta\Sigma_{\alpha}(\bm{r},E,\omega)/\Sigma_{\alpha}(\bm{r},E)=h_{\alpha}(\omega,\omega_{0},\varepsilon_{\alpha})\mathbbm{1}_{\mathcal{P}}(\bm{r})$.
To sample the noise source, action only needs to be taken when a neutron has a
collision within the region $\mathcal{P}$ exhibiting the cross section
oscillations. At this point, we need to separately sample the three components
of the noise source, that is, the copy term $\delta\Sigma_{t}\varphi_{c}$, the
fission term, and the scattering term. The noise source contribution from the
term $\delta\Sigma_{t}\varphi_{c}$ is sampled by creating an exact copy of the
particle and multiplying its weight by
$h_{t}(\omega,\omega_{0},\varepsilon_{t})$. For fission and scattering, a
nuclide is sampled from the material at $\bm{r}$, in the standard manner, and
we apply forced fission and implicit absorption. If the nuclide is fissile,
then a random variable $\xi\sim\mathcal{U}(0,1)$ is used to calculate the
number of fission noise particles to produce:

$n=\left\lfloor\frac{\nu_{f}\left(E\right)\sigma_{f}\left(E\right)}{\sigma_{t}\left(E\right)k_{\text{eff}}}+\xi\right\rfloor\text{,}$
(20)

where $\sigma_{t}(E)$ and $\sigma_{s}(E)$ are respectively the microscopic
total and scattering cross sections for the selected nuclide. Each fission
noise particle will inherit the weight of its parent, and have its energy and
direction sampled normally. The weights of these particles are multiplied by
$h_{f}(\omega,\omega_{0},\varepsilon_{f})$; furthermore, any delayed neutrons
will also have their weight multiplied by the complex yield
$\lambda_{j}/(\lambda_{j}+i\omega)$. Finally, noise particles from scattering
must be sampled; first we sample a scattering channel, and then we sample the
energy and direction of the scattered noise particle. In addition to
multiplying the noise particle’s weight by
$h_{s}(\omega,\omega_{0},\varepsilon_{s})$, it is also necessary to multiply
by the probability of scatter, $\sigma_{s}(E)/\sigma_{t}(E)$, where
$\sigma_{s}(E)$ is the microscopic scattering cross sections for the selected
nuclide.

###  III.2 Mechanical Vibrations

Mechanical vibrations are inherently different from the previous case of cross
section oscillations, due to the fact that they involve material displacements
as a function of time. Assuming that the displacement occurs along one spatial
dimension, we will have material $M_{L}$ at position $x<x_{0}$ and material
$M_{R}$ at position $x>x_{0}$ before the perturbation is applied. In
stationary conditions, the material as a function of position is thus
described by the function

$M(x)=\begin{cases}M_{L}&\text{for }x<x_{0}\\\ M_{R}&\text{for
}x>x_{0}\end{cases}\text{.}$ (21)

Now, if a perturbation is introduced into this system in the form of a
sinusoidal vibration of the interface between the two materials, with angular
frequency $\omega_{0}$ and amplitude $\varepsilon>0$, then the material found
at any given position at a given position in the perturbed region is also a
function of time. For positions where $\absolutevalue{x-x_{0}}<\varepsilon$,
the function which determines the material at position $x$ and time $t$ is
written as

$M(x,t)=M(x)+\\\ \Delta
M\left[H(x-x_{0})-H(x-x_{0}-\varepsilon\sin(\omega_{0}t))\right]\text{,}$ (22)

where $\Delta M=M_{L}-M_{R}$, and $H(x)$ is the Heaviside function [16]. This
function always returns either $M_{L}$ or $M_{R}$. Based on the material
perturbation $\delta M=M(x,t)-M(x)$ across the interface, the corresponding
perturbation of the operator kernels
$\delta[\nu_{\alpha}\Sigma_{\alpha}f_{\alpha}]$ must then be

$\delta[\nu_{\alpha}\Sigma_{\alpha}f_{\alpha}](x,E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}},t)=\Delta[\nu_{\alpha}\Sigma_{\alpha}f_{\alpha}]\times\\\
\left[H(x-x_{0})-H(x-x_{0}-\varepsilon\sin(\omega_{0}t))\right]\text{,}$ (23)

where

$\Delta[\nu_{\alpha}\Sigma_{\alpha}f_{\alpha}]=\\\
\nu_{\alpha}^{L}(E)\Sigma_{\alpha}^{L}(E)f_{\alpha}^{L}(E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}})-\\\
\nu_{\alpha}^{R}(E)\Sigma_{\alpha}^{R}(E)f_{\alpha}^{R}(E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}})\text{.}$ (24)

The Fourier transform of Eq. (23) has been shown to be [16]

$\delta[\nu_{\alpha}\Sigma_{\alpha}f_{\alpha}](x,E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}},\omega)=\\\
\Delta[\nu_{\alpha}\Sigma_{\alpha}f_{\alpha}]\bigg{\\{}c_{0}(x,x_{0})\delta(\omega)+\\\
\sum_{n=1}^{\infty}c_{n}(x,x_{0})\left[\delta(\omega-n\omega_{0})+\delta(\omega+n\omega_{0})(-1)^{n}\right]\bigg{\\}}\text{,}$
(25)

with

$c_{n}(x,x_{0})=\frac{2}{n}\sin(n\arcsin(\frac{x-x_{0}}{\varepsilon}))e^{-in\pi/2}$
(26)

for $n\geq 1$, and

$c_{0}(x,x_{0})=\begin{cases}\pi-2\arcsin(\frac{x-x_{0}}{\varepsilon})&\text{for
}x\geq x_{0}\\\ \\\ -\pi-2\arcsin(\frac{x-x_{0}}{\varepsilon})&\text{for
}x<x_{0}\end{cases}\text{.}$ (27)

Immediately apparent from Eq. (25) is that, while the noise source for cross
section oscillations is monochromatic, the noise source for mechanical
vibrations contains an infinite number of discrete harmonics at $n\omega_{0}$,
all multiple of the perturbation frequency $\omega_{0}$. Additionally, there
is a non-trivial spatial dependence indicated by Eq. (26). All even harmonics
are purely real and have even spatial parity about $x_{0}$, whereas all odd
harmonics are purely imaginary and have odd parity. We will now consider
possible methods to sample the noise source for mechanical vibrations in a
Monte Carlo code.

####  III.2.1 Approximate Noise Source Sampling

From Eq. (23), we see that the noise source is caused by a difference of the
product of the cross section, reaction yield, and scattering law. In the
initial approach implemented in the noise solver of the TRIPOLI-4® Monte Carlo
code, developed at CEA [39], it was assumed that

$\Delta[\nu_{\alpha}\Sigma_{\alpha}f_{\alpha}]\approx\nu_{\alpha}f_{\alpha}\Delta\Sigma_{\alpha}(x,x_{0},\omega)\text{,}$
(28)

where only the difference in the macroscopic cross sections was considered. It
must be noted that, while $\Delta\Sigma_{\alpha}$ requires evaluating the
cross section on both sides of the interface, $\nu_{\alpha}$ and $f_{\alpha}$
were both evaluated only for the material at the location where the collision
occurred. This was partially done for simplicity, as evaluating $f_{\alpha}$
on both sides of the interface to obtain
$\Delta[\nu_{\alpha}\Sigma_{\alpha}f_{\alpha}]$ could pose particular
challenges222While the data representations used in Monte Carlo codes allow
for easy sampling from the distribution $f_{\alpha}$, they often do not
facilitate the evaluation of the distribution $f_{\alpha}$ corresponding to a
given argument..

Figure 1: Depicted is the operation of copying noise particles from fission in
the fuel, and placing them in the moderator. The line at $x_{0}$ is the
boundary between the fuel and the moderator for the static problem. Dashed
lines at $x_{0}-\varepsilon$ and $x_{0}+\varepsilon$ represent the minimum and
maximum positions of the interface with the vibration. The particle in the
fuel is the sampled fission noise particle. An exact copy of this fission
noise particle (identical direction, energy, and weight) is placed at the
symmetric position within the moderator, corresponding to an $x$-coordinate of
$x_{0}+d$. The weight of the copied particle in the moderator must be
multiplied by $(-1)^{n+1}$, to account for the parity of the coefficients
$c_{n}(x,x_{0})$.

At first glance, it would appear that Eq. (28) could be handled with
importance sampling, as was done for cross section oscillations, where one
only needs information from the material at which the collision occurred, and
this is the main advantage of the proposed approximation. Sadly, further
assumptions are required. As a tangible example, let us consider the case of a
fuel pin vibrating in water moderator. Unlike for the case of oscillating
cross sections, where all of the nuclides participating in the perturbation
are present at the sampling site in the unperturbed system, in Eq. (28)
isotopes from the fuel participate in the perturbation in the water and vice
versa. It should therefore be possible to sample fission noise particles from
collisions occurring in the water surrounding the vibrating pin, which would
not be allowed with the importance sampling scheme from oscillations. A clever
method was devised to circumvent this issue in TRIPOLI-4®: when a fission
noise particle is sampled in the fuel region, a copy is made and then
translated across the boundary of the vibration, to the symmetric location
relative to the vibrating interface (taking into account the parity of the
given harmonic). So long as the flux is approximately constant over the
distance of the translation, this is likely a reasonable approximation,
allowing for the placement of fission noise source particles in the water.
This approximation of the placement of fission noise particles in the
moderator is depicted in Fig. 1. Otherwise, the scheme implemented in
TRIPOLI-4® samples noise source particles for vibrations in the same manner as
for cross section oscillations. TRIPOLI-4® therefore also ignores the
possibility of creating scatter noise particles from isotopes present on the
other side of the interface (i.e. a noise particle born from a scattering
event with a fuel isotope, when having a collision in the water).

####  III.2.2 Exact Noise Source Sampling

The noise equation is typically written using macroscopic material cross
sections, yields, and energy-angle distributions. Each of the factors in the
source term is actually a weighted average of the cross sections, yields and
energy-angle distributions of each individual isotope present in the material.
It is therefore possible to rewrite Eq. (3) for mechanical vibrations as

$\nu_{\alpha}(\bm{r},E,t)\Sigma_{\alpha}(\bm{r},E,t)f_{\alpha}(\bm{r},E\to
E^{\prime},\bm{\hat{\Omega}}\to\bm{\hat{\Omega}}^{\prime},t)=\\\
\sum_{k}N_{k}(\bm{r},t)\nu_{k,\alpha}(\bm{r},E)\sigma_{k,\alpha}(\bm{r},E)f_{k,\alpha}(\bm{r},E\to
E^{\prime},\bm{\hat{\Omega}}\to\bm{\hat{\Omega}}^{\prime})\text{.}$ (29)

Here, $N_{k}$ is the concentration of nuclide $k$, and $\sigma_{k,\alpha}$ is
the microscopic cross section of nuclide $k$ for reaction channel $\alpha$.
For an individual nuclide, the microscopic cross sections, reaction yields and
energy-angle distributions are time-independent. However, the macroscopic
material-averaged quantities become time-dependent through the concentration
$N_{k}(\bm{r},t)$ of nuclides in the perturbed region, which is clearly a
time-dependent quantity. In Monte Carlo codes, each individual nuclide is
modelled explicitly, and the ensemble material averages naturally arise from
neutrons undergoing collisions with all of the nuclides present in a given
material. This suggests a new sampling strategy for the noise source due to
mechanical vibrations. The concentration of nuclides in the perturbed region
can be decomposed as

$N_{k}(\bm{r},t)=N_{k}(\bm{r})+\delta N_{k}(\bm{r},t)\text{,}$ (30)

allowing us to write each term of the noise source as

$\iint\sum_{k}\delta
N_{k}(\bm{r},t)\nu_{k,\alpha}(\bm{r},E^{\prime})\sigma_{k,\alpha}(\bm{r},E^{\prime})\times\\\
f_{k,\alpha}(\bm{r},E^{\prime}\to
E,\bm{\hat{\Omega}}^{\prime}\to\bm{\hat{\Omega}})\varphi_{c}(\bm{r},\bm{\hat{\Omega}}^{\prime},E^{\prime})\differential
E^{\prime}\differential\bm{\hat{\Omega}}^{\prime}\text{.}$ (31)

The form of $\delta N_{k}(\bm{r},t)$ depends again on the material interface
behavior given in Eq. (22), and therefore leads to a similar form for the
Fourier transform $\delta N_{k}(\bm{r},\omega)$.

To sample Eq. (31), one might be inclined to multiply and divide by
$N_{k}(r)$, similar to the approach taken for oscillations, so to convert Eq.
(31) into a complex-weighted particle production rate. Sadly, such a strategy
is not possible for the case of vibrations. Consider the case of a UO2 fuel
fuel pin vibrating in water. When inside the perturbed region in the water
surrounding the pin, $\delta N_{{}^{235}\text{U}}\neq 0$, but
$N_{{}^{235}\text{U}}=0$, and it is therefore impossible to divide by
$N_{{}^{235}\text{U}}(\bm{r})$. Our proposed solution to this problem is to
define a new fictitious material which is only used for the sampling of the
noise source. This material must contain all nuclides which are present in the
two materials involved in the vibration, and all of their concentrations must
be nonzero. The exact choice of the concentration for each isotope is
arbitrary. In this work, we decided that if a nuclide is only present in the
left or the right material (but not both), then its concentration in the
material where it is present is taken for use in the fictitious material. If
an isotope exists in both the left and right material, then the average
concentration was used. The question of the optimal choice of these
concentrations is left for future work. We shall denote the isotope
concentrations of this fictitious material as $N^{*}_{k}(\bm{r})$;
correspondingly, the fictitious material has a total cross section of

$\Sigma^{*}_{t}(\bm{r},E)=\sum_{k}N^{*}_{k}(\bm{r})\sigma_{k,t}(E)\text{.}$
(32)

It is then possible to multiply and divide Eq. (31) by this fictitious
concentration, and perform importance sampling. The collision kernel which
should be used to sample the noise source is therefore

$\frac{\Sigma^{*}_{t}(\bm{r},E)}{\Sigma_{t}(\bm{r},E)}\sum_{k}\frac{\delta
N_{k}(\bm{r},\omega)}{N^{*}_{k}(\bm{r})}\frac{N^{*}_{k}(\bm{r})\sigma_{k,t}(E)}{\Sigma^{*}_{t}(\bm{r},E)}\times\\\
\sum_{m}\frac{\sigma_{k,m}(E)}{\sigma_{k,t}(E)}\nu_{k,m}(E)f_{k,m}(E\to
E^{\prime},\bm{\hat{\Omega}}\to\bm{\hat{\Omega}}^{\prime})\text{.}$ (33)

Now, the noise source sampling strategy is as follows. During power iteration,
when a neutron undergoes a collision in the region of a vibration, a new
fictitious material must be constructed. Nuclide $k$ is sampled with
probability $N^{*}_{k}(\bm{r})\sigma_{k,t}(E)/\Sigma^{*}_{t}(\bm{r},E)$, and
will be used only for the sampling of the noise source particles. At this
point, the sampling of the noise source is congruent to the approach used for
cross section oscillations. We have effectively reduced the problem of
sampling the vibration source to the problem of sampling an oscillating
(space-dependent) source. If nuclide $k$ is fissile, then fission noise
particles must be sampled. Next, the noise particles born from scattering are
sampled. A scattering channel for the nuclide can then be sampled nominally,
along with the outgoing energy and direction. All the sampled noise source
particles must have their weights multiplied by the complex factor

$\frac{\Sigma^{*}_{t}(\bm{r},E)\delta
N_{k}(\bm{r},\omega)}{\Sigma_{t}(\bm{r},E)N^{*}_{k}(\bm{r})}\text{.}$ (34)

Finally, a copy of the incident particle is also created and added to the
noise source, with its weight multiplied by
$\delta\Sigma_{t}(\bm{r},E,\omega)/\Sigma_{t}(\bm{r},E)$. The fictitious
material is not actually needed for producing the copy noise particle, which
can be handled in an identical manner to the approach taken earlier in Eq.
(14).

It should be noted that for realistic problems vibration regions could
overlap. For the case of a fuel pin with cladding, there would effectively be
two material interfaces (fuel/cladding and cladding/water): if the amplitude
of the vibration is larger than the thickness of the cladding, the two
vibration regions will overlap. This will have the effect of creating three
effective regions: one where a fuel-cladding mixture is visible, one where a
fuel-cladding-water mixture is visible, and another where a cladding-water
mixture is visible. Since the effects of vibrations can be combined linearly,
treatment of this situation is straightforward. One just needs to ensure that
all nuclides of all vibrations at a given position are included in the
fictitious material, and then sum all the contributions to $\delta
N_{k}(\bm{r},\omega)$ coming from the vibration of each interface. In our
implementation we chose to construct the fictitious material only when it is
known what vibration regions are influencing the noise sampling at the given
collision site. This avoids determining all possible region overlaps and
material combinations in advance, which might be rather involved depending on
the geometric form of each vibration.

##  IV Transport of Noise Particles

Once the noise source has been sampled as shown in the previous section, noise
particles must be transported according to the stochastic rules defined by the
operator $\mathcal{B}(\omega)$, in order to solve the fixed-source problem in
Eq. (11). For both oscillations and mechanical vibrations, the structure of
the noise source term is

$-\delta{\cal B}(\omega)\varphi_{c}=-\sum_{n}\delta{\cal
B}_{n}\varphi_{c}\delta(\omega-n\omega_{0}),$ (35)

where $\delta{\cal B}_{n}=\delta{\cal B}(n\omega_{0})$ is the perturbation
operator evaluated at frequency $n\omega_{0}$. (The case of oscillation is
special in that the source contains only a single frequency.) Correspondingly,
we therefore expect the noise field to be of the form

$\delta\varphi({\bf r},{\bm{\Omega}},E,\omega)=\sum_{n}\delta\varphi_{n}({\bf
r},{\bm{\Omega}},E)\delta(\omega-n\omega_{0}),$ (36)

where $\delta\varphi_{n}({\bf r},{\bm{\Omega}},E)$ is the solution of the
Fourier-transformed noise equation corresponding to the noise source component
at frequency $n\omega_{0}$, namely

${\cal B}_{n}\delta\varphi_{n}=-\delta{\cal B}_{n}\varphi_{c},$ (37)

with ${\cal B}_{n}={\cal B}(n\omega_{0})$ the noise operator evaluated at the
discrete frequencies $n\omega_{0}$ of the source. This yields an infinite
system of fully decoupled linear equations for the noise components
$\delta\varphi_{n}$. For negative frequencies where $n<0$, we make use of the
fact that $\delta\varphi_{-n}=\delta\varphi_{n}^{\dagger}$, the symbol
$\dagger$ indicating the complex conjugate. The null frequency at $n=0$
represents the the time-averaged effect on the static flux, as the
perturbation will introduce a change of reactivity into the system, which is
typically neglected in the orthodox linearization approach [16].

The operator $\delta\mathcal{B}$ is complex-valued, indicating that the noise
source particles must therefore carry complex statistical weights, and each
component can be positive or negative. The noise transport operator
$\mathcal{B}$ is also complex-valued, and differs from the standard Boltzmann
operator in three places, which therefore alter the regular flight and
collision sampling procedures. We first notice that the fission production
must be divided by $k_{\text{eff}}$, which must be known before starting the
noise simulation. The second modification is that delayed neutrons must have
their weights multiplied by the complex yield factor
$\lambda_{j}/(\lambda_{j}+i\omega)$, depending on the angular frequency
$\omega$ and on the precursor family $j$. Both of these changes are minor, and
straightforward to implement in most Monte Carlo codes. More problematic is
the fact that $\mathcal{B}$ involves an effective total cross section which is
complex:

$\Sigma_{t}(\bm{r},E)+i\frac{\omega}{v}\text{.}$ (38)

Two different methods have been proposed to address this peculiarity. Yamamoto
proposed to change the particle’s weight continuously along the flight, using
a complex exponential transform involving also the modification of the track-
length estimator [20]. In this work, we use an alternative implementation
proposed by Rouchon et al., which adds the real term

$\eta\frac{\omega}{v}\delta\varphi$ (39)

to both sides of Eq. (11), where $\eta$ is real and has the same sign as
$\omega$ [21]. Doing so effectively adds a new copy reaction with cross
section $\Sigma_{\omega}(E)=\eta\omega/v$, having a complex yield
$\nu_{\omega}=(\eta-i)/\eta$, and the effective total cross section becomes

$\tilde{\Sigma}(\bm{r},E)=\Sigma_{t}(\bm{r},E)+\Sigma_{\omega}(E)\text{,}$
(40)

which is a real positive quantity. Flight distances of noise particles are now
sampled using the effective total cross section $\tilde{\Sigma}$. At each
collision, a copy of the noise particle is made with probability
$\Sigma_{\omega}/\tilde{\Sigma}$, which has its complex weight multiplied by
$\nu_{\omega}$. This is consistent with treatment of $\delta\Sigma_{t}$ given
in Eq. (14). In general, $\Sigma_{t}$ must be replaced by $\tilde{\Sigma}$
throughout the transport algorithms. Otherwise, transport algorithms are left
unchanged. Variance reduction techniques such as implicit capture and roulette
can still be used. In the method proposed by Rouchon et al., implicit capture
is performed by multiplying the complex weight by $\Sigma_{s}/\tilde{\Sigma}$.
For the roulette, we follow the algorithm proposed by Yamamoto, where the real
and imaginary weight components undergo roulette independently [21, 20]. The
choice of $\eta=1$ was used, as this has been previously demonstrated to yield
good performance for most frequencies of interest for reactor physics [21].

Each component of the complex weight of noise particles can be positive or
negative; correspondingly, the noise field will be estimated through the sum
of these positive and negative contributions, which causes a very large
variance in the scores. This phenomenon has been observed in cross section
oscillation problems [30], but is exacerbated by the case of vibrating fuel
pins [31]. On the one hand, as the noise source for cross section oscillations
is constant within the vibration region, the particles of differing sign are
only produced through delayed fission, which is a weak effect in most
realistic systems. During noise transport, sign changes to particle weights
can only occur through the copy channel, and delayed fission. Both of these
reaction channels are typically weak, and this therefore leads to minimal
appearances of particles of differing sign components. On the other hand, as
is depicted in Sec. VI.1, and previously shown by Zoia et al. [16], the noise
source of a vibrating fuel rod has a positive contribution on one side, and a
negative contribution on the other, leading to nearly equal quantities of
positive and negative noise source particles at the beginning of the problem.
The resulting noise field is effectively the difference of two nearly equal
stochastic quantities, which leads to catastrophic convergence issues for the
statistical uncertainty. From these preliminary investigations, it is evident
that the development of specialized variance reduction techniques is required
if we want to obtain usable neutron noise results via Monte Carlo methods.

##  V Weight Cancellation for Variance Reduction

Monte Carlo transport problems involving a mixture of positive and negative
particles are not unique to neutron noise, and have appeared repeatedly in the
literature for e.g. performing delta tracking without a majorant cross section
[40, 41], for the calculation of the higher harmonics of $k$-eigenvalue
problems [42, 43], and for the computation of critical buckling [44, 45].
Based on his pioneering work on higher harmonics and critical buckling,
Yamamoto’s initial implementation for solving the noise equation by Monte
Carlo methods utilized an approximate weight cancellation technique, as it was
determined to be necessary when trying to solve the noise equation for
frequencies outside of the ‘plateau region’
$\lambda<\omega<\lambda+\beta_{\text{eff}}/\Lambda_{\text{eff}}$ [20]. When
solving the noise equation for frequencies within the plateau region, weight
cancellation was determined to be unnecessary, and was abandoned in Yamamoto’s
later implementation done in MCNP4C, due to the complexity of incorporating
such changes in a production-level Monte Carlo code [46]. Furthermore,
Yamamoto’s initial implementation using weight cancellation was not able to
estimate the variance of scores, and it was therefore impossible to determine
what effect, if any, weight cancellation had on the estimation of the noise
field, beside enabling the calculations to run in finite time [20]. The method
developed by Rouchon et al. is able to estimate the variance of the noise
field; previous results from simpler case-studies showed that convergence
within the plateau region could be achieved without weight cancellation, which
was therefore not implemented in the code [21]. Nonetheless, the convergence
difficulties mentioned above for the case of vibrations has led us to begin an
examination of the possible benefits of approximate or exact weight
cancellation. Our recent work has demonstrated that the use of weight
cancellation can provide a very large improvement in computational efficiency
for the case of neutron noise induced by cross section oscillations, leading
to an improvement by a factor of 90-100 on the figure of merit (FOM) when
compared to the implementation in TRIPOLI-4® [36]. In this section, we show
how weight cancellation can be performed on neutron noise problems for the
case of vibrations, and illustrate the computational gains which can be
achieved through the use of the algorithm.

###  V.1 When to Perform Cancellation

Previous investigations in the context of weight cancellation for delta
tracking methods have shown that it is most efficient to perform cancellation
on the fission source [34]: fission particles are unique in that their angular
distribution is almost always considered to be perfectly isotropic, and their
energy is almost independent of the incident energy. This effectively reduces
the dimensionality of weight cancellation, which makes exact regional
cancellation methods more efficient (in terms of the amount of weight being
cancelled), and reduces the bias inherent to approximate cancellation methods
[34, 35].

Most of our previous work on cancellation has been applied to eigenvalue
problems solved by using the power iteration, where one naturally has access
to the fission source between subsequent fission generations, providing the
ideal opportunity to apply the weight cancellation algorithm [34]. However,
neutron noise equations are effectively fixed-source problems, where source
particles are sampled and then distributed to all the available computational
nodes for the simulation. Since fixed-source problems must be always
subcritical, all simulated particles eventually die: once a node receives its
packet of particles to transport, it never needs to communicate with the other
nodes again, until it receives the next packet of particles for the subsequent
batch. Any fission particles born during a batch are typically added to the
secondary particle bank, so that they can be transported later on: the node
will eventually process all of the generated fission noise particles (because
of sub-criticality), completing the batch. A depiction of the MPI
communication scheme for this standard fixed-source algorithm is provided in
Figure 2. Since all of the fission particles are just added to the bank of
particles to be transported during the current batch, there is never a step in
the algorithm where we have access to the entire fission source during the
computation, contrary to what happens in power iteration schemes.

Figure 2: Standard MPI communication scheme for a single fixed-source batch.
Figure 3: Proposed MPI communication scheme where a single fixed-source batch
is broken into inner fission generations. While only three inner generations
are depicted here to demonstrate the concept, there will generally be hundreds
or thousands of inner generations per noise batch.

However, having access to the entire fission source is essential to maximizing
the efficiency of weight cancellation: if not all of the particles
participate, less weight will be canceled. We have therefore implemented a
different fixed-source iteration algorithm for noise problems, in order to
gain access to the entire fission source. The scheme presented here was
inspired by the decomposition proposed by Yamamoto, although differing in
implementation [20, 46]. For a single batch, the initial source particles are
sampled, and then sent to all the available nodes. Each node will transport
all of the assigned initial particles, and will store all of the produced
fission particles in a separate bank. Once it has transported the initial
assigned packet of particles, each node will send its bank of fission
particles back to the master node. Now, with all of the fission particles on
the master node, it is possible to perform weight cancellation on the entire
fission source. Once the cancellation operation is complete, the new fission
particles are redistributed to the nodes, where they will be transported. Any
resulting fission particles will again be stored in a separate bank; once each
node has finished, the fission particles are again sent to the master node,
where cancellation can occur again. This practically breaks the fixed-source
batch into ‘inner generations’, as depicted in Fig. 3. This MPI communication
scheme is nearly identical to that used in standard power iteration
implementations, with the exception that all of the particles will eventually
die out, signaling the end of the fixed-source batch.

###  V.2 Approximate Cancellation

In this work, we will consider both exact and approximate weight cancellation
methods. We will first examine approximate cancellation, based on the
technique initially proposed by Zhang et al., which is conceptually simple,
easy to implement in a Monte Carlo code, and results in very good cancellation
efficiency [47]. For this purpose, a rectilinear mesh is imposed on top of the
problem domain, and the particles which will undergo cancellation are then
sorted into this mesh, based on their position (and possibly energy, if
desired). The average weight of all particles in each mesh element is
calculated, and this average weight is then assigned to each of the particles
in the element. Summing the particle weights to compute the average is what
effectively leads to cancellation, as the total weight will be the sum of
positive and negative components. While this approach imposes a bias on the
resulting fission source, the bias can be made minimal by taking a
sufficiently fine mesh. However, refining the mesh is a double-edged sword:
the bias will be reduced, but the efficiency of cancellation will be reduced
as well, as there will be fewer particles present in each mesh element.
Despite this, our previous work has shown approximate cancellation to be much
more efficient at cancelling weight than the exact version (treated in Sec.
V.3), and we have have not observed any measurable bias for a reasonably
refined mesh [35].

###  V.3 Exact Cancellation

Exact regional cancellation was originally proposed by Booth and Gubernatis,
for use in calculating the higher harmonics of $k$-eigenvalue problems [48].
The initial formulation was limited to one-dimensional configurations, but we
have previously expanded this formalism to three-dimensional multi-group
settings. This cancellation technique is applied on the fission source, and in
practice can be used in combination with delta tracking [49] or negative-
weighted delta tracking [40, 41] transport methods: we have provided a
successful numerical demonstration of exact regional cancellation for a
modified version of the C5G7 reactor physics benchmark [34]. Fission particles
are sorted into geometric volumes, referred to as cancellation regions. Our
previous work used a regular rectilinear mesh, with each mesh element acting
as a cancellation region. Fission particles must store the position of their
parent’s previous collision ($\bm{r}^{\prime}$), their parent’s energy
($E^{\prime}$), and their parent’s penultimate direction
($\bm{\hat{\Omega}}^{\prime\prime}$), in addition to their own phase space
coordinates $(\bm{r},\bm{\hat{\Omega}},E)$. From this information, the fission
density function can be calculated, which for the case of delta tracking reads

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

where $\Sigma_{\text{maj}}(E^{\prime})$ is the majorant cross section used in
the delta tracking algorithm [49, 50, 34]. Using this fission density
function, each fission particle $i$ (having a weight $w_{i}$) is split into
two portions: a point-wise part with weight

$w_{p,i}=\frac{\zeta(\bm{r}|\bm{r}^{\prime}\bm{\hat{\Omega}}^{\prime\prime},E^{\prime})-\beta_{i}}{\zeta(\bm{r}|\bm{r}^{\prime}\bm{\hat{\Omega}}^{\prime\prime},E^{\prime})}\text{,}$
(42)

and a uniform part with weight

$w_{u,i}=\frac{\beta_{i}}{\zeta(\bm{r}|\bm{r}^{\prime}\bm{\hat{\Omega}}^{\prime\prime},E^{\prime})}\text{.}$
(43)

Here $\beta_{i}$ is a cancellation parameter, which can take any value, so
long as it does not depend on $(\bm{r},\bm{\hat{\Omega}},E)$ [48, 34, 51]. The
point-wise portion represents the piece of the fission particle which must be
placed exactly at $\bm{r}$, and the uniform portion represents the piece which
can be uniformly distributed across the entire cancellation region. A thorough
discussion on how to determine an appropriate value of $\beta_{i}$, which
ultimately affects the cancellation efficiency, is provided in Refs. 34 and
51. For each cancellation region $\mathcal{R}$, we iterate over all of the
particles in the region, summing their uniform parts

$U_{\mathcal{R}}=\sum_{i\in\mathcal{R}}w_{u,i}\text{.}$ (44)

We then set all of the fission particle weights to be equal to their point-
wise portion ($w_{i}:=w_{p,i}$). Lastly, we must sample fission particles with
a total weight equal to the remaining weight $U_{\mathcal{R}}$, which will be
uniformly distributed within $\mathcal{R}$. The exact cancellation algorithm
can be applied as is to neutron noise problems: the only required modification
concerns the number of uniform particles to be added to each region, which for
particles carrying real and imaginary statistical weights reads

$n=\left\lfloor\text{max}\big{(}\absolutevalue{\Re{U_{\mathcal{R}}}},\absolutevalue{\Im{U_{\mathcal{R}}}}\big{)}\right\rfloor\text{,}$
(45)

each with a weight of $w=U_{\mathcal{R}}/n$.

Compared to approximate cancellation, exact regional cancellation is much more
difficult to implement in a Monte Carlo code. Additionally, it results in much
less total weight cancellation as well, which can lead to only minimal
performance improvements in problems which require cancellation [35, 34, 51].
However, the exact regional method aligns with the Monte Carlo philosophy, in
that no approximation nor bias is introduced in the sampling scheme.

##  VI Simulation Results

The Monte Carlo noise equation solution method devised by Rouchon et al. was
implemented in the development version of TRIPOLI-4® [21, 39]. Since
TRIPOLI-4® is a mature general-purpose code with about 400 kSLOC (source lines
of code), implementing all of the changes necessary to perform the exact noise
source sampling technique outlined in Sec. III.2.2, or the cancellation
methods described in Sec. V, would require an exceptionally large rewriting
effort. Therefore, all of the methods described in the paper have been
implemented in a multi-group Monte Carlo mini-app called MGMC, with
approximately 13 kSLOC and thus allowing for quick implementation and testing
of new transport algorithms [52]. MGMC is able to perform transport in general
3D geometries, composed of surface-based volumes, universes, and lattices.
Shared-memory parallelism is implemented via OpenMP, and distributed-memory
parallelism is implemented via MPI. Basic quantities such as flux and reaction
rates can be scored using either track-length or collision estimators, across
regular rectilinear meshes, and are written to binary NumPy files, for easy
data processing in Python [53]. These features make MGMC entirely
representative of a larger production-level code; therefore, it can be assumed
that any improvement in calculation efficiency observed in MGMC would be
similar to the improvements one could expect from implementing the same
methods in a production-level code, such as TRIPOLI-4®. Delta tracking and
weight cancellation have previously been added to MGMC, along with neutron
noise transport [34, 36]. Newly implemented for this work was the sampling the
noise source corresponding to mechanical vibrations. MGMC has been released as
free software, and is available under the CeCILL-v2.1 license [52].

In order to illustrate the novel noise source sampling and the weight
cancellation methods implemented in MGMC, in this work we have selected a
benchmark problem that has been previously used to compare the results of
several different neutron noise solvers [30, 31]. It is a 2D, 2 group problem,
in the form of a reflected 17x17 fuel assembly, with square pins. A diagram of
this system is presented in Fig. 4. The assembly pitch is
$1.26\text{\,}\mathrm{cm}$, with fuel pins having a side length of
$0.7314\text{\,}\mathrm{cm}$. A $0.08\text{\,}\mathrm{cm}$ water blade
surrounds the entire assembly. The neutron noise is induced by a vibrating
fuel pin (marked in pink in Fig. 4), subject to a sinusoidal displacement
along the x-axis, with an angular frequency of
$\omega_{0}=2\pi\,$\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}$$ and an amplitude
of $\varepsilon=$0.2\text{\,}\mathrm{cm}$$, as per benchmark specifications.

Since exact regional weight cancellation can only be performed with delta
tracking, all simulations were run with this transport method, to ensure
fairness in the comparison of simulation runtimes. All simulations began with
$10^{6}$ particles in power iteration, distributed uniformly within the
assembly, all in the first group. The first $13$ generations were discarded to
allow for source convergence. Between each noise batch, three extra power
iteration inactive generations were performed, to ensure proper decorrelation
of the noise source between replicas. Simulations were set to run for either
$10^{4}$ noise batches, or a maximum run time of 2 days (whichever came
first). Simulations were run on a computing cluster at CEA, each run using 16
MPI processes, and 32 OpenMP threads per MPI process. The noise field was
estimated with a track-length estimator in all simulations. The real and
imaginary components of the field are scored in separate tallies. Breaking the
scoring of complex quantities into a real and imaginary component facilitates
the addition of neutron noise simulations to existing codes.

For the purpose of code-to-code comparison and verification, our Monte Carlo
simulations were compared to the results generated by the deterministic noise
solver previously implemented in APOLLO3® [54, 55]. The noise equation solver
has been added to the IDT lattice solver of APOLLO3®, which uses the method of
short characteristics in conjunction with the discrete ordinates method [55].
To accomplish this, the standard iteration loops are applied to the complex
fission source and scattering source, with the addition of an iteration loop
between the real and imaginary components. Thus, the standard one-group
transport solver methods can be used, and one can consequently benefit from
all numerical methods already implemented in APOLLO3® [56].

Figure 4: Depiction of the 2D reflected assembly for the neutron noise
benchmark problem proposed in [30]. The pink fuel pin (located two cells down
and two cells left from the center) experiences a small sinusoidal vibration
to the left and right. All other fuel pins are unperturbed.

To compare the noise source and noise field between codes, these quantities
were normalized by the estimated value of the static flux at the center of the
assembly in group 2, for the respective code.

###  VI.1 Sampling of the Noise Source

Figure 5: Real and imaginary components of the noise source, for the first
harmonic.

We first compare the first harmonic of the noise source, sampled with the
constant flux approximation and the exact method, against the noise source
obtained from APOLLO3®. By default, the deterministic noise solver of APOLLO3®
computes the harmonics of the exact noise source (Sec. III.2.2). The real and
imaginary components of the source are shown in Fig. 5. The real component of
group 2 is zero, as this portion of the source only comes from delayed
neutrons, which according to the benchmark model are all born in group 1.
Immediately apparent is the fact that the amplitude of the source in the
constant flux approximation is much smaller than the true source, for the
imaginary component. The approximate source sampling also has an asymmetry,
where the source has a slightly larger amplitude in the fuel pin than in the
water. Both of these effects are primarily due to the approximate source
sampling method ignoring the difference in the scattering law in the noise
source. This is evident upon examination of the real component in group 1,
which is purely due to delayed fission: very good agreement is observed
between the approximate method and APOLLO3®. Since the fission component
appears to have the right amplitude and form in the real component, it can be
determined that the discrepancy in the imaginary component is primarily due to
the inadequate treatment of the scattering laws in the approximate method. The
agreement observed in the real component also indicates that making a copy of
the fission noise source particles in the fuel and moving them into the water
is an adequate approximation for this system, where the flux is relatively
constant over the perturbed region. Other numerical investigations concerning
systems with larger flux gradients (such as the one-dimensional rod model
investigated in [16]) show that this approximation might become inappropriate
and lead to large differences in the resulting noise source.

The exact noise source sampling method has excellent agreement with the
APOLLO3® results, for all groups and components. In the imaginary components
shown in Fig. 5, the exact source has the correct amplitude and shape when
compared to APOLLO3®. By using the same fictitious material to sample the
noise source on both sides of the vibration interface, it is possible to
eliminate the discontinuity which was observed on in the approximate method.
This demonstrates that our proposed exact source sampling method is not just
correctly treating the perturbation in the macroscopic cross sections, but
also in the scattering laws.

Figure 6: Real and imaginary components of the noise source, for the second
harmonic.

Next, we look at the second harmonic of the noise source. Here, group 2 of the
imaginary component is zero, and group 1 of the imaginary component comes from
delayed fission, as shown in Fig. 6. Most of the noise source is now in groups
1 and 2 of the real component. Again, we observe that the proposed exact noise
source sampling technique is in excellent agreement with the source calculated
by APOLLO3®. The noise source sampled using the approximate method agrees
quite well with the APOLLO3® source and the exact sampling in the imaginary
component, which is only due to delayed fission. Examining the real component,
it is again observed that the amplitude of the source using the approximate
method is smaller than the one obtained with the exact method. The asymmetry
which was observed in the first harmonic is less visible in the second
harmonic. In the first harmonic this asymmetry was visible exactly at the
material interface. For the second harmonic, the amplitude of the real
component of the first group is slightly smaller in the water than in the
fuel.

As the exact noise source sampling technique has now been verified, and
demonstrated to have far superior agreement with APOLLO3®, only the exact
noise source sampling method is used for the remainder of the paper.

###  VI.2 Application of Weight Cancellation

Next, the effects of weight cancellation on the performance of noise
simulations is examined. Three different cancellation methods were used:
approximate cancellation with a coarse mesh, approximate cancellation with a
fine mesh, and exact cancellation. A $170\times 170$ regular rectilinear mesh
was used for both the approximate coarse and exact cancellation methods. For
approximate cancellation with a fine mesh, a $340\times 340$ rectilinear mesh
was used instead. A simulation was also run without cancellation, and without
inner generations (essentially the algorithm used in TRIPOLI-4®), which served
as a baseline for all comparisons. The resulting noise field was scored using
the track-length estimator over a $138\times 138$ regular rectilinear mesh. To
compare the efficiency of each simulation method, we have used the average
ratio of the figure of merit (FOM) for results obtained using cancellation to
results obtained using the baseline algorithm. The FOM is defined as

$\text{FOM}=\frac{1}{T\sigma^{2}}\text{,}$ (46)

where $T$ is the wall-clock run time, and $\sigma^{2}$ is the variance of the
quantity being examined. Cancellation is only applied to the noise particles,
and can therefore only change the amount of time spent transporting noise
particles. It will have no effect on the time spent performing power iteration
to sample the noise particles to begin each noise replica. Because of this, it
can be argued that the FOM should only be calculated using the time spent
transporting noise particles. We have chosen to look at the FOM for both the
total run time (time spent during power iteration, noise transport, and
cancellation), and the noise run time (time spent during noise transport and
cancellation), as we think both quantities are of interest. For each run time,
the ratio for the FOM was computed in each element of the noise field scoring
mesh, and then the average of that ratio was calculated independently for each
component (real or imaginary) of the noise field and each energy group. As a
result, for each cancellation method, 8 FOM ratios are provided, 4 calculated
with the total run time, and 4 calculated with the noise run time.

####  VI.2.1 Analysis of the First Harmonic

Figure 7: First harmonic of the noise field running through the center of the
vibrating fuel pin, obtained with APOLLO3® and MGMC using different weight
cancellation methods (or no weight cancellation).

For the first harmonic of the original benchmark problem, a slice of the noise
field at the axis through the middle of the perturbed fuel pin is provided in
Fig. 7. These plots compare all four the Monte Carlo simulations to the
deterministic results obtained with APOLLO3®. Immediately, one is drawn to the
results from the baseline algorithm where no cancellation was used. Not only
are the error bars quite large for both components and both groups, but it
also appears as though there was a bias in the result, since this noise field
systematically has a higher value than in all other simulations. This is most
noticeable in the real component, and at regions far away from the perturbed
pin in the imaginary component. Running the same simulation parameters with a
different seed for the random number generation indicated that this behavior
is not systematic, and that sometimes a similar phenomenon is observed, but
with a systematically lower value than the deterministic results. This seems
to suggest that the Monte Carlo results from the baseline algorithm presented
here are far from converged, and would require a much larger number of noise
batches to generate better statistics. Unfortunately, this was not feasible
given the poor performance of this algorithm. It could also be the case that
the error bars are under-estimated because of the correlations induced by
power iteration; to test this hypothesis, further investigations with an
increased number of decorrelation cycles would be needed. Looking at the noise
field estimated with exact regional cancellation, there is certainly an
observable improvement in the resulting variance, and the systematic bias
observed without cancellation has been reduced. For the imaginary component,
the error bars are now overlapping the deterministic results at some
positions. Despite the disagreements observed far from the pin, very decent
agreement is observed in the vicinity of the perturbation, where the amplitude
of the noise field is largest.

The two approximate cancellation methods had the best agreement with the
deterministic results. While there were always visible differences in the real
component, there was always agreement within the error bars. For the imaginary
component, both approximate methods had nearly perfect agreement with
APOLLO3®, although the fine mesh variant did have slightly better agreement.
These small differences could be due to the approximate nature of the
cancellation technique, as refining the cancellation mesh should reduce the
bias in the results. Based on the error bars and the level of agreement
observed between the two approximate methods and the deterministic results, it
is difficult to conclude if the slight disagreement in the coarse mesh results
are actually due to the bias induced by approximate cancellation. At the
observed level of agreement, one must also call into question the accuracy of
the deterministic results, which have a bias induced by the geometric and
angular discretization.

|  | Real | Imaginary
---|---|---|---
|  | Group 1 | Group 2 | Group 1 | Group 2
Approximate Coarse | Total Run Time | 304 | 304 | 234 | 235
Noise Run Time | 1623 | 1624 | 1252 | 1254
Approximate Fine | Total Run Time | 129 | 129 | 122 | 122
Noise Run Time | 487 | 487 | 461 | 461
Exact | Total Run Time | 8 | 8 | 8 | 8
Noise Run Time | 14 | 14 | 14 | 14
Table 1: First harmonic improvement factors for the FOM of given weight
cancellation methods when compared to the baseline solution strategy which
does not use weight cancellation.

The improvements in the FOM for the simulations using cancellation, when
compared to the baseline algorithm without cancellation, are presented in Tab.
1. The approximate cancellation with the coarse mesh had the best improvement
in performance. Looking at the total run time, the real component was improved
by a factor of about 300 and the imaginary component by a factor of about 230.
If the noise run time is used, these become factors of 1620 and 1250
respectively. Approximate cancellation with the fine mesh gave an improvement
of approximately 120 when looking at the total run time, and approximately 480
when looking at the noise run time. The performance gains when using the fine
mesh might have been smaller due to less weight being cancelled with a finer
mesh, as a finer mesh leads to fewer noise particles per mesh element. This
performance penalty leads to a theoretical reduction in any bias in the
results, although we were unfortunately not able to measure this phenomenon.
Exact regional cancellation leads to the lowest FOM improvements, with a
factor of 8 observed when using the total run time, and a factor of 14 for the
noise run time.

A word of caution must be added to the FOM ratios provided in Tab. 1. All of
these ratios used the baseline algorithm results as a denominator. Since the
FOM only considers the variance of the examined score, these reported factors
do not consider the fact that the average values obtained with the baseline
algorithm seemed to have a systematic bias. There is no good way to
quantitatively measure the improvements in the average that cancellation
provided, as demonstrated in Fig. 7. As previously mentioned, this apparent
bias could also indicate that the error bars from the results obtained with
the baseline algorithm were underestimated, and the performance improvements
might actually be even larger; to support this statement, one would need to
perform independent replicas and estimate the ‘true’ error bars. Additionally,
the observed improvements in the FOM merit are not solely due to the
application of weight cancellation. It was observed that breaking the fixed-
source problem into inner generations, without applying any weight
cancellation between inner generations, reduced the simulation time when
compared to the baseline method. This phenomenon will be the subject of future
work.

####  VI.2.2 Analysis of the Second Harmonic

Figure 8: Second harmonic of the noise field running through the center of the vibrating fuel pin, obtained with MGMC using different approximate weight cancellation meshes, or without weight cancellation. |  | Real | Imaginary
---|---|---|---
|  | Group 1 | Group 2 | Group 1 | Group 2
Approximate Coarse | Total Run Time | 244 | 244 | 344 | 344
Noise Run Time | 1322 | 1323 | 1864 | 1864
Approximate Fine | Total Run Time | 119 | 119 | 152 | 151
Noise Run Time | 446 | 446 | 570 | 569
Table 2: Second harmonic improvement factors for the FOM of given weight
cancellation methods when compared to the baseline solution strategy which
does not use weight cancellation.

For the sake of completeness, we have then examined how our proposed approach
behaves when applied to the calculation of the second harmonic of the noise
field. We chose to use the two different approximate cancellation methods and
the baseline method to calculate the second harmonic, to evaluate the
performance gains provided by weight cancellation. It was chosen to not use
exact weight cancellation here, as its performance was significantly poorer
than the two approximate methods. The resulting noise fields are depicted in
Fig. 8. While the amplitude of the second harmonic in the vicinity of the
vibration appears to be smaller than the amplitude observed in the first
harmonic, this is not the case at positions farther from the perturbation. Far
to the left or right of the vibration, the amplitude of the second harmonic is
much larger than that of the first, which is coherent with previous findings
by Zoia et al. [16]. This might happen because of two possibly concurrent
reasons: specific symmetries of the system could suppress the first harmonic
and thus promote the second; moreover, it has been shown that the orthodox
linearization of the noise equations can nonphysically amplify the second
harmonic (with respect to the exact solution of the exact noise equations)
while leaving the first harmonic almost unaffected. Examining the real
component, it then appears as though the results from the baseline method
demonstrate a large bias, and do not agree with the two other results. The
estimated value of the imaginary component has much better agreement with the
values of the solutions which used approximate cancellation, but there is more
than $100\%$ relative error.

Focusing on the two approximate weight cancellation methods, both solutions
are in excellent agreement in the immediate vicinity of the vibrating pin in
the real component. Moving farther away from the pin, the results obtained
with the fine cancellation mesh have a slightly lower value than those
obtained with the coarse mesh. It is possible that this difference is due to
the finer mesh imposing less of a bias when compared to the coarse mesh. Both
solutions agree within one standard deviation however, and the results would
require further convergence to confirm this hypothesis. In both the real and
imaginary components, the structure of the static flux is visible, which was
not the case for the first harmonic results, which is coherent with the
expected behaviour of local and global components of the noise field (see e.g.
[3] and references therein).

Table 2 displays the FOM ratios for the two approximate cancellation methods,
when compared to the baseline method. In general, the performance improvements
are very similar to the values seen with the first harmonic. In the first
harmonic, the imaginary component of the noise field was dominant, while real
component is dominant for the second harmonic. Comparing the FOM ratio of the
coarse mesh for the imaginary component of the first harmonic with the real
component of the first harmonic, we see that the coarse mesh cancellation
method was actually slightly more effective at estimating the noise field for
the dominant component of the second harmonic than of the first. This is also
true for the non-dominant component with the coarse mesh as well. Strangely,
this is not true for the fine cancellation mesh, as this method had slightly
poorer performance on the dominant component of the second harmonic, than on
the first.

###  VI.3 The Effects of Vibration Frequency and Amplitude

Given the remarkable improvements in performance that weight cancellation has
enabled, it is now possible to envision the use of Monte Carlo noise
simulations to perform reactor analysis and assess the impact of specific
parameters on the noise field. In this section, examine the first harmonic for
two altered versions of the benchmark: one where the angular frequency of the
vibration of the fuel pin has been increased from
$\omega_{0}=2\pi\,$\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}$$ to
$\omega_{0}=4\pi\,$\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}$$ (without changing
the amplitude), and a second one where the vibration amplitude has been
increased from $\varepsilon=$0.2\text{\,}\mathrm{cm}$$ to
$\varepsilon=$0.4\text{\,}\mathrm{cm}$$ (without changing the frequency). The
resulting noise amplitude and phase for these two problems (in addition to
modulus and phase of the original benchmark parameters) are presented in Fig.
9. Approximate cancellation with the fine mesh was used to obtain these
results. In the previous sections, we examined the complex noise field
$\delta\varphi$ because this is the quantity which is estimated in the Monte
Carlo simulation, and the associated error bars are available, permitting a
comparison in performance with the FOM. Reactor physicists, however, are
typically more interested in the modulus ($\absolutevalue{\delta\varphi}$) and
phase
($\operatorname{atan2}\left(\Im{\delta\varphi},\Re{\delta\varphi}\right)$),
which can be calculated from the complex field. Despite having error bars for
the complex field, it was not possible to exhibit error bars for the modulus
and phase, as the covariance between the real and imaginary part was not
scored. Due to this limitation, no error bars are given in Fig. 9.

Figure 9: Modulus and phase of the noise field through the center of the
vibrating fuel pin, for different vibration parameters. All results were
obtained using MGMC, with approximate cancellation and the fine cancellation
mesh.

First, we will consider the effects of the proposed variations to the
vibration parameters on the noise modulus. From Fig. 9, we first note that in
the immediate vicinity of the perturbation, changing only the frequency seems
to have little impact on the modulus. Increasing the amplitude of the
vibration however, leads to an increase in magnitude and width of the peaks.
One would certainly expect a widening of the peak in this case, as the area of
the perturbation has been increased. With a larger perturbation region, more
noise source particles will be sampled in the vicinity of each interface,
contributing to the larger magnitude. Farther away from the perturbed pin, we
see that both the $\omega_{0}=4\pi\,$\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}$$
and the $\varepsilon=$0.4\text{\,}\mathrm{cm}$$ cases have a slightly larger
amplitude, although this effect is slightly less apparent in group 2. These
results could indicate that for this system the global component of the noise
is more sensitive to the frequency of the vibration, while the local component
is more sensitive to the amplitude of the vibration, in agreement with the
spectral analysis previously performed by Rouchon [54].

Examining the noise phase, it appears as though augmenting the perturbation
(either in frequency or amplitude) has little effect on the depth of the
observed well. Changing the vibration parameters did seem to change the width
of the phase well in both cases, and this effect is apparent in both energy
groups. The increased frequency reduced the well width, while the increased
amplitude enlarged it. Curiously, the observed changes in the well width only
ever occur on $+x$ side of the well in the first group, and on the $-x$ side
of the well for the second group. A physical explanation for these observed
phenomena would require further investigations.

##  VII Conclusions

In this work, we have considerably extended our previous findings concerning
the use of Monte Carlo methods for neutron noise simulations. The
contributions presented in this paper are twofold. First, we have proposed a
novel method to sample in an exact manner the noise source particles produced
from mechanical vibrations within a system. Contrary to the previous
implementation in TRIPOLI-4®, which for the sake of simplicity utilized many
approximations, we have shown that the noise source due to mechanical
vibrations can be represented without any approximation by mapping it to a
perturbation of the isotopic concentrations. For this purpose, a fictitious
material must be defined, containing the union of the nuclides on either side
of the vibrating interface. A nuclide is then randomly chosen from this
fictitious material, which allows sampling the noise source particles in a
manner which is nearly identical to the method used in the much simpler case
of cross section oscillations. This new method was implemented in the multi-
group Monte Carlo mini-app MGMC and verified against the noise source
calculated with the deterministic noise solver in APOLLO3®: the two noise
sources were in excellent agreement.

Second, we have applied weight cancellation methods as a variance reduction
technique for neutron noise simulations. For the case of mechanical
vibrations, the noise source contains nearly equal positive and negative
contributions in both the real and imaginary components, which makes it very
difficult to estimate the resulting noise field, because of very large
variances in scores. Preliminary investigations on the case of cross section
oscillations have shown that these issues can be overcome using weight
cancellation methods [36]. In view of probing the effects of weight
cancellation for the more involved case of mechanical vibrations, we have
implemented a new Monte Carlo noise solver in MGMC. We have started from the
Monte Carlo neutron noise method which was previously demonstrated in
TRIPOLI-4®, as it is able to produce error bars for the noise field. This
solution scheme was then modified, breaking the fixed-source noise batches
into inner fission generations, in order to accommodate weight cancellation.
To test our new solution strategy with weight cancellation, we have selected a
recent benchmark for neutron noise problems, consisting in a two-dimensional,
two-group reflected fuel assembly with a single vibrating fuel pin.
Approximate weight cancellation and exact regional weight cancellation were
tested on this benchmark problem: use of a relatively coarse mesh with
approximate cancellation leads to improvements in the FOM by factors as high
as 1624 when compared to the original neutron noise scheme implemented in
TRIPOLI-4®, at the expense of a slight but undetected bias in the noise field.
Use of a finer mesh for approximate cancellation, which is required in order
to quench the bias inherent to this method, yielded an improvement in the FOM
by factors as high as 487. By comparison, the exact weight cancellation
technique only improved the FOM by a factor of 14 and showed less potential
for implementation in production Monte Carlo codes, due to its algorithmic
complexity.

The application of approximate weight cancellation improves the quality of the
resulting noise field due to mechanical vibrations to the point that it
becomes feasible to use Monte Carlo methods to analyse changes in the noise
field stemming from small changes to the vibration frequency and amplitude;
previously, without the use of weight cancellation, such a comparison would
have been extremely challenging, if not impossible, for the case of a fuel pin
vibration [31]. Our results would indicate that the application of weight
cancellation to Monte Carlo neutron noise simulations is essentially mandatory
if one is to obtain usable results for the study of neutron noise in power
reactors. Although further investigations are required, our analysis shows
that the bias imposed on the resulting noise field by the approximate weight
cancellation is typically small, and more than acceptable since it can be
reduced by refining the cancellation mesh. While refining the mesh will reduce
the efficiency of cancellation (and also the efficiency of the noise
calculation), the observed improvements in the FOM for the examined benchmark
problem suggest that this is likely not a problem. In the future, we hope to
demonstrate this new Monte Carlo neutron noise solution scheme on larger
problems, such as vibrating fuel assemblies, or vibrating clusters of fuel
pins, to further demonstrate the potential for Monte Carlo method to be
applied to neutron noise analysis.

## Acknowledgements

TRIPOLI-4® and APOLLO3${}^{\text{\textregistered}}$ are registered trademarks
of CEA. The authors thank EDF and Framatome for partial financial support.

## References

  * [1] M. M. R. Williams, _Random Processes in Nuclear Reactors_. Pergamon Press, Oxford, 1974.
  * [2] J. A. Thie, _Power Reactor Noise_. American Nuclear Society, 1981.
  * [3] I. Pázsit and C. Demazière, “Noise techniques in nuclear systems,” in _Handbook of nuclear engineering_ , D. Cacuci, Ed. Springer Verlag, 2010, ch. 14, pp. 1629–1737.
  * [4] M. Seidl, K. Kosowski, U. Schüler, and L. Belblidia, “Review of the historic neutron noise behavior in German KWU built PWRs,” _Progress in Nuclear Energy_ , vol. 85, pp. 668–675, 2015.
  * [5] K. Behriguer, G. Kosály, and L. Kostic, “Reactivity and neutron density noise excited by random rod vibration,” _Annals of Nuclear Energy_ , vol. 63, pp. 306–318, 1977.
  * [6] D. N. Fry, J. March-Leuba, and F. J. Sweeney, “Use of neutron noise of diagnosis of in-vessel anomalies in light-water reactors,” Oak Ridge National Laboratory, Tech. Rep. ORNL/TM-8774, 1986.
  * [7] I. Pázsit and G. T. Analytis, “Theoretical investigation of the neutron noise diagnostics of two-dimensional control rod vibrations in a PWR,” _Annals of Nuclear Energy_ , vol. 7, pp. 171–183, 1980.
  * [8] G. Kosály, “Noise investigations in boiling-water and pressurized-water reactors,” _Progress in Nuclear Energy_ , vol. 5, pp. 145–199, 1980.
  * [9] I. Pázsit, L. A. Torres, M. Hursin, H. Nylén, V. Dykin, and C. Montalvo, “Development of a new method to determine the axial void velocity profile in BWRs from measurements of the in-core neutron noise,” _Progress in Nuclear Energy_ , vol. 138, p. 103805, 2021.
  * [10] V. Arzhanov and I. Pázsit, “Diagnostics of core barrel vibrations by in-core and ex-core neutron noise,” _Progress in Nuclear Energy_ , vol. 43, no. 1, pp. 151–158, 2003.
  * [11] I. Pázsit, “The linearization of vibration-induced noise,” _Annals of Nuclear Energy_ , vol. 11, pp. 441–454, 1984.
  * [12] I. Pázsit and J. Karlsson, “On the perturbative calculation of the vibration noise by strong absorber,” _Annals of Nuclear Energy_ , vol. 24, pp. 449–466, 1996.
  * [13] A. Jonsson, H. N. Tran, V. Dykin, and I. Pázsit, “Analytical investigation of the properties of the neutron noise induced by vibrationg absorber and fuel rods,” _Kerntechnik_ , vol. 77, pp. 371–380, 2012.
  * [14] T. Yamamoto and H. Sakamoto, “Decomposition of neutron noise in a reactor into higher-order mode components and investigation of the space and frequency dependence,” _Progress in Nuclear Energy_ , vol. 117, p. 103098, 2019.
  * [15] I. Pázsit and V. Dykin, “The role of the eigenvalue separation in reactor dynamics and neutron noise theory,” _J. Nucl. Sci. Technol._ , vol. 55, pp. 484–495, 2018.
  * [16] A. Zoia, A. Rouchon, B. Gasse, C. Demazière, and P. Vinai, “Analysis of the neutron noise induced by fuel assembly vibrations,” _Annals of Nuclear Energy_ , vol. 154, p. 108061, 2021.
  * [17] C. Demazière, A. Rouchon, and A. Zoia, “Understanding the neutron noise induced by fuel assembly vibrations in linear theory,” _Annals of Nuclear Energy_ , vol. 175, p. 109169, 2022.
  * [18] U. Rohde, M. Seidl, S. Kliem, and Y. Bilodid, “Neutron noise observations in German KWU built PWRs and analyses with the reactor dynamics code DYN3D,” _Annals of Nuclear Energy_ , vol. 112, pp. 715 – 734, 2018.
  * [19] C. Demazière, “CORESIM: a multi-purpose neutronic tool for research and education,” _Annals of Nuclear Energy_ , vol. 38, p. 2698–2718, 2011.
  * [20] T. Yamamoto, “Monte carlo method with complex-valued weights for frequency domain analyses of neutron noise,” _Annals of Nuclear Energy_ , vol. 58, p. 72–79, 2013.
  * [21] A. Rouchon, A. Zoia, and R. Sanchez, “A new Monte Carlo method for neutron noise calculations in the frequency domain,” _Annals of Nuclear Energy_ , vol. 102, p. 465–475, 2017.
  * [22] N. Olmo-Juan, C. Demazière, T. Barrachina, R. Miró, and G. Verdú, “PARCS vs CORE SIM neutron noise simulations,” _Progress in Nuclear Energy_ , vol. 115, pp. 169 – 180, 2019.
  * [23] A. Vidal-Ferrándiz, A. Carreno, D. Ginestar, C. Demazière, and G. Verdú, “A time and frequency domain analysis of the effect of vibrating fuel assemblies on the neutron noise,” _Annals of Nuclear Energy_ , vol. 137, p. 107076, 2020.
  * [24] D. Chionis, A. Dokhane, L. Belblidia, H. Ferroukhi, G. Girardin, and A. Pautz, “Development and verification of a methodology for neutron noise response to fuel assembly vibrations,” _Annals of Nuclear Energy_ , vol. 147, p. 107669, 2020.
  * [25] C. Demazière, P. Vinai, M. Hursin, S. Kollias, and J. Herb, “Overview of the CORTEX project,” in _Proceedings of Physor2018_ , Cancun, Mexico, 2018\.
  * [26] H. Yi, P. Vinai, and C. Demazière, “On the simulation of neutron noise using a discrete ordinates method,” _Annals of Nuclear Energy_ , vol. 164, p. 108570, 2021.
  * [27] T. Yamamoto, “Implementation of a frequency-domain neutron noise analysis method in a production-level continuous energy Monte Carlo code: verification and application to a BWR,” _Annals of Nuclear Energy_ , vol. 115, pp. 494–501, 2018.
  * [28] A. Rouchon, W. Jarrah, and A. Zoia, “The new neutron noise solver of the Monte Carlo code TRIPOLI-4,” in _Proceedings of M &C2019_, Portland USA, 2019.
  * [29] A. Mylonakis, H. Yi, P. Vinai, and C. Demazière, “Neutron noise simulations in a heterogeneous system: A comparison between a diffusion-based and a discrete ordinates solver,” in _Proceedings of M &C2019_, Portland USA, 2019.
  * [30] P. Vinai, H. Yi, A. Mylonakis, C. Demazière, B. Gasse, A. Rouchon, A. Zoia, A. Vidal-Ferràndiz, D. Ginestar, G. Verdú, and T. Yamamoto, “Comparison of Neutron Noise Solvers Based on Numerical Benchmarks in a 2-D Simplified UOX Fuel Assembly,” in _M &C 2021_, 2021.
  * [31] P. Vinai, H. Yi, C. Demazière, A. Rouchon, A. Zoia, A. Vidal-Ferràndiz, A. Carreño, D. Ginestar, and G. Verdú, “On the simulation of neutron noise induced by vibrations of fuel pins in a simplified fuel assembly,” _Annals of Nuclear Energy_ , Submitted.
  * [32] Lamirand, Vincent _et al._ , “Neutron noise experiments in the akr-2 and crocus reactors for the european project cortex,” _EPJ Web Conf._ , vol. 225, p. 04023, 2020.
  * [33] Mylonakis, A. G. _et al._ , “CORE SIM+ simulations of COLIBRI fuel rods oscillation experiments and comparison with measurements,” _EPJ Web Conf._ , vol. 247, p. 21006, 2021.
  * [34] H. Belanger, D. Mancusi, and A. Zoia, “Exact weight cancellation in Monte Carlo eigenvalue transport problems,” _Physical Review E_ , vol. 104, no. 1, p. 015306, 2021.
  * [35] ——, “Solving Eigenvalue Transport Problems with Negative Weights and Regional Cancellation,” in _M &C 2021_, Oct 2021.
  * [36] ——, “Variance Reduction Techniques for Monte Carlo Neutron Noise Simulations,” in _PHYSOR 2022_ , 2022.
  * [37] I. Pázsit, “Investigation of the space-dependent noise induced by a vibrating absorber,” _Atomkernenergie_ , vol. 30, pp. 29–35, 1977.
  * [38] M. Faucher, D. Mancusi, and A. Zoia, “New kinetic simulation capabilities for Tripoli-4®: Methods and applications,” _Annals of Nuclear Energy_ , vol. 120, p. 74–88, 2018.
  * [39] E. Brun, F. Damian, C. Diop, E. Dumonteil, F. Hugot, C. Jouanne, Y. Lee, F. Malvagi, A. Mazzolo, O. Petit, J. Trama, T. Visonneau, and A. Zoia, “TRIPOLI-4®, CEA, EDF and AREVA reference Monte Carlo code,” _Annals of Nuclear Energy_ , vol. 82, p. 151–160, 2015.
  * [40] L. L. Carter, E. D. Cashwell, and W. M. Taylor, “Monte Carlo Sampling with Continuously Varying Cross Sections Along Flight Paths,” _Nuclear Science and Engineering_ , vol. 48, no. 4, p. 403–411, 1972.
  * [41] D. Legrady, B. Molnar, M. Klausz, and T. Major, “Woodcock tracking with arbitrary sampling cross section using negative weights,” _Annals of Nuclear Energy_ , vol. 102, p. 116–123, 2017.
  * [42] T. E. Booth, “Computing the Higher k-Eigenfunctions by Monte Carlo Power Iteration: A Conjecture,” _Nuclear Science and Engineering_ , vol. 143, no. 3, p. 291–300, 2003.
  * [43] T. Yamamoto, “Convergence of the second eigenfunction in Monte Carlo power iteration,” _Annals of Nuclear Energy_ , vol. 36, no. 1, p. 7–14, 2009\.
  * [44] ——, “Monte Carlo algorithm for buckling search and neutron leakage-corrected calculations,” _Annals of Nuclear Energy_ , vol. 47, p. 14–20, 2012.
  * [45] ——, “Monte Carlo method with complex weights for neutron leakage-corrected calculations and anisotropic diffusion coefficient generations,” _Annals of Nuclear Energy_ , vol. 50, p. 141–149, 2012.
  * [46] ——, “Implementation of a frequency-domain neutron noise analysis method in a production-level continuous energy Monte Carlo code: Verification and application in a BWR,” _Annals of Nuclear Energy_ , vol. 115, p. 494–501, 2018.
  * [47] P. Zhang, H. Lee, and D. Lee, “A general solution strategy of modified power method for higher mode solutions,” _Journal of Computational Physics_ , vol. 305, p. 387–402, 2016.
  * [48] T. E. Booth and J. E. Gubernatis, “Exact Regional Monte Carlo Weight Cancellation for Second Eigenfunction Calculations,” _Nuclear Science and Engineering_ , vol. 165, no. 3, p. 283–291, 2010.
  * [49] E. R. Woodcock, T. Murphy, P. J. Hemmings, and T. C. Longworth, “Techniques used in the GEM code for Monte Carlo neutronics calculations in reactors and other systems of complex geometry,” Argonne National Laboratory, Tech. Rep., 1965, ANL-7050.
  * [50] J. Leppänen, “On the use of delta-tracking and the collision flux estimator in the Serpent 2 Monte Carlo particle transport code,” _Annals of Nuclear Energy_ , vol. 105, p. 161–167, 2017.
  * [51] H. Belanger, D. Mancusi, A. Rouchon, and A. Zoia, “Unbiasedness and Optimization of Regional Weight Cancellation,” _Physical Review E_ , Submitted.
  * [52] H. Belanger, “MGMC,” 2022, v0.2.0. [Online]. Available: https://github.com/HunterBelanger/mgmc
  * [53] C. R. Harris _et al._ , “Array programming with NumPy,” _Nature_ , vol. 585, no. 7825, p. 357–362, 2020.
  * [54] A. Rouchon, “Analyse et développement d’outils numériques déterministes et stochastiques résolvant les équations du bruit neutronique et applications aux réacteurs thermiques et rapides,” Ph.D. dissertation, Université Paris-Saclay, 2016.
  * [55] D. Schneider, F. Dolci, F. Gabriel, J. Palau, M. Guillo, B. Pothet, P. Archier, K. Ammar, F. Auffret, R. Baron _et al._ , “APOLLO3: CEA/DEN deterministic multi-purpose code for reactor physics analysis,” in _PHYSOR 2016_ , 2016.
  * [56] A. Rouchon, M. V. Le Brun, and A. Zoia, “Analysis and Comparison of APOLLO3 and TRIPOLI-4 Neutron Noise Solvers,” in _PHYSOR 2020_ , vol. 247, 2021, p. 21002. [Online]. Available: https://doi.org/10.1051/epjconf/202124721002