# Limit theorems for the neutron transport equation

Eric Dumonteil, Emma Horton , Andreas E. Kyprianou† and Andrea Zoia CEA, Paris
Saclay. E-mail<EMAIL_ADDRESS>Department of Statistics, University of
Warwick, Coventry, CV4 7AL, UK. E-mail: {emma.horton},
<EMAIL_ADDRESS>Paris Saclay. E-mail<EMAIL_ADDRESS>
###### Abstract

Over the last decade, ingenuous developments in Monte Carlo methods have
enabled the unbiased estimation of adjoint-weighted reactor parameters
expressed as bilinear forms, such as kinetics parameters and sensitivity
coefficients. A prominent example is the Iterated Fission Probability method,
which relies on the simulation of the fission chains descending from an
ancestor neutron: the neutron population at an asymptotic fission generation
yields an estimate of the importance function (and hence of the adjoint
fundamental eigenmode) at the phase-space coordinates of the ancestor neutron.
In this paper we first establish rigorous results concerning the moments of
the asymptotic neutron population stemming from a single initial particle,
with special focus on the average and the variance. Then, we propose a simple
benchmark configuration where exact solutions are derived for these moments,
which can be used for the verification of new functionalities of production
Monte Carlo codes involving the Iterated Fission Probability method.

Key words: neutron transport, k-effective, moments, branching processes.

MSC 2020: 60J80, 60J85, 82D75, 60J05

##  1 Introduction

Monte Carlo simulation is the gold standard method for radiation transport
applications, since it involves a minimal number of approximations. In
particular, the phase space does not need to be discretized and all the
particle-nuclei interaction physics contained in the nuclear data libraries
can be used [19]. Thanks to these features, particle-transport codes based on
the Monte Carlo method allow one to establish reference solutions against
which those produced by faster but approximate deterministic solvers (which
rely on the discretization of the phase space) can be benchmarked for
accuracy.

The stochastic version of the standard power iteration algorithm is the
workhorse of Monte Carlo codes for $k$-eigenvalue problems: a collection of
neutrons is followed over a sufficiently large number of fission generations,
and statistics are recorded on the successive generations once the population
has settled into a stationary state. At equilibrium, the ratio between the
statistical weights of the neutrons at two consecutive fission generations
converges to the fundamental eigenvalue $k_{0}$, and the neutron flux within a
generation correspondingly converges to the fundamental $k$-eigenmode,
$\varphi_{0}$ [19]. Over the last decade, a major breakthrough based on the
rediscovery of the Iterated Fission Probability (IFP) method has enabled the
use of Monte Carlo to compute adjoint-weighted quantities [21, 13]. The key
idea is that the adjoint fundamental eigenmode of a $k$-eigenvalue calculation
is proportional to the neutron importance, which can be estimated by tracking
the descendants of an ancestor neutron over a sufficiently large number of
power iteration generations and collecting the statistical weights of the
surviving neutrons. This technique has been implemented in several production
Monte Carlo codes, enabling unbiased estimates for the fundamental adjoint
flux [28, 26] or adjoint-weighted parameters (bilinear forms) such as kinetics
parameters or perturbations and sensitivities to nuclear data [30, 12, 27, 2,
23, 14].

When developing new algorithms and functionalities in production Monte Carlo
codes such as MCNP6® [18], MONK® [24] or TRIPOLI-4® [1], which involve several
hundred thousands of lines of code, it is of utmost importance to rely on
exact solutions whenever possible for verification purposes. Several such sets
of analytical solutions have been established for regular (forward) transport
problems, and in particular $k$-eigenvalue calculations [7, 25, 17]. For
adjoint eigenvalue problems and related adjoint-weighted parameters,
comparatively fewer exact solutions have been derived. A prominent example is
the two-group, infinite medium model benchmark [15], which has been
successfully used to verify recent developments of Monte Carlo codes [30, 13,
2]. In view of these considerations, in this work we set out to establish
benchmark solutions for adjoint $k$-eigenvalue problems in spatially finite
media, which can usefully complement those of the infinite medium model.

Since IFP methods are generally concerned with the asymptotic behaviour of a
stochastic neutron population, in this paper we provide general statements
about the moments of the neutron population as a function of fission
generations. More precisely, we provide exact asymptotics for the moments (of
any order) of the fission generation populations for the sub-, super- and
critical cases. In particular, this implies that we can obtain precise
asymptotics for the average and variance of the neutron population. We further
introduce a benchmark configuration where exact results can be derived. For
this purpose, we resort to the ‘rod-model’, a simple transport configuration
where neutrons are constrained to move along a line, the only permissible
directions being those in the increasing or decreasing spatial coordinate
[31]. Within this framework, we establish reference solutions for the
asymptotic average number of neutrons stemming from an ancestor particle, as
well as for the second moment of this counting process, for sub-critical,
critical and super-critical configurations. These findings can be used as an
ideal verification test-bed for Monte Carlo code developers interested in IFP-
based algorithms.

This paper is organized as follows. In Sec. 2, we introduce the general
neutron transport model that we will work with and present a stochastic
representation of the model in the form of a branching process. In Sec. 3, we
give an informal statement of the main results. In Sec. 4 we consider the
specific example of the rod model to illustrate the agreement between our
theoretical results and Monte Carlo simulations in some simple cases.
Conclusions will finally be drawn in Sec. 5. The precise statements of the
main results described in Sec. 3 are provided in the Appendix (Sec. A) with
the proofs being given in the supplementary material of Ref. 5.

##  2 The stochastic interpretation of the neutron transport equation

Consider a collection of neutrons evolving in the position (${\bf r}$),
direction (${\mathbf{\Omega}}$) and energy ($E$) phase space
$\mathcal{S}:=D\times\mathbb{S}_{2}\times(E_{\mathtt{min}},E_{\mathtt{max}})$,
where the spatial domain $D\subset\mathbb{R}^{3}$ is open and bounded, the
direction domain is the unit sphere $\mathbb{S}_{2}\subset\mathbb{R}^{3}$, and
the energy domain is $(E_{\mathtt{min}},E_{\mathtt{max}})$, where
$0<E_{\mathtt{min}}\leq E_{\mathtt{max}}<\infty$.

It is customarily assumed that neutron trajectories behave stochastically,
their randomness being due to the interactions with the surrounding nuclei;
furthermore, thanks to the very low density of neutrons in matter, their
transport process is inherently linear, in that the probability of neutron-
neutron interactions is negligible for all practical purposes. Between
collisions, particles move along straight lines whose length obeys a non-
homogeneous exponential probability distribution with parameter
$\Sigma_{\mathtt{t}}({\bf r},E)$, with units given by inverse length. The
total macroscopic cross section $\Sigma_{\mathtt{t}}({\bf r},E)$ defines the
probability per unit length that the neutron has a collision in the following
infinitesimal displacement about its current position ${\bf r}$, in the
direction ${\mathbf{\Omega}}$ and with energy $E$. At the scale of neutron
paths, materials are ideally isotropic, so that $\Sigma_{\mathtt{t}}$ does not
depend on ${\mathbf{\Omega}}$. Displacements are formally associated to the
backward streaming operator

(1) $T^{\dagger}[g]({\bf r},{\mathbf{\Omega}},E)=-{\mathbf{\Omega}}\cdot\nabla
g({\bf r},{\mathbf{\Omega}},E)+\Sigma_{\mathtt{t}}({\bf r},E)g({\bf
r},{\mathbf{\Omega}},E),\quad({\bf r},{\mathbf{\Omega}},E)\in\mathcal{S},$

with $g\in L^{\infty}_{+}(\mathcal{S})$ the collection of non-negative,
uniformly bounded, measurable functions on $\mathcal{S}$.

We assume that the domain $D$ has leakage boundary conditions, which means
that neutrons leaving the domain are lost. At the end of each flight, the
neutron undergoes a collision event, provided it is still within the viable
domain. We have
$\Sigma_{\mathtt{t}}=\Sigma_{\mathtt{c}}+\Sigma_{\mathtt{s}}+\Sigma_{\mathtt{f}}$,
where $\Sigma_{\mathtt{c}}$ (resp. $\Sigma_{\mathtt{s}}$,
$\Sigma_{\mathtt{f}}$) denotes the capture (resp. scattering, fission) cross
section. The ratios $\Sigma_{\mathtt{r}}/\Sigma_{\mathtt{t}}$, with
$\mathtt{r}=\mathtt{c}$ (resp. $\mathtt{s}$, $\mathtt{f}$), yield the
probability that the collision event is capture (resp. scattering or fission).
If the neutron is captured, its history is terminated; if it undergoes
scattering, its direction and energy coordinates are randomly distributed
according to a (normalized) probability density $f_{\mathtt{s}}({\bf
r},{\mathbf{\Omega}}\to{\mathbf{\Omega}}^{\prime},E\to E^{\prime})$; if it
undergoes fission, it disappears and is replaced by a random number of new
neutrons, whose average number is $\nu_{\mathtt{f}}({\bf r},E)$ and whose
direction and energy coordinates are randomly distributed according to a
(normalized) probability density $\chi_{\mathtt{f}}({\bf
r},{\mathbf{\Omega}}\to{\mathbf{\Omega}}^{\prime},E\to E^{\prime})$.
Scattering and fission are formally associated to the backward operators

(2) $\displaystyle S^{\dagger}[g]({\bf r},{\mathbf{\Omega}},E)$
$\displaystyle=\Sigma_{\mathtt{s}}({\bf
r},E)\int_{E_{\mathtt{min}}}^{E_{\mathtt{max}}}\int_{\mathbb{S}_{2}}g({\bf
r},{\mathbf{\Omega}}^{\prime},E^{\prime})f_{\mathtt{s}}({\bf
r},{\mathbf{\Omega}}\to{\mathbf{\Omega}}^{\prime},E\to E^{\prime}){\rm
d}{\mathbf{\Omega}}^{\prime}{\rm d}E^{\prime},$ (3) $\displaystyle
F^{\dagger}[g]({\bf r},{\mathbf{\Omega}},E)$
$\displaystyle=\nu_{\mathtt{f}}({\bf r},E)\Sigma_{\mathtt{f}}({\bf
r},E)\int_{E_{\mathtt{min}}}^{E_{\mathtt{max}}}\int_{\mathbb{S}_{2}}g({\bf
r},{\mathbf{\Omega}}^{\prime},E^{\prime})\chi_{\mathtt{f}}({\bf
r},{\mathbf{\Omega}}\to{\mathbf{\Omega}}^{\prime},E\to E^{\prime}){\rm
d}{\mathbf{\Omega}}^{\prime}{\rm d}E^{\prime}.$

For now, in order to keep the notation simple, we only focus on prompt
neutrons however, we refer the reader to Remark 1 for an explanation of how
this is easily extended to the case where delayed neutrons are also accounted
for.

The backward operators $T^{\dagger}$, $S^{\dagger}$ and $F^{\dagger}$ are
related to each other via the $k$-eigenvalue formulation of the adjoint
neutron transport equation (NTE) [19], which reads

(4)
$(T^{\dagger}-S^{\dagger})[\varphi^{\dagger}]=\frac{1}{k}F^{\dagger}[\varphi^{\dagger}]\quad\Leftrightarrow\quad(T^{\dagger}-S^{\dagger})^{-1}F^{\dagger}[\varphi^{\dagger}]=k\varphi^{\dagger}.$

Taking the adjoint of all the operators in Eq. (4), we have the customary
forward formulation

(5)
$(T-S)[\varphi]=\frac{1}{k}F[\varphi]\quad\Leftrightarrow\quad(T-S)^{-1}F[\varphi]=k\varphi,$

where $T$, $S$ and $F$ are the formal adjoints of the operators $T^{\dagger}$,
$S^{\dagger}$ and $F^{\dagger}$, respectively. The fundamental eigenmodes
$\varphi^{\dagger}_{0}({\bf r},{\mathbf{\Omega}},E)$ and $\varphi_{0}({\bf
r},{\mathbf{\Omega}},E)$ of Eqs. (4) and (5), respectively, can be given a
physical meaning: $\varphi^{\dagger}_{0}({\bf r},{\mathbf{\Omega}},E)$ is the
importance of a neutron injected into the system with phase-space coordinates
$({\bf r},{\mathbf{\Omega}},E)$, and $Q_{\mathtt{f}}=F[\varphi_{0}]$ is the
stationary post-fission neutron distribution (with $\varphi_{0}$ being the
pre-fission steady state). We will expand on this further in due course. Since
$\varphi^{\dagger}_{0}$ and $Q_{\mathtt{f}}$ are unique up to multiplicative
constants, we specify their normalisation as

(6) $\langle Q_{\mathtt{f}},\mathbf{1}\rangle=1\quad\text{ and }\quad\langle
Q_{\mathtt{f}},\varphi^{\dagger}_{0}\rangle=1,$

where $\mathbf{1}$ is the function that takes the value $1$ everywhere and the
angle brackets denote the scalar product. The physical picture related to Eqs.
(4) and (5) is the following: given a population of pre-fission particles
distributed according to $\varphi_{0}$, $F$ initiates the next generation of
fission particles and $(T-S)^{-1}$ transports them to the next set of
absorption (capture plus fission) sites, resulting in another population of
particles distributed according to $\varphi_{0}$, albeit multiplied by a
factor of $k_{0}$. In this respect, it is natural to introduce the concept of
fission generation, and rewrite Eq. (5) as

(7) $F(T-S)^{-1}Q_{\mathtt{f}}=k_{0}Q_{\mathtt{f}},$

where the fundamental eigenvalue $k_{0}$ physically represents the ratio
between the number of neutrons in two successive fission generations. If
$k_{0}>1$ the system is supercritical, corresponding to an exponential growth
of the number of particles in system with respect to generations; if $k_{0}<1$
the system is subcritical, corresponding to an exponential decay in the number
of particles; and if $k_{0}=1$ there is an equilibrium between neutron loss
(absorption and leakage) and gain (fission).

The probabilistic interpretation of Eqs. (4) and (5) is now made explicit. Let
$N_{n}$ denote the number of neutrons that stem from the $n$-th fission event
in their genealogical line of descent, with $n\geq 1$; furthermore, let
$\\{({\bf
r}_{i}^{(n)},{\mathbf{\Omega}}_{i}^{(n)},E_{i}^{(n)}),\,i=1,\dots,N_{n}\\}$
denote their phase space configurations. Our goal is to characterize the
statistical behaviour of the collection of neutrons in the $n$-th fission
generation, which, in view of the previous remarks, is clearly a discrete-time
branching process. Define

(8) $\mathcal{X}_{n}(A):=\sum_{i=1}^{N_{n}}\delta_{({\bf
r}_{i}^{(n)},{\mathbf{\Omega}}_{i}^{(n)},E_{i}^{(n)})}(A),\qquad
A\in\mathcal{B}(\mathcal{S}),n\geq 1,$

where $\mathcal{B}(\mathcal{S})$ is the set of Borel subsets of $\mathcal{S}$.
Then $\mathcal{X}_{n}(A)$ denotes the number of post-fission neutrons in
$A\subset\mathcal{S}$. Moreover, $\mathcal{X}=(\mathcal{X}_{n})_{n\geq 1}$ is
a branching process with expectation semigroup

(9) $\Psi_{n}[f]({\bf r},{\mathbf{\Omega}},E):=\mathbb{E}_{\delta_{({\bf
r},{\mathbf{\Omega}},E)}}\left[\sum_{i=1}^{N_{n}}f({\bf
r}_{i}^{(n)},{\mathbf{\Omega}}_{i}^{(n)},E_{i}^{(n)})\right],$

where $\mathbb{E}_{\delta_{({\bf r},{\mathbf{\Omega}},E)}}$ is the expectation
operator associated with the law $\mathbb{P}_{\delta_{({\bf
r},{\mathbf{\Omega}},E)}}$ of the branching process when initiated from a
single particle at $({\bf r},{\mathbf{\Omega}},E)\in\mathcal{S}$ and $f\in
L_{+}^{\infty}(\mathcal{S})$, the set of non-negative, bounded measurable
functions on $\mathcal{S}$. The physical interpretation of $\Psi_{n}[f]({\bf
r},{\bf\Omega},E)$ is the expected behaviour of $\mathcal{X}_{n}$ when
initiated from a single neutron with configuration $({\bf
r},{\mathbf{\Omega}},E)\in\mathcal{S}$. For example, taking $f=1$, then
$\Psi_{n}[f]({\bf r},{\bf\Omega},E)$ is equal to $\mathbb{E}_{\delta_{({\bf
r},{\bf\Omega},E)}}[N_{n}]$, the expected number of $n$-th generation post-
fission neutrons, which has direct bearing on the IFP algorithm.

##  3 Statistics of fission chains: main results

It is well known, see e.g. Ref. 3, that $\Psi_{n}[f]$ is related to $k_{0}$,
$\varphi_{0}^{\dagger}$ and $Q_{\mathtt{f}}$ by

(10)
$\Psi_{n}[\varphi_{0}^{\dagger}]=k_{0}^{n}\varphi_{0}^{\dagger},\quad\text{
and }\quad\langle Q_{\mathtt{f}},\Psi_{n}[f]\rangle=k_{0}^{n}\langle
Q_{\mathtt{f}},f\rangle,$

where brackets denote again the scalar product. Furthermore, it is also known
that, see also Ref. 3,

(11) $\Psi_{n}[f]({\bf r},{\mathbf{\Omega}},E)\sim k_{0}^{n}\langle
Q_{\mathtt{f}},f\rangle\varphi_{0}^{\dagger}({\bf r},{\mathbf{\Omega}},E),$

for large $n$. In particular, taking $f=1$ in the above asymptotic and using
Eq. (6), we obtain

(12) $\mathbb{E}_{\delta_{({\bf r,\bf\Omega,E})}}[N_{n}]\sim
k_{0}^{n}\varphi_{0}^{\dagger}({\bf r},{\bf\Omega},E),$

which shows that the expected number of particles in the system for large $n$
is proportional to $\varphi_{0}^{\dagger}({\bf r},{\bf\Omega},E)$,
demonstrating that $\varphi_{0}^{\dagger}$ indeed physically represents the
importance of an initial particle at $({\bf r},{\bf\Omega},E)$ to the growth
of the neutron population in the system.

In addition to the asymptotic behaviour of the first moment given in Eq. (11),
in the last decade several results pertaining to the asymptotic behaviour of
moments of functionals of $(\mathcal{X}_{n})_{n\geq 1}$ have emerged, see for
example Refs. 33, 34, 35 and Refs. 9, 6 for general branching processes. In
the present article, we unify and extend the aforementioned results in the
spirit of Ref. 9, which deals with the evolution of a neutron population
described in terms of a branching process in a continuous-time setting. More
precisely, for $n,\ell\geq 1,f\in L^{\infty}_{+}(\mathcal{S})$ and $({\bf
r},{\mathbf{\Omega}},E)\in\mathcal{S}$, define

(13) $\Psi_{n}^{(\ell)}[f]({\bf
r},{\mathbf{\Omega}},E):=\mathbb{E}_{\delta_{({\bf
r},{\mathbf{\Omega}},E)}}\left[\left(\sum_{i=1}^{N_{n}}f({\bf
r}_{i}^{(n)},{\mathbf{\Omega}}_{i}^{(n)},E_{i}^{(n)})\right)^{\ell}\right].$

We will show that, regardless of the value of $k_{0}$, for each $\ell\geq 1$
there exist $g_{\ell}(n)$ and $h_{\ell}$ such that

(14) $\Psi_{n}^{(\ell)}[f]({\bf
r},{\mathbf{\Omega}},E)\sim\frac{h_{\ell}}{g_{\ell}(n)}\varphi_{0}^{\dagger}({\bf
r},{\mathbf{\Omega}},E)$

for large $n$. Note that, thanks to Eq. (11), in the case $\ell=1$ we have
$g_{1}(n)=k_{0}^{-n}$ and $h_{1}=\langle Q_{\mathtt{f}},f\rangle$.

Our second contribution will be to show that, in the critical case when
$k_{0}=1$, one can obtain a precise asymptotic for the survival probability,
$\mathbb{P}_{\delta_{({\bf r},{\mathbf{\Omega}},E)}}(N_{n}>0)$ for a
population descending from a single particle starting at $({\bf
r},{\mathbf{\Omega}},E)$. This result, combined with the moment asymptotics in
the critical case, implies that the limiting distribution of
$\frac{1}{n}\sum_{i=1}^{N_{n}}f({\bf
r}_{i}^{(n)},{\mathbf{\Omega}}_{i}^{(n)},E_{i}^{(n)})$, conditional on
$N_{n}>0$, is exponential with a rate that can be explicitly determined by the
model parameters. The latter two results are classical in the branching
processes literature: the survival probability dates back to Kolmogorov, see
Ref. 16, and the second result is due to Yaglom, see Ref. 32. We also refer
the reader to Refs. 10, 22 for continuous-time results, and to Refs. 8, 20 for
discrete-time results. The proofs of our derivations are provided in Sec. A
and in the supplementary material given in Ref. 5. In the case where one
assumes reflective boundary conditions instead of leakage, the moment
asymptotic results (for all values of $k_{0}$) would still hold as long as Eq.
(11) is true: this is due to the fact that the proofs are inductive, so only
require the result to be true for the first moment.

###  3.1 The critical case

For the critical case, with $k_{0}=1$, it is possible to show (see Sec. A)
that

(15) $\Psi_{n}^{(\ell)}[f]({\bf r},{\mathbf{\Omega}},E)\sim
n^{\ell-1}\frac{\ell!}{2^{\ell-1}}\,\langle
Q_{\mathtt{f}},f\rangle^{\ell}\langle\varphi_{0},\Sigma_{\mathtt{f}}\mathcal{V}[\varphi_{0}^{\dagger}]\rangle^{\ell-1}\varphi_{0}^{\dagger}({\bf
r},{\mathbf{\Omega}},E)$

for sufficiently large $n$, where for $g\in L^{\infty}_{+}(\mathcal{S})$ we
have defined

(16) $\mathcal{V}[g]({\bf r},{\mathbf{\Omega}},E)=\mathcal{E}_{({\bf
r},{\mathbf{\Omega}},E)}\left[\left(\sum_{i=1}^{M}g({\bf
r},{\mathbf{\Omega}}_{i},E_{i})\right)^{2}\right]-\mathcal{E}_{({\bf
r},{\mathbf{\Omega}},E)}\left[\sum_{i=1}^{M}g^{2}({\bf
r},{\mathbf{\Omega}}_{i},E_{i})\right],$

where $M$ denotes the random number of neutrons produced from a fission event
and $\mathcal{E}_{({\bf r},{\mathbf{\Omega}},E)}$ is the average over the
fission offspring (number and configurations) produced from a fission event at
$({\bf r},{\mathbf{\Omega}},E)$.

To illustrate Eq. (15) further, we consider some simples cases. For example,
when $\ell=1$, we see that we recover Eq. (11) and hence, in particular,
taking $f=1$, we recover Eq. (12). Now let us consider the case when $f=1$ and
$\ell=2$. In this case, thanks to Eq. (6), Eq. (15) becomes

(17) $\mathbb{E}_{\delta_{({\bf r},{\mathbf{\Omega}},E)}}[N_{n}^{2}]\sim
n\langle\varphi_{0},\Sigma_{\mathtt{f}}\mathcal{V}[\varphi_{0}^{\dagger}]\rangle\varphi_{0}^{\dagger}({\bf
r},{\mathbf{\Omega}},E),$

for large $n$, which means that the second moment diverges linearly with
respect to the fission generations.

As promised, for the critical case, we also give a precise asymptotic for the
survival probability. Theorem 2 of Sec. A shows that for any $({\bf
r},{\mathbf{\Omega}},E)\in\mathcal{S}$ we have

(18) $\mathbb{P}_{\delta_{({\bf
r},{\mathbf{\Omega}},E)}}(N_{n}>0)\sim\frac{2}{n}\frac{\varphi_{0}^{\dagger}({\bf
r},{\mathbf{\Omega}},E)}{\langle\varphi_{0},\Sigma_{\mathtt{f}}\mathcal{V}[\varphi_{0}^{\dagger}]\rangle}$

as $n\to\infty$, where $\mathcal{V}$ was defined in Eq. (16). Here we see
another perspective of the aforementioned interpretation of
$\varphi_{0}^{\dagger}$ in the critical case, in that the survival probability
is proportional to $\varphi_{0}^{\dagger}$.

Combining Eqs. (15) and (18), it is straightforward to show that for any
$({\bf r},{\mathbf{\Omega}},E)\in\mathcal{S}$ and any $f\in
L^{\infty}_{+}(\mathcal{S})$, conditional on survival, the distribution of the
process normalised by the current generation $n$ is asymptotically
exponential. That is,

(19) $\left(\frac{1}{n}\sum_{i=1}^{N_{n}}f({\bf
r}_{i}^{(n)},{\mathbf{\Omega}}_{i}^{(n)},E_{i}^{(n)})\Bigg{|}N_{n}>0\right)\sim
Y$

for large $n$, where $Y$ is an exponential random variable with rate
$2/\langle
Q_{\mathtt{f}},f\rangle\langle\varphi_{0},\Sigma_{\mathtt{f}}\mathcal{V}[\varphi_{0}^{\dagger}]\rangle$.
We refer the reader to Ref. [5] for the details of the proof.

###  3.2 The supercritical case

For the supercritical case, with $k_{0}>1$, the findings of Sec. A show that

(20) $\Psi_{n}^{(\ell)}[f]({\bf r},{\mathbf{\Omega}},E)\sim
k_{0}^{n\ell}\,\langle Q_{\mathtt{f}},f\rangle^{\ell}L_{\ell}({\bf
r},{\mathbf{\Omega}},E)\,\varphi_{0}^{\dagger}({\bf r},{\mathbf{\Omega}},E)$

for large $n$. Here $L_{\ell}({\bf r},{\mathbf{\Omega}},E)$ is defined
recursively with $L_{1}=1$ and for $\ell\geq 2$,

(21) $L_{\ell}({\bf
r},{\mathbf{\Omega}},E)=\frac{{\ell!}}{\varphi_{0}^{\dagger}({\bf
r},{\bf\Omega},E)}\sum_{m=1}^{\infty}\frac{1}{k_{0}^{m\ell}}\Psi_{m}^{-}\Bigg{[}\mathcal{E}_{\cdot}\Bigg{[}\sum_{[\ell_{1},\dots,\ell_{M}]_{\ell}^{2+}}\prod_{\begin{subarray}{c}j=1\\\
\ell_{j}>0\end{subarray}}^{M}L_{\ell_{j}}(\cdot,{\mathbf{\Omega}}_{j},E_{j})\varphi_{0}^{\dagger}(\cdot,{\mathbf{\Omega}}_{j},E_{j})\Bigg{]}\Bigg{]}({\bf
r},{\mathbf{\Omega}},E),$

where $[\ell_{1},\dots,\ell_{M}]_{\ell}^{2+}$ denotes the set of non-negative
tuples $(\ell_{1},\dots,\ell_{M})$, such that $\sum_{j=1}^{M}\ell_{j}=\ell$
and at least two of the $\ell_{j}$ are strictly positive. We have used the
notation $\Psi_{n}^{-}$ for the semigroup associated to the collection of
particles stopped just before the $n$-th fission event in their genealogical
lines of descent.

Again, in the case where $\ell=1$, we recover Eq. (12). Now considering the
case where $\ell=2$ and $f=1$, we have

(22) $\mathbb{E}_{\delta_{({\bf r},{\bf\Omega},E)}}[N_{n}^{2}]\sim
k_{0}^{2n}L_{2}({\bf r},{\bf\Omega},E)\varphi_{0}^{\dagger}(\bf
r,\bf\Omega,E),$

where we have used the fact that $\langle Q_{\mathtt{f}},{\bf 1}\rangle=1$
from Eq. (6). Since $L_{1}=1$, it follows that $L_{2}$ is given by

(23) $L_{2}({\bf r},{\bf\Omega},E)=\frac{2}{\varphi_{0}^{\dagger}({\bf
r},{\bf\Omega},E)}\sum_{m=1}^{\infty}\frac{1}{k_{0}^{2m}}\Psi_{m}^{-}\Bigg{[}\mathcal{E}_{\cdot}\Bigg{[}\sum_{[\ell_{1},\dots,\ell_{M}]_{2}^{2+}}\prod_{\begin{subarray}{c}j=1\\\
\ell_{j}>0\end{subarray}}^{M}\varphi_{0}^{\dagger}(\cdot,{\mathbf{\Omega}}_{j},E_{j})\Bigg{]}\Bigg{]}({\bf
r},{\mathbf{\Omega}},E).$

###  3.3 The subcritical case

For the subcritical case, with $k_{0}<1$, the findings of Sec. A show that

(24) $\Psi_{n}^{(\ell)}[f]({\bf r},{\mathbf{\Omega}},E)\sim
k_{0}^{n}\,L_{\ell}\,\varphi_{0}^{\dagger}({\bf r},{\mathbf{\Omega}},E)$

for large $n$, where again $L_{\ell}$ is defined recursively, with
$L_{1}=\langle Q_{\mathtt{f}},f\rangle$ and

(25) $L_{\ell}=\langle
Q_{\mathtt{f}},f^{\ell}\rangle+\sum_{m=0}^{\infty}\frac{1}{k_{0}^{m}}\left\langle\varphi_{0},\Sigma_{\mathtt{f}}(\cdot)\mathcal{E}_{\cdot}\Bigg{[}\sum_{[\ell_{1},\dots,\ell_{M}]_{\ell}^{2+}}{\ell\choose\ell_{1}\dots\ell_{M}}\prod_{{\begin{subarray}{c}j=1\\\
\ell_{j}>0\end{subarray}}}^{M}\Psi_{m}^{(\ell_{j})}[f](\cdot,{\mathbf{\Omega}}_{j},E_{j})\Bigg{]}\right\rangle,$

for $\ell\geq 2$. The set $[\ell_{1},\dots,\ell_{M}]_{\ell}^{2+}$ is defined
as in the previous section.

As with the other two cases, setting $\ell=1$ we recover Eq. (11). Setting
$f=1$ and $\ell=2$ and using Eq. (6), we have

(26) $\mathbb{E}_{\delta_{({\bf r},{\bf\Omega},E)}}[N_{n}^{2}]\sim
k_{0}^{n}L_{2}\varphi_{0}^{\dagger}({\bf r},{\bf\Omega},E),$

where

(27)
$L_{2}=1+2\sum_{m=0}^{\infty}\frac{1}{k_{0}^{m}}\left\langle\varphi_{0},\Sigma_{\mathtt{f}}(\cdot)\mathcal{E}_{\cdot}\Bigg{[}\sum_{[\ell_{1},\dots,\ell_{M}]_{2}^{2+}}\prod_{{\begin{subarray}{c}j=1\\\
\ell_{j}>0\end{subarray}}}^{M}\mathbb{E}_{\delta_{(\cdot,{\mathbf{\Omega}}_{j},E_{j})}}[N_{m}]\Bigg{]}\right\rangle.$

######  Remark 1.

In the case where we consider delayed neutron production, as well as prompt
neutrons, the fission operator $F^{\dagger}$ would be replaced by

$\displaystyle F^{\dagger}[g]({\bf r},{\mathbf{\Omega}},E)$
$\displaystyle=\nu_{\mathtt{f,p}}({\bf r},E)\Sigma_{\mathtt{f}}({\bf
r},E)\int_{E_{\mathtt{min}}}^{E_{\mathtt{max}}}\int_{\mathbb{S}_{2}}g({\bf
r},{\mathbf{\Omega}}^{\prime},E^{\prime})\chi_{\mathtt{f,p}}({\bf
r},{\mathbf{\Omega}}\to{\mathbf{\Omega}}^{\prime},E\to E^{\prime}){\rm
d}{\mathbf{\Omega}}^{\prime}{\rm d}E^{\prime}$
$\displaystyle\quad+\sum_{i=1}^{D}\nu_{\mathtt{f,d,i}}({\bf
r},E)\Sigma_{\mathtt{f}}({\bf
r},E)\int_{E_{\mathtt{min}}}^{E_{\mathtt{max}}}\int_{\mathbb{S}_{2}}g({\bf
r},{\mathbf{\Omega}}^{\prime},E^{\prime})\chi_{\mathtt{f,d,i}}({\bf
r},{\mathbf{\Omega}}\to{\mathbf{\Omega}}^{\prime},E\to E^{\prime}){\rm
d}{\mathbf{\Omega}}^{\prime}{\rm d}E^{\prime},$

where $\nu_{\mathtt{f,p}}({\bf r},E)$ denotes the average number of prompt
neutrons produced at a fission event whose directions and energies are
randomly distributed according to a (normalised) probability density
$\chi_{\mathtt{f,p}}({\bf
r},{\mathbf{\Omega}}\to{\mathbf{\Omega}}^{\prime},E\to E^{\prime})$, and
similarly $\nu_{\mathtt{f,d,i}}({\bf r},E)$ denotes the average number of
delayed neutrons coming from the $i$-th precursor, whose directions and
energies are distributed according to $\chi_{\mathtt{f,d,i}}({\bf
r},{\mathbf{\Omega}}\to{\mathbf{\Omega}}^{\prime},E\to E^{\prime})$, for
$i=1,\dots,D$. Since, in generational time, both prompt and delayed neutrons
are produced on the same time scales, our results remain unchanged. More
precisely, we could let $M_{i}$ denote the random number of delayed neutrons
produced from the $i$-th precursor and $M_{0}$ denote the number of prompt
neutrons produced, so that our hold in the setting of delayed neutrons by
taking $M=M_{0}+\dots+M_{D}$.

##  4 Rod model

In this section we illustrate the main statements of Sec. 3 by considering the
one-dimensional rod model, which is a highly simplified neutron transport
configuration where analytical results can be easily established [36]. As
such, the rod model provides an ideal benchmark framework for the verification
of Monte Carlo simulations; contrary to the infinite-medium two-group setup,
the rod model in particular offers the possibility of effectively probing
spatial and angular effects.

The rod model, which was introduced by the pioneering work by G. M. Wing [31],
assumes that neutron displacements are restricted to a straight line, with the
only permissible directions being $\\{-1,+1\\}$. For the benchmark considered
here, we take the viable spatial domain to be the bounded segment $D=(-R,R)$
for some $R>0$, so that the spatial position of the particle can be described
using the scalar coordinate $x\in(-R,R)$. For the sake of simplicity, we will
further assume that all neutrons have unit energy. Thus, in this case, the
phase space is reduced to $\mathcal{S}=D\times\\{-1,+1\\}$. Moreover, we
impose leakage boundary conditions at $x=-R$ and $x=R$.

As for the nuclear data, we take spatially-homogeneous cross sections

(28)
$\displaystyle\qquad\Sigma_{\mathtt{s}}(x)\equiv\Sigma_{\mathtt{s}}>0,\quad\Sigma_{\mathtt{f}}(x)\equiv\Sigma_{\mathtt{f}}>0,\quad\Sigma_{\mathtt{c}}(x)\equiv\Sigma_{\mathtt{c}}>0.$

The scattering and fission distributions are assumed to be isotropic:

(29) $\displaystyle
f_{\mathtt{s}}(x,\Omega\to\Omega^{\prime})=\chi_{\mathtt{f}}(x,\Omega\to\Omega^{\prime})=\frac{1}{2}\left(\delta(\Omega-\Omega^{\prime})+\delta(\Omega+\Omega^{\prime})\right).$

Finally, we assume that exactly two particles are emitted at fission events,
which imposes the average fission multiplicity $\nu_{\mathtt{f}}=2$.

For this choice of the physical parameters, the adjoint operators
$S^{\dagger}$ and $F^{\dagger}$ are given by

(30) $\displaystyle S^{\dagger}[g](x,\Omega)$
$\displaystyle=\frac{\Sigma_{\mathtt{s}}}{2}\big{(}g(x,\Omega)+g(x,-\Omega)\big{)}$
(31) $\displaystyle F^{\dagger}[g](x,\Omega)$
$\displaystyle=\frac{\nu_{\mathtt{f}}\Sigma_{\mathtt{f}}}{2}\big{(}g(x,\Omega)+g(x,-\Omega)\big{)}=\Sigma_{\mathtt{f}}\big{(}g(x,\Omega)+g(x,-\Omega)\big{)}.$

Correspondingly, the eigenvalue problem in Eq. (4) yields the set of two
coupled ordinary differential equations (ODE)

(32) $\displaystyle-\Omega$
$\displaystyle\frac{d}{dx}\varphi^{\dagger}(x,\Omega)+\Sigma_{\mathtt{t}}\varphi^{\dagger}(x,\Omega)-\frac{\Sigma_{\mathtt{s}}}{2}\big{(}\varphi^{\dagger}(x,\Omega)+\varphi^{\dagger}(x,-\Omega)\big{)}=\frac{\Sigma_{\mathtt{f}}}{k}(\varphi^{\dagger}(x,\Omega)+\varphi^{\dagger}(x,-\Omega)\big{)},$

for $\Omega=\pm 1$, with the boundary conditions

(33)
$\displaystyle\varphi^{\dagger}(R,\Omega=+1)=\varphi^{\dagger}(-R,\Omega=-1)=0.$

Standard ODE techniques (see e.g. Refs. 4, 36) show that $k_{0}$ is the
largest $k$ that satisfies the dispersion relation

(34)
$2\cos(2R\alpha_{k})=\left(\frac{\alpha_{k}}{\Sigma_{\mathtt{t}}}-\frac{\Sigma_{\mathtt{t}}}{\alpha_{k}}\right)\sin(2R\alpha_{k}),$

where

(35)
$\alpha_{k}=\sqrt{\Sigma_{\mathtt{t}}(\nu_{\mathtt{f}}\Sigma_{\mathtt{f}}/k+\Sigma_{\mathtt{s}}-\Sigma_{\mathtt{t}})}.$

Based on the dispersion law for $k$, it is possible to choose a combination of
nuclear data and system size such that $k_{0}$ is equal to some target value:
for a given $k_{0}$, we define

(36)
${\bar{c}}_{0}=\frac{\frac{\nu_{\mathtt{f}}\Sigma_{\mathtt{f}}}{k_{0}}+\Sigma_{\mathtt{s}}}{\Sigma_{\mathtt{t}}},$

with
$\Sigma_{\mathtt{t}}=\Sigma_{\mathtt{c}}+\Sigma_{\mathtt{s}}+\Sigma_{\mathtt{f}}$,
and the corresponding value of the system size must satisfy

(37)
$R_{0}=\frac{\arctan\left(\frac{1}{\sqrt{{\bar{c}}_{0}-1}}\right)}{\Sigma_{\mathtt{t}}\sqrt{{\bar{c}}_{0}-1}}.$

In particular, when $k_{0}=1$, the choice of $R_{0}$ corresponds to the
critical half-size of the rod. Since
$k_{0}<k_{\infty}:=\nu_{\mathtt{f}}\Sigma_{\mathtt{f}}/(\Sigma_{\mathtt{c}}+\Sigma_{\mathtt{f}})$,
$k_{\infty}$ being the infinite multiplication factor, a necessary condition
for the existence of a $R_{0}$ ensuring $k_{0}\geq 1$ is that the nuclear data
are chosen so that $k_{\infty}>1$.

The adjoint and forward dominant eigenmodes are readily obtained, recalling
that we impose the normalisation $\langle Q_{\mathtt{f}},\mathbf{1}\rangle=1$
and $\langle Q_{\mathtt{f}},\varphi_{0}^{\dagger}\rangle=1$, which fixes the
multiplicative constants. For the adjoint dominant eigenfunction we have

(38) $\varphi_{0}^{\dagger}(x,\Omega=\pm
1)=\frac{8\sin(\alpha_{k_{0}}R_{0})}{4R_{0}\alpha_{k_{0}}+2\sin(2\alpha_{k_{0}}R_{0})}\left[\cos(\alpha_{k_{0}}x)\mp\frac{\sin(\alpha_{k_{0}}x)}{\tan(\alpha_{k_{0}}R_{0})}\right],$

whereas for the forward dominant eigenfunction we have

(39) $\varphi_{0}(x,\Omega=\pm
1)=\frac{\alpha_{k_{0}}}{4\nu_{\mathtt{f}}\Sigma_{\mathtt{f}}\sin(\alpha_{k_{0}}R_{0})}\left[\cos(\alpha_{k_{0}}x)\pm\frac{\sin(\alpha_{k_{0}}x)}{\tan(\alpha_{k_{0}}R_{0})}\right],$

with

(40)
$Q_{\mathtt{f}}(x,\Omega)=\frac{\alpha_{k_{0}}}{4\sin(\alpha_{k_{0}}R_{0})}\cos(\alpha_{k_{0}}x).$

###  4.1 The statistical behavior of the fission chains

Knowledge of the dominant eigenfunctions allows one to obtain explicit results
for the asymptotic moments of the neutron population at large $n$. For the
single-speed rod model, the semigroup of the branching process is given by

(41)
$\Psi_{n}[f](x,\Omega)=\mathbb{E}_{\delta_{(x,\Omega)}}\left[\sum_{i=1}^{N_{n}}f(x_{i}^{(n)},\Omega_{i}^{(n)})\right].$

In the following we will focus in particular on the case $f=1$, i.e., the
counting process for the fission neutrons being in the rod at a given
generation $n$.

####  4.1.1 The critical case

Let us consider the critical case, with $k_{0}=1$. Recalling that $N_{n}$
denotes the number of neutrons generated within the rod in the $n$-th fission
generation, Eq. (11) shows that

(42) $\Psi_{n}[{\bf
1}](x,\Omega)=\mathbb{E}_{\delta_{(x,\Omega)}}[N_{n}]\sim\langle
Q_{\mathtt{f}},\mathbf{1}\rangle\varphi_{0}^{\dagger}(x,\Omega),$

as $n\to\infty$. Since the normalization has been chosen so that $\langle
Q_{\mathtt{f}},\mathbf{1}\rangle=1$, we therefore have

(43)
$\mathbb{E}_{\delta_{(x,\Omega)}}[N_{n}]\sim\varphi_{0}^{\dagger}(x,\Omega).$

Figure 1: Rod model benchmark for critical conditions, with $k_{0}=1$:
analysis of the average number of fission neutrons. Comparison between the
exact expression for $\varphi_{0}^{\dagger}(x,\Omega)$ in Eq. (42) (solid
lines) and Monte Carlo simulations for
$\mathbb{E}_{\delta_{(x,\Omega)}}[N_{n}]$ (symbols). Error bars correspond to
one sigma uncertainty. Shades of blue: case $\Omega=+1$, with increasing
$n=10$, $n=50$, $n=100$, and $n=300$; shades of red: case $\Omega=-1$, with
increasing $n=10$, $n=50$, $n=100$, and $n=300$.

For illustration, in Fig. 1 we display the comparison between the exact result
in Eq. (42) and Monte Carlo simulations. We have set the cross sections as

(44)
$\Sigma_{\mathtt{s}}=0.3,\quad\Sigma_{\mathtt{c}}=0.2,\quad\Sigma_{\mathtt{f}}=0.7,$

with $2$ secondary neutrons per fission events, and we have imposed
$D=(-R_{0},R_{0})$ where $R_{0}$ was given in Eq. (37) in order to ensure
$k_{0}=1$. Monte Carlo simulations have been run using $10^{6}$ particles,
with $n$ varying from $n=10$ to $n=300$. The Monte Carlo tally corresponds to
the ensemble-averaged number of fission neutrons
$\mathbb{E}_{\delta_{(x,\Omega)}}[N_{n}]$ being in the rod at generation $n$,
for a single ancestor neutron starting with coordinates $(x,\Omega)$. The
statistical agreement between Monte Carlo results and Eq. (42) is very good.
The asymptotic behavior is attained for relatively small values of $n$: at
$n=10$, the average number of fission neutrons has already settled to the
spatial shape given by Eq. (42).

Now let us consider the second moment of the number of fission neutrons in the
rod. Recall from Eq. (17) that

(45) $\mathbb{E}_{\delta_{(x,\Omega)}}[N_{n}^{2}]\sim
n\Sigma_{\mathtt{f}}\langle\varphi_{0},\mathcal{V}[\varphi_{0}^{\dagger}]\rangle\varphi_{0}^{\dagger}(x,\Omega),$

as $n\to\infty$. Here we have used the fact that the system is spatially
homogeneous, which implies that
$\langle\Sigma_{\mathtt{f}}\varphi_{0},\mathcal{V}[\varphi_{0}^{\dagger}]\rangle=\Sigma_{\mathtt{f}}\langle\varphi_{0},\mathcal{V}[\varphi_{0}^{\dagger}]\rangle$.
Moreover, combining Eq. (16) and the fact that fission is isotropic, we have

$\displaystyle\mathcal{V}[g](x,\Omega)$
$\displaystyle=\mathcal{E}_{(x,\Omega)}[g^{2}(x,\Omega_{1})+2g(x,\Omega_{1})g(x,\Omega_{2})+g^{2}(x,\Omega_{2})]-\mathcal{E}_{(x,\Omega)}[g^{2}(x,\Omega_{1})+g^{2}(x,\Omega_{2})]$
$\displaystyle=2\mathcal{E}_{(x,\Omega)}[g(x,\Omega_{1})g(x,\Omega_{2})]$
$\displaystyle=2\mathcal{E}_{(x,\Omega)}[g(x,\Omega_{1})]^{2}$ (46)
$\displaystyle=\frac{1}{2}\big{[}g(x,\Omega)+g(x,-\Omega)\big{]}^{2}.$

Based on this expression, we can easily derive

(47)
$\langle\varphi_{0},\mathcal{V}[\varphi_{0}^{\dagger}]\rangle=\frac{16\left[5+\cos(2\alpha_{k_{0}}R_{0})\right]\sin^{2}(\alpha_{k_{0}}R_{0})}{6\Sigma_{\mathtt{f}}\left[2\alpha_{k_{0}}R_{0}+\sin(2\alpha_{k_{0}}R_{0})\right]^{2}}.$

For illustration, in Fig. 2 we display the comparison between the exact result
in Eq. (45), normalized to $n$, and Monte Carlo simulations, also normalized
to $n$. Cross sections and $R_{0}$ are the same as in the previous numerical
example. Monte Carlo simulations have been run using $10^{6}$ particles, with
$n$ varying from $n=10$ to $n=300$. The Monte Carlo tally corresponds to the
ensemble-averaged normalized second moment
$\mathbb{E}_{\delta_{(x,\Omega)}}[N^{2}_{n}]/n$ of the number of fission
neutrons being in the rod at generation $n$, for a single ancestor neutron
starting with coordinates $(x,\Omega)$. The statistical agreement between
Monte Carlo results and Eq. (45) is again very good. The asymptotic behavior
is attained later than the case of the average value: the second moment of
fission neutrons is settled at $n=50$.

Figure 2: Rod model benchmark for critical conditions, with $k_{0}=1$:
analysis of the second moment of the number of fission neutrons. Comparison
between the exact expression for
$\Sigma_{\mathtt{f}}\langle\varphi_{0},\mathcal{V}[\varphi_{0}^{\dagger}]\rangle\varphi_{0}^{\dagger}(x,\Omega)$
in Eq. (45) (solid lines) and Monte Carlo simulations for
$\mathbb{E}_{\delta_{(x,\Omega)}}[N_{n}^{2}]/n$ (symbols). Error bars
correspond to one sigma uncertainty. Shades of blue: case $\Omega=+1$, with
increasing $n=10$, $n=50$, $n=100$, and $n=300$; shades of red: case
$\Omega=-1$, with increasing $n=10$, $n=50$, $n=100$, and $n=300$.

Finally, the survival probability for the critical rod model can be explicitly
computed using Eq. (18). Indeed, due to spatial homogeneity of the cross-
sections, we have

(48)
$\mathbb{P}_{\delta_{(x,\Omega)}}(N_{n}>0)\sim\frac{2}{n}\frac{\varphi_{0}^{\dagger}(x,\Omega)}{\Sigma_{\mathtt{f}}\langle\varphi_{0},\mathcal{V}[\varphi_{0}^{\dagger}]\rangle},$

which can be computed explicitly using Eqs. (38) and (47). The comparison
between Monte Carlo simulations and the exact formula in Eq. (48) is
illustrated in Fig. 3 as a function of the number of generations $n$. The
convergence to the asymptotic shape requires a larger $n$ compared to the
cases of the first and second moment of the number of fission neutrons: Monte
Carlo simulations attain the shape of Eq. (48) at about $n=300$.

Figure 3: Rod model benchmark for critical conditions, with $k_{0}=1$:
analysis of the survival probability. Comparison between the exact expression
for
$2\varphi_{0}^{\dagger}(x,\Omega)/(\Sigma_{\mathtt{f}}\langle\varphi_{0},\mathcal{V}[\varphi_{0}^{\dagger}]\rangle)$
in Eq. (48) (solid lines) and Monte Carlo simulations for
$n\mathbb{P}_{\delta_{(x,\Omega)}}(N_{n}>0)$ (symbols). Error bars correspond
to one sigma uncertainty. Shades of blue: case $\Omega=+1$, with increasing
$n=10$, $n=50$, $n=100$, and $n=300$; shades of red: case $\Omega=-1$, with
increasing $n=10$, $n=50$, $n=100$, and $n=300$.

####  4.1.2 The supercritical case

Next let us consider the supercritical case, with $k_{0}>1$. From Eqs. (11)
and (6), we have

(49) $\mathbb{E}_{\delta_{(x,\Omega)}}[N_{n}]\sim
k_{0}^{n}\varphi_{0}^{\dagger}(x,\Omega),$

with $\varphi_{0}^{\dagger}(x,\Omega)$ given explicitly in Eq. (38). A
comparison with respect to Monte Carlo simulations is illustrated in Fig. 4,
with the same cross sections as in the critical case; the rod half-length
$R_{0}$ has been adjusted according to Eq. (37) in order to ensure
$k_{0}=1.01$. Monte Carlo simulations are again run using $10^{6}$ particles,
with $n$ varying from $n=10$ to $n=100$. The Monte Carlo tally corresponds to
the normalized ensemble-averaged number of fission neutrons
$\mathbb{E}_{\delta_{(x,\Omega)}}[N_{n}]/k_{0}^{n}$ being in the rod at
generation $n$, for a single ancestor neutron starting with coordinates
$(x,\Omega)$. The statistical agreement between Monte Carlo results (also
normalized by the factor $k_{0}^{n}$) and Eq. (49) is very good. The
asymptotic behavior is attained for relatively small values of $n$: at $n=10$,
the average number of fission neutrons has already settled to the spatial
shape given by Eq. (49).

Figure 4: Rod model benchmark for supercritical conditions, with $k_{0}=1.01$:
analysis of the average number of fission neutrons. Comparison between the
exact expression for $\varphi_{0}^{\dagger}(x,\Omega)$ in Eq. (49) (solid
lines) and Monte Carlo simulations for
$\mathbb{E}_{\delta_{(x,\Omega)}}[N_{n}]/k_{0}^{n}$ (symbols). Error bars
correspond to one sigma uncertainty. Shades of blue: case $\Omega=+1$, with
increasing $n=10$, $n=50$, and $n=100$; shades of red: case $\Omega=-1$, with
increasing $n=10$, $n=50$, and $n=100$.

For the second moment, recall from Eq. (23) that

(50) $\mathbb{E}_{\delta_{(x,\Omega)}}[N_{n}^{2}]\sim
k_{0}^{2n}L_{2}(x,\Omega)\varphi_{0}^{\dagger}(x,\Omega)$

for sufficiently large $n$. Now, note that the combinatorial sum over the set
$[\ell_{1},\dots,\ell_{M}]_{2}^{2+}$ in Eq. (23) contains only one element,
$(1,1)$. Due to this and using the fact that fission is isotropic, the
aforementioned combinatorial sum collapses down to
$\frac{1}{2}\mathcal{V}[\varphi_{0}^{\dagger}]$, as given in Eq. (46). It
therefore follows that

(51)
$L_{2}(x,\Omega)=\sum_{m=1}^{\infty}\frac{\Psi_{m}^{-}[\mathcal{V}[\varphi_{0}^{\dagger}]](x,\Omega)}{k_{0}^{2m}\varphi_{0}^{\dagger}(x,\Omega)}.$

From a practical point of view, $L_{2}$ contains terms
$\Psi_{m}^{-}[\mathcal{V}[\varphi_{0}^{\dagger}]](x,\Omega)$, whose analytical
expression is known only for large $n$. However, since
$\Psi_{n}^{-}[\mathcal{V}[\varphi_{0}^{\dagger}]](x,\Omega)$ is the expected
value of $\mathcal{V}[\varphi_{0}^{\dagger}]$ evaluated over the coordinates
$(x_{i},\Omega_{i})$ of neutrons about to undergo a fission event at
generation $n\geq 1$ (see Eq. (21)), its expression can be estimated by Monte
Carlo simulation. For this purpose, we set

(52)
$L_{2}^{n}(x,\Omega)=\sum_{m=1}^{n}\frac{\Psi_{m}^{-}[\mathcal{V}[\varphi_{0}^{\dagger}]](x,\Omega)}{k_{0}^{2m}\varphi_{0}^{\dagger}(x,\Omega)},$

so that Eq. (50) can written in a more practical form:

(53)
$\frac{\mathbb{E}_{\delta_{(x,\Omega)}}[N_{n}^{2}]}{L_{2}^{n}(x,\Omega)}\sim
k_{0}^{2n}\varphi_{0}^{\dagger}(x,\Omega)$

for sufficiently large $n$.

Figure 5: Rod model benchmark for supercritical conditions, with $k_{0}=1.01$:
analysis of the second moment of the number of fission neutrons. Comparison
between the exact expression for $\varphi_{0}^{\dagger}(x,\Omega)$ in Eq. (53)
(solid lines) and Monte Carlo simulations for
$\mathbb{E}_{\delta_{(x,\Omega)}}[N_{n}^{2}]/(k_{0}^{2n}L_{2}^{n}(x,\Omega))$
(symbols). Error bars correspond to one sigma uncertainty. Shades of blue:
case $\Omega=+1$, with increasing $n=10$, $n=50$, and $n=100$; shades of red:
case $\Omega=-1$, with increasing $n=10$, $n=50$, and $n=100$.

A comparison between the analytical result in Eq. (53) and Monte Carlo
simulations is presented in Fig. 5, as a function of $n$. The system
parameters are the same as for the case of the first moment. Convergence to
the asymptotic shape is attained at about $n=50$, which is again later than
for the first moment. A good statistical agreement is found between the Monte
Carlo findings and the exact formula.

####  4.1.3 The subcritical case

We finally consider the subcritical case, with $k_{0}<1$. From Eqs. (11) and
(6), we have

(54) $\mathbb{E}_{\delta_{(x,\Omega)}}[N_{n}]\sim
k_{0}^{n}\varphi_{0}^{\dagger}(x,\Omega),$

where, as usual, $\varphi_{0}^{\dagger}(x,\Omega)$ is given explicitly in Eq.
(38). A comparison with respect to Monte Carlo simulations is illustrated in
Fig. 6, with the same cross sections as in the critical case; the rod half-
length $R_{0}$ has been adjusted according to Eq. (37) in order to ensure
$k_{0}=0.99$. Monte Carlo simulations are again run using $10^{6}$ particles,
with $n$ varying from $n=10$ to $n=100$. The Monte Carlo tally corresponds to
the normalized ensemble-averaged number of fission neutrons
$\mathbb{E}_{\delta_{(x,\Omega)}}[N_{n}]/k_{0}^{n}$ being in the rod at
generation $n$, for a single ancestor neutron starting with coordinates
$(x,\Omega)$. The statistical agreement between Monte Carlo results (also
normalized by the factor $k_{0}^{n}$) and Eq. (54) is very good. Similarly to
the supercritical case, the asymptotic behavior is attained for relatively
small values of $n$: at $n=10$, the average number of fission neutrons has
already settled to the spatial shape given by Eq. (54).

Figure 6: Rod model benchmark for subcritical conditions, with $k_{0}=0.99$:
analysis of the average number of fission neutrons. Comparison between the
exact expression for $\varphi_{0}^{\dagger}(x,\Omega)$ in Eq. (54) (solid
lines) and Monte Carlo simulations for
$\mathbb{E}_{\delta_{(x,\Omega)}}[N_{n}]/k_{0}^{n}$ (symbols). Error bars
correspond to one sigma uncertainty. Shades of blue: case $\Omega=+1$, with
increasing $n=10$, $n=50$, and $n=100$; shades of red: case $\Omega=-1$, with
increasing $n=10$, $n=50$, and $n=100$.

The second moment stems from Eq. (26):

(55) $\mathbb{E}_{\delta_{(x,\Omega)}}[N_{n}^{2}]\sim
k_{0}^{n}L_{2}\varphi_{0}^{\dagger}(x,\Omega).$

Again, noting that the combinatorial sum in the definition of $L_{2}$ given in
(27) contains one element, $(1,1)$, and using the fact that fission is
isotropic, yields

(56)
$L_{2}=1+2\Sigma_{\mathtt{f}}\sum_{m=0}^{\infty}\frac{\langle\varphi_{0},\mathcal{E}_{\cdot}[\Psi_{m}[\mathbf{1}](x,\Omega_{1})]^{2}\rangle}{k_{0}^{m}},$

where the average $\mathcal{E}$ is taken over the random post-fission
directions $\Omega_{1}$. Similarly to the case of the supercritical
configuration, the term $L_{2}$ contains the expectations
$\Psi_{m}[\mathbf{1}](x,\Omega_{1})$, whose analytical expression is only
known for large $m$. Again, this term can be estimated by Monte Carlo. To this
end, we set

(57)
$L_{2}^{n}=1+2\Sigma_{\mathtt{f}}\sum_{m=0}^{n}\frac{\langle\varphi_{0},\mathcal{E}_{\cdot}[\Psi_{m}[\mathbf{1}](x,\Omega_{1})]^{2}\rangle}{k_{0}^{m}}$

for a given fission generation $n$. The inner product appearing in the
numerator can be written more explicitly:

(58)
$\displaystyle\langle\varphi_{0},\mathcal{E}_{\cdot}[\Psi_{m}[\mathbf{1}](x,\Omega_{1})]^{2}\rangle$
$\displaystyle=\frac{1}{4}\int_{-L_{0}}^{L_{0}}(\varphi_{0}(x,+1)+\varphi_{0}(x,-1))\left(\mathbb{E}_{\delta_{(x,+1)}}[N_{m}]+\mathbb{E}_{\delta_{(x,-1)}}[N_{m}]\right)^{2}dx,$

where now the expected values in the integrand can be easily estimated by
Monte Carlo and then weighted by the forward eigenfunction given in Eq. (39).
With this definition, we therefore have

(59) $\frac{\mathbb{E}_{\delta_{(x,\Omega)}}[N_{n}^{2}]}{L_{2}^{n}}\sim
k_{0}^{n}\varphi_{0}^{\dagger}(x,\Omega).$

Figure 7: Rod model benchmark for subcritical conditions, with $k_{0}=0.99$:
analysis of the second moment of the number of fission neutrons. Comparison
between the exact expression for $\varphi_{0}^{\dagger}(x,\Omega)$ in Eq. (59)
(solid lines) and Monte Carlo simulations for
$\mathbb{E}_{\delta_{(x,\Omega)}}[N_{n}^{2}]/(k_{0}^{n}L_{2}^{n})$ (symbols).
Error bars correspond to one sigma uncertainty. Shades of blue: case
$\Omega=+1$, with increasing $n=10$, $n=50$, and $n=100$; shades of red: case
$\Omega=-1$, with increasing $n=10$, $n=50$, and $n=100$.

A comparison between the analytical result in Eq. (59) and Monte Carlo
simulations is presented in Fig. 7, as a function of $n$. The system
parameters are the same as for the case of the first moment. Convergence to
the asymptotic shape is attained at about $n=50$, which is again later than
for the first moment. A good statistical agreement is found between the Monte
Carlo findings and the exact formula.

##  5 Conclusions

In this work we have proposed a general formalism to characterize the moments
of the neutron population as a function of fission generations, with a
particular emphasis on determining reference solutions for Monte Carlo
simulations involving adjoint-weighted parameters in $k$-eigenvalue
calculations. In particular, precise asymptotics have been established for the
average and variance of the neutron population for sub-, super- and critical
configurations.

In view of using these findings as a test-bed for the verification of Monte
Carlo developments based on the Iterated Fission Probability (IFP) method,
which is the most widely adopted strategy to estimate adjoint-weighted
quantities in production particle-transport codes, we have introduced a
benchmark using the one-dimensional and single-speed ‘rod-model’. Within this
framework, we have determined reference solutions for the asymptotic average
number of neutrons stemming from an ancestor particle, as well as for the
second moment of this counting process. While the proposed benchmark is
necessarily simple in order to allow for analytical solutions, the statements
concerning the statistical moments are fairly broad and can be thus extended
to more complex configurations.

Future work will concern the generalization of the results proposed in this
paper in several directions: exact results for a benchmark configuration using
continuous-energy nuclear data might be derived, at least introducing suitable
hypotheses on the spatial dependence of the system; furthermore, our strategy
might be successfully applied with a limited amount of modifications to the
investigation of the generalized IFP method associated to $\alpha$-eigenvalue
problems [29].

##  Appendix A General statements of the main results

In this section we present the precise statements for the moment asymptotics
given informally in Eqs. (15), (20) and (24), as well as the theorem for the
survival probability given in Eq. (18) and the Yaglom limit in Eq. (19). The
proofs can be found in Ref. 5.

Before stating the theorems, we first introduce some assumptions on the model,
using the notation defined in Sec. 2.

  * (H1)

The cross sections $\Sigma_{\mathtt{s}}$, $\Sigma_{\mathtt{f}}$,
$\Sigma_{\mathtt{a}}$ and $\nu_{\mathtt{f}}$ are uniformly bounded from above.

  * (H2)

$\inf\big{(}\Sigma_{\mathtt{f}}({\bf r},E)\chi_{\mathtt{f}}({\bf
r},{\mathbf{\Omega}}\to{\mathbf{\Omega}}^{\prime},E\to
E^{\prime})+\Sigma_{\mathtt{s}}({\bf r},E)f_{s}({\bf
r},{\mathbf{\Omega}}\to{\mathbf{\Omega}}^{\prime},E\to E^{\prime})\big{)}>0$,
where the infimum is taken over all ${\bf r}\in D$,
${\mathbf{\Omega}},{\mathbf{\Omega}}^{\prime}\in\mathbb{S}_{2}$ and
$E,E^{\prime}\in(E_{\mathtt{min}},E_{\mathtt{max}})$.

  * (H3)

There exists a constant $C>0$, such that for any $g\in
L_{+,1}^{\infty}(\mathcal{S})$, the space of non-negative, measurable
functions on $\mathcal{S}$ that are uniformly bounded by unity, we have

$\langle\varphi_{0},\Sigma_{\mathtt{f}}\mathcal{V}[g]\rangle\geq
C\langle\varphi_{0},{\color[rgb]{0,0,0}F[g]}\rangle^{2},$

where

(60) $\mathcal{V}[g]({\bf r},{\mathbf{\Omega}},E)=\mathcal{E}_{({\bf
r},{\mathbf{\Omega}},E)}\left[\left(\sum_{i=1}^{M}g({\bf
r},{\mathbf{\Omega}}_{i},E_{i})\right)^{2}\right]-\mathcal{E}_{({\bf
r},{\mathbf{\Omega}},E)}\left[\sum_{i=1}^{M}g^{2}({\bf
r},{\mathbf{\Omega}}_{i},E_{i})\right].$

Assumption (H1) ensures finite activity in finite time, and (H2) is an
irreducibility-type condition that ensures at least one of fission or scatter
can occur everywhere in the domain. It was shown in [3] that under these two
assumptions, (11) holds. Assumption (H3) can also be thought of as an
irreducibility or spread-out-ness condition at branching events. We note that
it is possible to weaken assumption (H3) at the expense of a more complicated
proof. However, this is outside the scope of this article and we refer the
reader to [11] for details on how this can be done in continuous time. Finally
we note that conditions (H1) and (H2) are clearly satisfied for the rod model
presented in Sec. 4, and we refer the reader to [10, §9] for verification of
(H3).

###  A.1 The critical case

Suppose $k_{0}=1$. We state three theorems, one each pertaining to the
asymptotic behaviour described in Eqs. (15), (18) and (19).

######  Theorem 1 (Critical case, $k_{0}=1$).

Suppose that (H1) and (H2) hold and $k_{0}=1$. Define

$\Delta_{n}^{(\ell)}=\sup_{({\bf r},{\mathbf{\Omega}},E)\in\mathcal{S},g\in
L_{+,1}^{\infty}(\mathcal{S})}\left|n^{-(\ell-1)}\varphi_{0}^{\dagger}({\bf
r},{\mathbf{\Omega}},E)^{-1}\Psi_{n}^{(\ell)}[g]({\bf
r},{\mathbf{\Omega}},E)-2^{-(\ell-1)}\ell!\,\langle
Q_{\mathtt{f}},g\rangle^{\ell}\langle\varphi_{0},\Sigma_{\mathtt{f}}\mathcal{V}[\varphi_{0}^{\dagger}]\rangle^{\ell-1}\right|,$

where $\mathcal{V}$ was defined in (60). Then, for all $\ell\geq 1$,

(61) $\sup_{n\geq 0}\Delta_{n}^{(\ell)}<\infty\text{ and
}\lim_{n\rightarrow\infty}\Delta_{n}^{(\ell)}=0.$

The next theorem shows that the survival probability behaves asymptotically
like $c/n$, for some constant $c>0$.

######  Theorem 2.

Suppose (H1), (H2) and (H3) hold, and that $k_{0}=1$. Then, for any $({\bf
r},{\mathbf{\Omega}},E)\in\mathcal{S}$, we have

(62) $\lim_{n\to\infty}n\mathbb{P}_{\delta_{({\bf
r},{\mathbf{\Omega}},E)}}(N_{n}>0)=\frac{2\varphi_{0}^{\dagger}({\bf
r},{\mathbf{\Omega}},E)}{\langle\varphi_{0},\Sigma_{\mathtt{f}}\mathcal{V}[\varphi_{0}^{\dagger}]\rangle},$

where $\mathcal{V}$ was defined in (60).

Finally, we show that the limiting distribution of the process normalised by
the current generation converges to an exponential random variable, when
conditioned on survival.

######  Theorem 3.

Suppose (H1), (H2) and (H3) hold, and that $k_{0}=1$. Then, for any $({\bf
r},{\mathbf{\Omega}},E)\in\mathcal{S}$ and any $f\in
L_{+,1}^{\infty}(\mathcal{S})$,

(63) $\lim_{n\to\infty}\left(\frac{1}{n}\sum_{i=1}^{N_{n}}f({\bf
r}_{i}^{(n)},{\mathbf{\Omega}}_{i}^{(n)},E_{i}^{(n)})\Bigg{|}N_{n}>0\right)=Y,$

where the above convergence is in distribution and $Y$ is an exponential
random variable with rate $2/\langle
Q_{\mathtt{f}},f\rangle\langle\varphi_{0},\Sigma_{\mathtt{f}}\mathcal{V}[\varphi_{0}^{\dagger}]\rangle$.

###  A.2 Non-critical cases

We will now consider the moment asymptotics when $k_{0}\neq 1$.

######  Theorem 4 (Supercritical case, $k_{0}>1$).

Suppose that (H1) and (H2) hold and $k_{0}>1$. Define

$\Delta_{n}^{(\ell)}=\sup_{({\bf r},{\mathbf{\Omega}},E)\in\mathcal{S},g\in
L_{+,1}^{\infty}(\mathcal{S})}\left|k_{0}^{-n\ell}\,\varphi_{0}^{\dagger}({\bf
r},{\mathbf{\Omega}},E)^{-1}\Psi_{n}^{(\ell)}[g]({\bf
r},{\mathbf{\Omega}},E)-\langle Q_{\mathtt{f}},g\rangle^{\ell}L_{\ell}({\bf
r},{\mathbf{\Omega}},E)\right|,$

where $L_{1}=1$ and, for $\ell\geq 2$, $L_{\ell}({\bf r},{\mathbf{\Omega}},E)$
is given by the recursion

(64) $L_{\ell}({\bf
r},{\mathbf{\Omega}},E)=\frac{\ell!}{\varphi_{0}^{\dagger}({\bf
r},{\mathbf{\Omega}},E)}\sum_{m=1}^{\infty}k_{0}^{-{\color[rgb]{0,0,0}\ell
m}}\Psi_{m}^{-}\Bigg{[}\mathcal{E}_{\cdot}\Bigg{[}\sum_{[\ell_{1},\dots,\ell_{M}]_{\ell}^{2+}}\prod_{\begin{subarray}{c}j=1\\\
\ell_{j}>0\end{subarray}}^{M}L_{\ell_{j}}(\cdot,{\mathbf{\Omega}}_{j},E_{j})\varphi_{0}^{\dagger}(\cdot,{\mathbf{\Omega}}_{j},E_{j})\Bigg{]}\Bigg{]}({\bf
r},{\mathbf{\Omega}},E),$

with $[\ell_{1},\dots,\ell_{M}]_{\ell}^{2+}$ defining the set of non-negative
tuples $(\ell_{1},\dots,\ell_{N})$, such that $\sum_{j=1}^{N}\ell_{j}=\ell$
and at least two of the $\ell_{j}$ are strictly positive.

Then, for all $\ell\geq 1$,

$\sup_{n\geq 0}\Delta_{n}^{(\ell)}<\infty\text{ and
}\lim_{n\rightarrow\infty}\Delta_{n}^{(\ell)}=0.$

Finally, we consider the subcritical case.

######  Theorem 5 (Subcritical case, $k_{0}<1$).

Suppose that (H1) and (H2) hold and $k_{0}<1$. Define

$\Delta_{n}^{(\ell)}=\sup_{({\bf r},{\mathbf{\Omega}},E)\in\mathcal{S},g\in
L_{+,1}^{\infty}(\mathcal{S})}\left|k_{0}^{-n}\,\varphi_{0}^{\dagger}({\bf
r},{\mathbf{\Omega}},E)^{-1}\Psi_{n}^{(\ell)}[g]({\bf
r},{\mathbf{\Omega}},E)-L_{\ell}\right|,$

where $L_{1}=\langle Q_{\mathtt{f}},g\rangle$ and for $\ell\geq 2$, $L_{\ell}$
is given by

$L_{\ell}=\langle
Q_{\mathtt{f}},g^{\ell}\rangle+\sum_{m=0}^{\infty}k_{0}^{-m}\left\langle\varphi_{0},\Sigma_{\mathtt{f}}(\cdot)\mathcal{E}_{\cdot}\Bigg{[}\sum_{[\ell_{1},\dots,\ell_{N}]_{\ell}^{2+}}{\ell\choose\ell_{1}\dots\ell_{N}}\prod_{{\begin{subarray}{c}j=1\\\
\ell_{j}>0\end{subarray}}}^{N}\Psi_{m}^{(\ell_{j})}[f](\cdot,{\mathbf{\Omega}}_{j},E_{j})\Bigg{]}\right\rangle,$

where $[\ell_{1},\dots,\ell_{N}]_{\ell}^{2+}$ is as in the previous theorem.

Then, for all $\ell\geq 1$,

$\sup_{n\geq 0}\Delta_{n}^{(\ell)}<\infty\text{ and
}\lim_{n\rightarrow\infty}\Delta_{n}^{(\ell)}=0.$

## Acknowledgments

This work was partially supported by the EPSRC grant EP/W026899/1. The authors
report there are no competing interests to declare.

## References

  * [1] E. Brun, F. Damian, C. M. Diop, E. Dumonteil, F.-X. Hugot, C. Jouanne, Y. K. Lee, F. Malvagi, A. Mazzolo, O. Petit, J.C. Trama, T. Visonneau, and A. Zoia.  TRIPOLI-4®, CEA, EDF and AREVA reference Monte Carlo code.  Ann. Nucl. Energy, 82:151-160, 2015.
  * [2] Timothy P. Burke and Brian C. Kiedrowski.  Monte Carlo perturbation theory estimates of sensitivities to system dimensions.  Nuclear Science and Engineering, 189(3):199–223, 2018.
  * [3] Alexander M. G. Cox, Emma Horton, Andreas E. Kyprianou, and Denis Villemonais.  Stochastic methods for neutron transport equation III: Generational many-to-one and k_eff.  SIAM Journal on Applied Mathematics, 81(3):982–1001, 2021.
  * [4] Alexander MG Cox, Simon C Harris, Andreas E Kyprianou, and Minmin Wang.  Monte carlo methods for the neutron transport equation.  SIAM/ASA Journal on Uncertainty Quantification, 10(2):775–825, 2022\.
  * [5] Eric Dumonteil, Emma Horton, Andreas E. Kyprianou, and Andrea Zoia.  Moment asymptotics for the neutron transport equation: proofs of the main results.  Preprint, 2022.
  * [6] Félix Foutel-Rodier and Emmanuel Schertzer.  Convergence of genealogies through spinal decomposition with an application to population genetics.  Probability Theory and Related Fields, 187(3):697–751, 2023.
  * [7] B.D. Ganapol.  Analytical benchmarks for nuclear engineering applications case studies in neutron transport theory.  Technical Report NEA/DB/DOC(2008)1, NEA No. 6292, OECD/NEA, 2008.
  * [8] Jochen Geiger.  Elementary new proofs of classical limit theorems for galton–watson processes.  Journal of Applied Probability, 36(2):301–309, 1999.
  * [9] Isaac Gonzalez, Emma Horton, and Andreas E. Kyprianou.  Asymptotic moments of spatial branching processes.  Probab. Theory Rel. Fields, 184(2):805–858, 2022.
  * [10] Simon C Harris, Emma Horton, Andreas E Kyprianou, and Minmin Wang.  Yaglom limit for critical nonlocal branching markov processes.  The Annals of Probability, 50:2373–2408, 2022.
  * [11] Emma Horton, Andreas E. Kyprianou, Martín-Chavez Pedro, Ellen Powell, and Victor Rivero.  Stability of (sub)critical non-local spatial branching processes with and without immigration.  Preprint, 2024.
  * [12] B. C. Kiedrowski and F. B. Brown.  Adjoint-Based k-Eigenvalue Sensitivity Coefficients to Nuclear Data Using Continuous-Energy Monte Carlo.  Nuclear Science and Engineering, 174(3):227–244, 2013.
  * [13] B. C. Kiedrowski, F. B. Brown, and P. P. H. Wilson.  Adjoint-Weighted Tallies for k-Eigenvalue Calculations with Continuous-Energy Monte Carlo.  Nuclear Science and Engineering, 168(3):226–241, 2011.
  * [14] Brian C. Kiedrowski.  Review of early 21st-century monte carlo perturbation and sensitivity techniques for k-eigenvalue radiation transport calculations.  Nuclear Science and Engineering, 185(3):426–444, 2017.
  * [15] Brian Christopher Kiedrowski.  Analytic, Infinite-medium Solutions for Point Reactor Kinetics Parameters and Reactivity Perturbations.  Technical Report LA-UR-10-01803, Los Alamos National Laboratory, Los Alamos, NM, USA, 2010.
  * [16] AN Kolmogorov.  Zur lösung einer biologischen aufgabe.  Comm. Math. Mech. Chebyshev Univ. Tomsk, 2(1):1–12, 1938.
  * [17] Gabriel Kooreman and David P. Griesheimer.  A fully analytic space-time depletion benchmark for fixed-source and eigenvalue calculations and its application to monte carlo simulation.  Annals of Nuclear Energy, 127:99–110, 2019.
  * [18] J.A. Kulesza et al.  MCNP code version 6.3.0 Theory & user manual.  Technical Report LA-UR-22-30006, Los Alamos National Laboratory, USA, 2022\.
  * [19] Iván Lux and László Koblinger.  Monte Carlo Particle Transport Methods: Neutron and Photon Calculations.  CRC Press, USA, 1 edition, 1991.
  * [20] Thomas W Mullikin.  Limiting distributions for critical multitype branching processes with discrete time.  Transactions of the American Mathematical Society, 106(3):469–494, 1963.
  * [21] Y. Nauchi and T. Kameyama.  Development of Calculation Technique for Iterated Fission Probability and Reactor Kinetic Parameters Using Continuous-Energy Monte Carlo Method.  Journal of Nuclear Science and Technology, 47:977–990, 2010.
  * [22] Ellen Powell.  An invariance principle for branching diffusions in bounded domains.  Probability Theory and Related Fields, 173(3-4):999–1062, 2019\.
  * [23] S.D. Richards et al.  MONK11A: status and plans for the MONK Monte Carlo code for criticality and reactor physics analyses.  In M&C 2019, Portland, USA, 2019.
  * [24] Simon D. Richards, Chris M.J. Baker, Adam J. Bird, Pat Cowan, Nigel Davies, Geoff P. Dobson, Tim C. Fry, Albrecht Kyrieleis, and Paul N. Smith.  Monk and mcbend: Current status and recent developments.  Annals of Nuclear Energy, 82:63–73, 2015.
  * [25] Avneet Sood, R.Arthur Forster, and D. Kent Parsons.  Analytical benchmark test set for criticality code verification.  Progress in Nuclear Energy, 42(1):55–106, 2003.
  * [26] Francesco Tantillo and Simon Richards.  Adjoint neutron flux estimator implementation and verification in the continuous energy monte carlo code monk.  Nuclear Engineering and Design, 356:110368, 2020.
  * [27] N. Terranova, D. Mancusi, and A. Zoia.  New perturbation and sensitivity capabilities in TRIPOLI-4.  Annals of Nuclear Energy, 121:335–349, 2018.
  * [28] Nicholas Terranova, Guillaume Truchet, Igor Zmijarevic, and Andrea Zoia.  Adjoint neutron flux calculations with tripoli-4: Verification and comparison to deterministic codes.  Annals of Nuclear Energy, 114:136–148, 2018.
  * [29] Nicholas Terranova and Andrea Zoia.  Generalized iterated fission probability for monte carlo eigenvalue calculations.  Annals of Nuclear Energy, 108:57–66, 2017.
  * [30] Guillaume Truchet, Pierre Leconte, Alain Santamarina, Emeric Brun, Frédéric Damian, and Andrea Zoia.  Computing adjoint-weighted kinetics parameters in Tripoli-4® by the Iterated Fission Probability method.  Annals of Nuclear Energy, 85:17–26, 2015.
  * [31] G.M. Wing.  An Introduction to Transport Theory.  Wiley, 1962.
  * [32] Akiva M Yaglom et al.  Certain limit theorems of the theory of branching random processes.  In Doklady Akad. Nauk SSSR (NS), volume 56, pages 795–798, 1947\.
  * [33] Andrea Zoia, Eric Dumonteil, and Alain Mazzolo.  Residence time and collision statistics for exponential flights: the rod problem revisited.  Physical Review E, 84(2):021139, 2011.
  * [34] Andrea Zoia, Eric Dumonteil, and Alain Mazzolo.  Counting statistics: a feynman-kac perspective.  Physical Review E, 85(1):011132, 2012.
  * [35] Andrea Zoia, Eric Dumonteil, and Alain Mazzolo.  Discrete feynman-kac formulas for branching random walks.  EPL (Europhysics Letters), 98(4):40012, 2012.
  * [36] Andrea Zoia, Eric Dumonteil, Alain Mazzolo, and Sameh Mohamed.  Branching exponential flights: travelled lengths and collision statistics.  Journal of Physics A: Mathematical and Theoretical, 45(42):425002, 2012.