# Percolation properties of the neutron population in nuclear reactors

Benjamin Dechenaux<EMAIL_ADDRESS>Thomas Delcambre  Institut de
Radioprotection et de Sûreté Nucléaire (IRSN)
PSN-RES/SNC/LN, F-92260, Fontenay-aux-Roses, France.  Eric Dumonteil  Institut
de Recherche sur les Lois Fondamentales de l’Univers
CEA, Université Paris-Saclay, 91191 Gif-sur-Yvette, France

###### Abstract

Abstract

Reactor physics aims at studying the neutron population in a reactor core
under the influence of feedback mechanisms, such as the Doppler temperature
effect. Numerical schemes to calculate macroscopic properties emerging from
such coupled stochastic systems however require to define intermediate
quantities (e.g. the temperature field), which are bridging the gap between
the stochastic neutron field and the deterministic feedback. By interpreting
the branching random walk of neutrons in fissile media under the influence of
a feedback mechanism as a directed percolation process and by leveraging on
the statistical field theory of birth death processes, we will build a
stochastic model of neutron transport theory and of reactor physics. The
critical exponents of this model, combined to the analysis of the resulting
field equation involving a fractional Laplacian will show that the critical
diffusion equation cannot adequately describe the spatial distribution of the
neutron population and shifts instead to a critical super-diffusion equation.
The analysis of this equation will reveal that non-negligible departure from
mean field behavior might develop in reactor cores, questioning the attainable
accuracy of the numerical schemes currently used by the nuclear industry.

††preprint: APS/123-QED

##  I  Introduction

At their inception in the 30’s, Monte Carlo algorithms intended to describe
the neutron transport in fissile media [1, 2] in support of the first
applications of nuclear power, might it be for the design of the first nuclear
weapons or in support of nuclear energy production. Since then, the nuclear
energy industry has used the so-called Monte Carlo criticality codes as
reference numerical schemes to solve the linear Boltzmann equation. Nuclear
safety demonstrations in particular rely on a Monte Carlo solving of this
equation with few, if any, hypotheses, while the use of deterministic codes is
preferred for an approximate but fast -and even online- solving of a
simplified version of this equation: the two (energy) group critical diffusion
equation. In both cases, the aim of such neutronics codes 111neutron transport
theory, also known as neutronics, aims at characterizing the properties of the
neutron gas in fissile media is to calculate the neutron spatial distribution
in the context of mean-field hypotheses. Indeed, even at the startup of a
nuclear reactor (i.e. when the neutron densities are small), such mean-field
equations are believed to accurately describe the average behavior of the
stochastic regime. Though, in the power regime of the reactor, stabilizing
effects come into action. For example, local power excursions and potential
sudden rise of local neutron densities are tempered by the Doppler broadening
of the neutron cross-sections (giving the probabilities that a neutron at a
given energy induce a particular reaction on a nuclei): the quantum resonances
of the neutron-nuclei system are broadened by the thermal agitation of the
nuclei from which results a modification of the neutron capture cross section.
The study of the properties of neutron transport in fissile media subject to
feedback mechanisms is called ’reactor physics’. Reactor physics numerical
solvers hence couple the mean-field equations of neutronics to thermal or
thermal-hydraulics solvers via the definition of intermediate quantities (such
as the temperature scalar field) to extract macroscopic measurable quantities
characterizing the neutron population (such as local neutron fluxes measured
by fission rate chambers for example).

The fluctuations of the neutron population were also largely investigated in
the past fifty years [4, 5, 6], while being assumed to be of marginal concern
for nuclear safety, since they vanish as the core power increases. However, a
few years ago, this paradigm was questioned [7, 8, 9] by the study of a
neutronics toy-model based on the celebrated branching Brownian motion [10],
that couples the diffusive random walk of neutrons [11], mimicking their
transport in the Brownian regime, to a Galton-Watson birth-death process,
reproducing the induced fission reaction with variable number of outgoing
neutrons. Indeed, the first moments of the master equation of this model
revealed that the spatial correlations within the neutron population might in
some cases diverge, invalidating the use of the mean-field equation. From a
physical point of view, this clustering phenomenon (already characterized at
the time in the theory of population ecology [12, 13, 14]) causes the
emergence of spatial patterns in the neutron gas. While always present and
increasing with time in 1 and 2 dimensional systems, the spatial correlations
however saturate in 3 dimensions. A dedicated experiment taking place at the
Rensselaer Poytechnic Institute was designed and took place in 2018 to
characterize such spatio-temporal correlations affecting the neutron
distributions. This experiment revealed the existence of such correlations and
triggered questions relative to their persistence in the power regime of
reactor physics where feedbacks occur [15].

In this paper, leveraging on the seminal approach of Janssen and De Dominicis
[16, 17], we will build a simplified model of reactor physics to investigate
the validity of mean field equations in the power regime of nuclear reactors.
In that aim, a field theoretic formulation of neutron physics will be built in
Section 2. While being 0 dimensional, this formulation will then be exploited
in Section 3 so as to take into account the principal features of neutronics
such as the calculation of fluctuations and correlations of the neutron field
or the effects of an external neutron source. It will also be extended in
order to add the effect of delayed neutrons in Section 4 and to take into
account spatial phenomena in Section 5. Along the path, we will highlight the
close relationship between the response functional formalism of Janssen and De
Dominics with the Doi and Peliti field theory [18, 19] which lend itself to an
exact microscopic interpretation, thus strengthening the use of both
approaches in neutronics. Finally, in Section 6, we will add a simplified
model of thermal feedback, showing that the neutron gas in a nuclear reactor
operated at criticality in the power regime can be described as a phase
transition belonging to the directed percolation universality class. The
calculation of the critical exponents at the percolation threshold using the
renormalization group will in return allow to specify the mean field equation,
formulated in terms of fractional Laplacian. We will in particular show that
in the directed percolation model, the critical diffusion equation shifts to a
critical super-diffusion equation, whose stochastic generators are Lévy
flights, ultimately dynamically modifying the medium properties in which
neutrons are propagating. The analysis of the fundamental mode of the
associated fractional Laplacian using an approximate formula will allow to
quantify the differences between this formal approach and classical coupling
schemes, thus setting an accuracy limit on all actual numerical solvers of
reactor physics.

##  II  A Field theoretic formulation of neutronics

Branching processes are among the simplest models capable of characterizing
the phenomenology of neutrons evolving in fissile materials [6, 5]. In their
simplest expression, these models approximate the path of neutrons in matter
as a stochastic process, in which they are subject to random extinction or
reproduction events, occurring at constant average rates. They are also known
to display a second order phase transition between an active and an inactive
state [10]. This transition is at the core root of nuclear reactors operations
and corresponds, in this context, to the appearance of self-sustained chain
reactions: the ’percolation threshold’ where neutron chains survive
indefinitely is precisely the point at which nuclear reactors do operate.

Field theoretic method have long been used in statistical and condensed matter
physics for their efficiency in the study of the intricate behavior of systems
undergoing phase transitions. They indeed prove particularly adapted to
unravel the universal properties and scaling behaviors in the vicinity of the
phase transition’s critical point. Along this line of thought, the present
paper ought to demonstrate that field theory also provides an adapted
framework for neutronics, allowing notably to unravel highly non-trivial
behaviors of the neutron population close to criticality.

The derivation of a field-theoretic version of neutronics followed in this
paper is based on the seminal work of Janssen [16] and De Dominicis [17].
Starting from the mean field equation of the problem, to which is supplemented
a random noise source term (as is commonly done to study the stability of
nuclear reactors [4]), the neutron population evolving in a nuclear reactor
can be decoupled into two components:

  * •

a smoothly varying neutron field $N(t)$. It corresponds to the instantaneous
number of neutrons that can be measured in a reactor (proportional to the
operating power for instance) and it is only sensitive to the global
reactivity $\rho$ of the reactor;

  * •

a random noise source $\eta(t)$, arising from various perturbations (such as
the stochasticity of the induced fission reactions that have a random number
of outgoing neutrons, or such as assembly vibrations) which are generally
driven by unknown or at least unspecified phenomena.

In the simplest model of the branching process associated to neutrons evolving
in an idealized, infinite and homogeneous reactor, such a decoupling
transforms the mean field equation into a Stochastic Differential Equation
(SDE), whose form is given by

$\cfrac{dN(t)}{dt}=\rho N(t)+\eta(t).$ (1)

Any macroscopic observable quantity $\mathcal{O}[N]$ is insensitive to the
rapidly varying and low lying noise term. As such, physical observables must
be averaged over all possible configurations of the noise. This is formally
achieved by performing a functional integration over the noise probability
functional $P[\eta(t)]$,

$\langle\mathcal{O}[N]\rangle\,\propto\,\int\mathcal{D}\eta\,\mathcal{O}[N]\,P[\eta].$
(2)

The development of the field-theoretic model starts by noting that the $N$
appearing under the integral sign of Eq. 2 is constrained to be a solution of
Eq. 1. The constraint (valid for all $t$),

$\mathcal{C}[N]\,=\,\left\\{\frac{dN(t)}{dt}\,-\,\rho\,N(t)\,-\,\eta(t)\right\\}=0,$
(3)

can be enforced with the help of a functional version of the resolution of the
identity [20],

$1=\int\mathcal{D}N\,\prod_{t}\delta\left(\mathcal{C}[N]\right).$ (4)

Performing a (functional) Fourier transform, one can write

$\displaystyle 1$ $\displaystyle=$
$\displaystyle\int\mathcal{D}[i\widetilde{N}]\int\mathcal{D}N\,e^{-\int
dt\widetilde{N}(t)\,\mathcal{C}[N]}$ (5) $\displaystyle=$
$\displaystyle\int\mathcal{D}[i\widetilde{N}]\int\mathcal{D}N\,e^{-\int
dt\widetilde{N}(t)\,\left(\frac{dN(t)}{dt}\,-\,\rho\,N(t)\,-\,\eta(t)\right)},$

where the purely imaginary auxiliary field $\widetilde{N}(t)$, known as the
Martin-Siggia-Rose response field [21], has been introduced. Inserting this
Fourier transformed identity back into Eq. 2, one obtains

$\displaystyle\langle\mathcal{O}[N]\rangle$ $\displaystyle\propto$
$\displaystyle\int\mathcal{D}[i\widetilde{N}]\int\mathcal{D}N\left\\{e^{-\int\widetilde{N}(t)\left(\frac{d}{dt}-\rho\right)N(t)\,dt}\right.$
(6) $\displaystyle\times$
$\displaystyle\mathcal{O}[N]\left.\times\int\mathcal{D}\eta\,e^{\int\widetilde{N}(t)\eta(t)dt}P[\eta]\right\\}.$

To go further, one needs to integrate over $\eta$ and thus specify the
probability functional to which it is associated. In the following, it will be
supposed that at any given moment in time, the noise follows a Gaussian
probability distribution

$P[\eta]\propto e^{-\frac{1}{4}\int\eta(t)\Gamma^{-1}\eta(t)\,dt}$ (7)

with zero mean $\langle\eta(t)\rangle=0$ and a covariance given by,

$\langle\eta(t)\eta(t^{\prime})\rangle=2\Gamma(N)\delta(t-t^{\prime}).$ (8)

Note that an explicit dependence of the variance $\Gamma(N)$ on the global
power level (i.e. $N(t)$) is to be expected on physical ground : the expected
noise in a reactor indeed depends on the operating power level. The integral
over the noise can now be straightforwardly evaluated as

$\int\mathcal{D}\eta\,e^{\int\widetilde{N}(t)\eta(t)dt}P[\eta]\propto
e^{\Gamma\widetilde{N}^{2}}.$ (9)

Inserting this result back into Eq. 6, one finally arrives at an expression
that solely depends on the variance of the noise probability distribution

$\langle\mathcal{O}[N]\rangle\,=\mathcal{N}\,\int\mathcal{D}[i\widetilde{N}]\int\mathcal{D}N\,\mathcal{O}[N]\,e^{-\mathcal{S}[N,\widetilde{N}]},$
(10)

where both a normalisation factor $\mathcal{N}$ and the so-called response
functional have been introduced, this last functional being defined through

$\mathcal{S}[N,\widetilde{N}]=\int\left\\{\widetilde{N}(t)\left(\frac{d}{dt}-\rho\right)N(t)-\widetilde{N}(t)\Gamma\widetilde{N}(t)\,\right\\}dt.$
(11)

It often reveals convenient to include in Eq. 3 an initial condition of the
type $N(t_{0})=N_{0}$. The initial condition is then included directly into
the response functional and it will be taken into account naturally in the
further calculation of observables. One can show that the inclusion of an
initial condition in the response functional gives

$\displaystyle\mathcal{S}[N,\widetilde{N}]$
$\displaystyle=\int\left\\{\widetilde{N}(t)\left(\frac{d}{dt}-\rho\right)N(t)\right.$
(12)
$\displaystyle\left.-\widetilde{N}(t)\,\delta(t-t_{0})\,N_{0}-\widetilde{N}(t)\Gamma\widetilde{N}(t)\,\right\\}dt.$

It can then be demonstrated that the net effect of this supplementary term
appearing in the response functional is to multiply each quantity evaluated at
the initial time by the factor $N_{0}$. In the remainder of the paper, only
the latter requirement will be retained and the terms corresponding to the
initial conditions will be systematically omitted in all of the response
functionals that will be considered.

Equations 10–12 form the field-theoretic version of the branching process. It
presents itself under the form of a (euclidean) path integral. The further
calculation of observable quantities can now be addressed thanks to the
powerful machinery developed in Quantum Field Theory (QFT) and statistical
mechanics.

##  III Calculation of observables

It can generally safely be assumed that in a nuclear reactor, the random power
fluctuations driven by random noise are much weaker than the actual operating
power level, i.e. the noise term is small. This suggests a strategy to extract
sensible results from the response functional formalism. One can expand the
function $\Gamma(N)$ into the Taylor series

$\Gamma(N)=\lambda_{0}+\lambda_{1}N+\lambda_{2}N^{2}\,+...,$ (13)

with small $\lambda_{i}$’s. Under this form, the noise term appearing in the
response functional can be treated perturbatively, just as one does with an
interaction potential in quantum mechanics. This perturbative treatment and
the calculation of observables of the type given in Eq. 10 is most
conveniently represented under the form of Feynman diagrams [22, 23].

The observables derived from the response functional will all depend on the
parameters $\lambda_{i}$, which play here the same role as the interaction
coupling constants in conventional QFT. The precise meaning of these
parameters and in particular their microscopic origin is, however, left
unspecified. In this respect, the response functional formalism can be viewed
as a purely phenomenological approach: the precise interpretation of the
coupling parameters is to be sought either a posteriori or from an external
knowledge of the problem.

One might thus rightly be worried that the calculation of observables relies
upon an infinite power series expansion, whose structure and convergence
depends on unknown parameters. In fact, the inclusion of powers higher than
one in Eq. 13 leads in some cases to diverging expressions, rendering the
perturbative expansion potentially inoperant (an example of which will be
encountered in the last section). This type of divergence is a common problem
in QFT and can generally be cured if one insists that all observables
quantities truly describe the macro-scale behaviour of the problem, i.e. that
all small scale fluctuations are truly being ’integrated over’. In a nutshell,
the observables must be derived at a fixed point of a spatio-temporal scaling
transformation of the system. By simple dimensional arguments [24], one finds
that high order terms in the Taylor expansion of the function $\Gamma$ dies
off when such a scaling operation is performed. In the Renormalization Group
(RG) jargon, they represent irrelevant parameters which formally vanish in the
vicinity of the RG fixed point. These higher order terms can thus be directly
discarded at this stage and $\Gamma$ may thus be simply expressed as

$\Gamma(N)=\lambda_{0}+\lambda_{1}N.$ (14)

###  III.1  Moments of the neutron population

Because the noise term has zero mean, it cannot impact the mean number of
neutrons $\langle N(t)\rangle=N(t)$. As long as mean values are concerned, one
can thus reasonably forget about the noise term in the response functional.
Using Eq. 10, the mean number of neutrons can be written as,

$\langle
N\rangle\,=\mathcal{N}\,\int\mathcal{D}[i\widetilde{N}]\int\mathcal{D}N\,N\,e^{-\int\left\\{\widetilde{N}(t)\left(\frac{d}{dt}-\rho\right)N(t)\right\\}dt}.$
(15)

This expression can be simply evaluated with the help of the functional
version of Wick’s theorem [23]. Taking into account for the initial condition
$N(0)=N_{0}$, one has,

$\langle N(t)\rangle=N_{0}\times\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman \vertex(a) {$t$};
\vertex[right=1.5cm of a](o1) {$t_{0}$}; \diagram* { (a) -- (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}=N_{0}\,\mathcal{G}(t-t_{0}).$
(16)

$\mathcal{G}(t-t_{0})=\langle N(t)\widetilde{N}(t_{0})\rangle$ is the retarded
Green function, solution of the equation

$\left(\frac{d}{dt}-\rho\right)\mathcal{G}(t-t^{\prime})=\delta(t-t^{\prime}),$
(17)

so that,

$\mathcal{G}(t-t^{\prime})=e^{\rho(t-t^{\prime})}\,\theta(t-t^{\prime}),$ (18)

with $\theta(t-t^{\prime})$ being the Heaviside function, equals to zero
whenever $t<t^{\prime}$ and one otherwise. The mean number of neutrons then
reads

$\langle N(t)\rangle=N_{0}\,e^{\rho(t-t_{0})}\,\theta(t-t_{0}).$ (19)

The second moment $\langle N^{2}(t)\rangle$ is less trivial, because it
depends on the noise variance $\Gamma(N)$. Forgetting temporarily the
$\lambda_{0}$ term in the series expansion of $\Gamma$, and resorting again to
Wick’s theorem, one can symbolically write

$\displaystyle\langle N^{2}(t)\rangle\,$ $\displaystyle=$ $\displaystyle
2\lambda_{1}\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[above left=of a](f1); \vertex[below left=of a](f2); \vertex[right=of
a](o1); \diagram* { (f1) -- (a) -- (f2), (a) -- (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+2\;\;\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[above= -0.3cm of a](c) ; \vertex[above= 0.4cm of a](b) ;
\vertex[right=of c](o1) ; \vertex[right=of b](o2) ; \diagram* { (c) -- (o1),
(b) -- (o2), }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}$ (20)
$\displaystyle=$ $\displaystyle\langle N^{2}(t)\rangle_{c}+\langle
N(t)\rangle^{2},$

where the notation $\langle\rangle_{c}$ is used to represent the first, so-
called connected, diagram. The variance can then be simply evaluated as
($\rho\neq 0$)

$\displaystyle\langle N^{2}\rangle-\langle N\rangle^{2}$ $\displaystyle=$
$\displaystyle\langle
N^{2}\rangle_{c}=2\lambda_{1}N_{0}\int_{t_{0}}^{t}d\tau\mathcal{G}^{2}(t-\tau)\mathcal{G}(\tau-
t_{0})$ (21) $\displaystyle=$
$\displaystyle\cfrac{2\lambda_{1}N_{0}}{\rho}\,e^{\rho(t-t_{0})}\left[e^{\rho(t-t_{0})}-1\right].$

This expression exactly matches the variance one would derive from a detailed
microscopic balance of the neutron population (in an infinite, homogeneous
medium) [6], provided that $\lambda_{1}$ be identified with the rescaled
variance $(\bar{\nu}^{(2)}-\rho)/2$ of the offspring distribution of the
underlying branching stochastic process. This result is a good illustration of
the general outreach of the response functional formalism. Relying only on the
SDE built upon the ’coarse-grained’ description of the problem, the formalism
is able to take advantage of the phenomenological and microscopic random
source term to extract residual information about the structure of the
fluctuations around the mean field behaviour of the system.

A strong limitation of this approach is the assumption that the noise term
follows a Gaussian profile. Indeed, the calculation of higher moments $\langle
N^{k}\rangle$ for the neutron population would formally require information
about the same moments of the offspring distribution. In the response
functional formalism, there is no possibility to include this type of
information in the calculations. This limitation does not however applies in
the present case, at least in the vicinity of criticality. It has indeed been
demonstrated in [25] that the microscopic branching process under study
displays a universal behaviour in the vicinity of the critical point $\rho=0$:
it then solely depends on the reactivity $\rho$ and the second moment of its
offspring distribution. Thus the Gaussian distribution imposed on the noise
source does not limit the extraction of higher moments of the neutron
distribution.

The approach proposed in [25] also relies on a field-theoretic formulation of
the stochastic branching process, pioneered in the work of Doi [18] and Peliti
[19]. Starting from a formal identification of the ladder operators
$a^{\dagger}$ and $a$ with the multiplication or derivation by $z$ of the
generating function $\psi(z,t)=\sum_{n}P(n,t)z^{n}$, one can give a structure
to the stochastic process that closely resemble the structure of Fock spaces
used in quantum mechanics [26]. A path integral formulation of the process is
then built upon the coherent states of the annihilation operator $a$ [27].

Both the Doi-Peliti and Janssen-De Dominicis approaches allow for the
derivation of a field theory, from two different starting points. In this
respect, they can be viewed as complementary:

  * •

the Doi-Peliti approach relies on the explicit modelling of the underlying
microscopic stochastic process. Every term in the theory is given a sound
microscopical interpretation from the start. The drawback is the difficulty to
include terms that cannot be described as stochastic processes (because their
microscopic origin is unclear for instance);

  * •

the Janssen-De Dominicis approach is much more phenomenological. Its starting
point is the mean field equation associated to the problem. This enables more
flexibility in the modelling of phenomena that could not be given a
microscopic origin (e.g. the modelling of counter-reaction mechanisms in
neutronics for instance). Relying on Gaussian noise distributions and
unspecified parameters, this approach is nevertheless more limited in scope
when a microscopic model is available.

Thanks to the universal behaviour of branching processes close to their
critical point, both approaches coincide in the present case: they both lead
to the same path integral formulation and the same diagrammatic expansion. As
such, the extraction of higher moments for the neutron population in the
response functional formalism exactly matches the derivation presented in
[25].

###  III.2 Addition of an external neutron source

At low power, external neutron sources, such as the spontaneous fission on
some heavy nuclei or $(\alpha,n)$ reactions might have a substantial impact on
the moments of the neutron population [5, 15]. The response functional
formalism lends itself easily to describe this phenomenon. Considering an
external neutron source, emitting an average of $S$ neutrons per second, the
SDE describing the problem now reads

$\cfrac{dN(t)}{dt}=\rho N(t)+S+\eta(t),$ (22)

as well as the associated response functional,

$\displaystyle\mathcal{S}[N,\widetilde{N}]$ $\displaystyle=$
$\displaystyle\int\left\\{\widetilde{N}(t)\left(\frac{d}{dt}-\rho\right)N(t)-\widetilde{N}S\right.$
(23)
$\displaystyle-\left.\widetilde{N}(t)\Gamma\widetilde{N}(t)\,\right\\}dt.$

Considering that the source term also comes with an intrinsic noise, we now
can set $\lambda_{0}\neq 0$ in the Taylor expansion of
$\Gamma(N)=\lambda_{0}+\lambda_{1}N$.

The moments of the neutron population must now take into account the external
source. In the perturbative calculation, this is done in the form of new
Feynman diagrams appearing in the calculations. For the mean number of
neutrons at time $t$, one has (setting $t_{0}=0$ in the following),

$\displaystyle\langle N(t)\rangle$ $\displaystyle=\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[right=of a](o1); \diagram* { (a) -- (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+S\,\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[right=of a,crossed dot] (o1) { }; \diagram* { (a) -- (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}$
$\displaystyle=N_{0}e^{\rho t}+S\int_{0}^{t}e^{\rho(t-t^{\prime})}dt^{\prime}$
(24) $\displaystyle=N_{0}e^{\rho t}+\frac{S}{\rho}\left(e^{\rho t}-1\right),$

in agreement with the literature [5]. The second propagator that now appears
in the diagram directly account for the creation of neutrons produced at a
constant rate $S$ by the source.

As for the variance of the neutron population, one has,

$\langle N^{2}(t)\rangle_{c}=2\lambda_{1}\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[above left=of a](f1); \vertex[below left=of a](f2); \vertex[right=of
a](o1); \diagram* { (f1) -- (a) -- (f2), (a) -- (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+2\lambda_{1}S\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[above left=of a](f1); \vertex[below left=of a](f2); \vertex[right=of
a, crossed dot](o1) { }; \diagram* { (f1) -- (a) -- (f2), (a) -- (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+2\lambda_{0}\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex[empty
dot] (a) { }; \vertex[above left=of a](f1); \vertex[below left=of a](f2);
\diagram* { (f1) -- (a) -- (f2), }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}$ (25)

The first diagram has already been calculated in the previous section. The
following two diagrams take into account respectively:

  * •

the variance of the neutrons’ branching process after their creation by the
source $S$;

  * •

the intrinsic variance $\lambda_{0}$ of the neutron source.

The calculation of the two unknown diagrams leads to the final expression
($\rho\neq 0$):

$\displaystyle\langle N^{2}(t)\rangle_{c}$ $\displaystyle=$
$\displaystyle\frac{2\lambda_{1}N_{0}}{\rho}e^{\rho t}\left(e^{\rho
t}-1\right)+\frac{2\lambda_{1}}{\rho}\left(e^{\rho t}-1\right)^{2}$ (26)
$\displaystyle+$ $\displaystyle\frac{\lambda_{0}}{\rho}\left(e^{2\rho
t}-1\right).$

The calculation of moments of higher orders now poses a real problem. Although
the universal behavior of the neutrons’ branching process could, here also, be
invoked, one is now limited by the noisy component of the external source
term. Only if the source possesses a Gaussian distribution can the calculation
be taken further and the moments $\langle N^{k}\rangle_{c}$ with $k>2$ be
exactly calculated.

##  IV Two species model

Another key aspect in the modelling of the temporal behavior of the neutron
population in a nuclear reactor is the role played by the precursors of
delayed neutrons $M(t)$. In rare cases indeed, an induced fission does not
emit all the outgoing neutrons instantaneously as the fission fragments can
also emit neutrons on longer time scale controlled by their electroweak decay,
ranging from a second to a few minutes, orders of magnitudes longer than
prompt neutrons and occurring in the present model at a rate $\lambda_{D}$.
Thus, while emitting a small fraction $\beta$of the overall number of
neutrons, precursors control the kinetics of nuclear reactor which could not
be operated on human time scales otherwise.

Including the coupling between the neutron and precursor populations, the
point-like model transforms to (with $r=\rho-\beta)$

$\begin{cases}\cfrac{dN(t)}{dt}&=rN(t)+\lambda_{D}M(t)+\eta_{1}(t),\\\ \\\
\cfrac{dM(t)}{dt}&=\beta N(t)-\lambda_{D}M(t)+\eta_{2}(t),\end{cases}$ (27)

where external sources have been neglected, for simplicity. As before, the two
random source terms are supposed to follow Gaussian probability distributions
with zero mean and a covariance given by the relation ($a,b=1,2$),

$\langle\eta_{a}(t)\eta_{b}(t^{\prime})\rangle=2\,\Gamma_{ab}(N,M)\,\delta(t-t^{\prime}),$
(28)

The derivation of the response functional associated to this coupled system of
equations still follows the same recipe as the one described in section II.
Four independent fields $\widetilde{N},\widetilde{M},M$ and $N$ are now
introduced into the theory and the response functional can be written
$\mathcal{S}=\mathcal{S}_{\text{neut.}}+\mathcal{S}_{\text{prec.}}+\mathcal{S}_{\Gamma}$,
with,

$\mathcal{S}_{\text{neut.}}=\int\left\\{\widetilde{N}(t)\left(\frac{d}{dt}-r\right)N(t)-\lambda_{D}\widetilde{N}(t)M(t)\right\\}dt,$
(29)

$\mathcal{S}_{\text{prec.}}=\int\left\\{\,\widetilde{M}(t)\left(\frac{d}{dt}+\lambda_{D}\right)M(t)-\beta\widetilde{M}(t)N(t)\right\\}dt,$
(30)

and, in its most general form,

$\displaystyle\mathcal{S}_{\Gamma}=\int\Gamma_{11}\widetilde{N}^{2}+\Gamma_{22}\widetilde{M}^{2}+2\Gamma_{12}\widetilde{N}\widetilde{M}$
(31) $\displaystyle=$
$\displaystyle\int\left\\{\lambda_{1}\widetilde{N}^{2}N+\lambda_{2}\widetilde{N}^{2}M+\lambda_{3}\widetilde{M}^{2}N\right.$
$\displaystyle+\left.\lambda_{4}\widetilde{M}^{2}M+2\lambda_{5}\widetilde{N}\widetilde{M}N+2\lambda_{6}\widetilde{N}\widetilde{M}M\right\\}dt.$

Two Green functions are to be associated to the free parts of the neutron and
precursor functionals,

$\mathcal{G}_{N}(t-t^{\prime})=\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman \vertex(a) {$t$};
\vertex[right=1.5cm of a](o1) {$t^{\prime}$}; \diagram* { (a) -- (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}=e^{r(t-t^{\prime})}\theta(t-t^{\prime})$
(32)

and

$\mathcal{G}_{M}(t-t^{\prime})=\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman \vertex(a) {$t$};
\vertex[right=1.5cm of a](o1) {$t^{\prime}$}; \diagram* { (a) -- [scalar] (o1)
}; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}=e^{-\lambda_{D}(t-t^{\prime})}\theta(t-t^{\prime}).$
(33)

Because of the couplings between the neutron and precursor populations
appearing in Eq. 29 and Eq. 30, observables must however be ultimately
constructed from a generalized set of propagators, which solve the full
coupled set of equations Eq. 27. These generalized propagators may be
symbolically written,

$\displaystyle\Sigma_{NN}(t,t^{\prime})=\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[small, inline=(a)]
\vertex(a); \vertex[right=0.5 cm of a,blob](b) {}; \vertex[right=0.7 cm of
b](o1); \diagram* { (a) -- (b) -- (o1) }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}=\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[right=of a](o1); \diagram* { (a) -- (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+\beta\,\lambda_{D}\,\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[small, inline=(a)]
\vertex(a); \vertex[right=0.5 cm of a,dot](b) {}; \vertex[right=0.7 cm of
b,dot](c) {}; \vertex[right=0.5 cm of c](o1); \diagram* { (a) -- (b) ,
(b)--[scalar] (c), (c) -- (o1) }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+\left(\beta\,\lambda_{D}\right)^{2}\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[small, inline=(a)]
\vertex(a); \vertex[right=0.5 cm of a,dot](b) {}; \vertex[right=0.7 cm of
b,dot](c) {}; \vertex[right=0.7 cm of c,dot](d) {}; \vertex[right=0.7 cm of
d,dot](e) {}; \vertex[right=0.5 cm of e](o1); \diagram* { (a) -- (b) ,
(b)--[scalar] (c), (c) -- (d) , (d)--[scalar] (e), (e) -- (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+...\,,\,$ (34)
$\displaystyle\Sigma_{MM}(t,t^{\prime})=\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[small, inline=(a)]
\vertex(a); \vertex[right=0.5 cm of a,blob](b) {}; \vertex[right=0.7 cm of
b](o1); \diagram* { (a) -- [scalar] (b) --[scalar] (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}=\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[right=of a](o1); \diagram* { (a) --[scalar] (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+\beta\,\lambda_{D}\,\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[small, inline=(a)]
\vertex(a); \vertex[right=0.5 cm of a,dot](b) {}; \vertex[right=0.7 cm of
b,dot](c) {}; \vertex[right=0.5 cm of c](o1); \diagram* { (a) --[scalar] (b) ,
(b)-- (c), (c) --[scalar] (o1) }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+\left(\beta\,\lambda_{D}\right)^{2}\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[small, inline=(a)]
\vertex(a); \vertex[right=0.5 cm of a,dot](b) {}; \vertex[right=0.7 cm of
b,dot](c) {}; \vertex[right=0.7 cm of c,dot](d) {}; \vertex[right=0.7 cm of
d,dot](e) {}; \vertex[right=0.5 cm of e](o1); \diagram* { (a) --[scalar] (b) ,
(b)-- (c), (c) --[scalar] (d) , (d)-- (e), (e) --[scalar] (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+...\,,\,$ (35)
$\displaystyle\Sigma_{NM}(t,t^{\prime})=\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[small, inline=(a)]
\vertex(a); \vertex[right=0.5 cm of a,blob](b) {}; \vertex[right=0.7 cm of
b](o1); \diagram* { (a) -- [scalar] (b) --(o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}=\beta\left\\{\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[right=of a,dot](b) {}; \vertex[right=of b](o1); \diagram* { (a)
--[scalar] (b), (b) -- (o1) }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+\beta\,\lambda_{D}\,\,\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[right=of a,dot](b) {}; \vertex[right=of b,dot](c) {}; \vertex[right=of
c,dot](d) {}; \vertex[right=of d](o1); \diagram* { (a) --[scalar] (b) , (b)--
(c), (c)-- [scalar] (d), (d) -- (o1) }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+...\right\\},$ (36)
$\displaystyle\Sigma_{MN}(t,t^{\prime})=\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[small, inline=(a)]
\vertex(a); \vertex[right=0.5 cm of a,blob](b) {}; \vertex[right=0.8 cm of
b](o1); \diagram* { (a) -- (b) -- [scalar] (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}=\lambda_{D}\left\\{\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[right=of a,dot](b) {}; \vertex[right=of b](o1); \diagram* { (a) --
(b), (b) --[scalar] (o1) }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+\beta\,\lambda_{D}\,\,\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[right=of a,dot](b) {}; \vertex[right=of b,dot](c) {}; \vertex[right=of
c,dot](d) {}; \vertex[right=of d](o1); \diagram* { (a) -- (b) , (b)--
[scalar](c), (c)-- [] (d), (d) --[scalar] (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys<EMAIL_ADDRESS>(37)

The calculation of each one of the $\Sigma_{IJ}$ involves an infinite series
of nested integrals. Fortunately, these infinite series can be summed over and
it is possible to obtain analytical results in all four cases (see annex A). A
much more convenient but indirect way to extract an analytical expression out
of these propagators is to realize that they all directly correspond to the
mean value of either the neutron or precursor populations for different
initial conditions. Indeed, setting $(N(t_{0}),M(t_{0}))=(N_{0},0)$, one can
show that,

$\langle N(t)\rangle_{(N_{0},0)}=N_{0}\times\Sigma_{NN}(t,t_{0}),$ (38)
$\langle M(t)\rangle_{(N_{0},0)}=N_{0}\times\Sigma_{NM}(t,t_{0}),$ (39)

whereas, for $(N(t_{0}),M(t_{0}))=(0,M_{0})$,

$\langle N(t)\rangle_{(0,M_{0})}=M_{0}\times\Sigma_{MN}(t,t_{0}),$ (40)
$\langle M(t)\rangle_{(0,M_{0})}=M_{0}\times\Sigma_{MM}(t,t_{0}),$ (41)

so that each of these propagators is the solution of the system of equations
Eq. 27 for a specific set of initial conditions (and with the stochastic terms
set to zero). This system of equations can be straightforwardly integrated and
solutions extracted. The explicit form of each generalized propagator is
provided in annex A.

Equipped with an analytical expression for the propagators of the theory, the
variance of the process may then be evaluated. Imposing the initial condition
$(N(0),M(0))=(N_{0},0)$ , the variance of the neutron and precursor
populations may be diagrammatically written,

$\displaystyle\langle N^{2}(t)\rangle_{c}$ $\displaystyle=$ $\displaystyle
2N_{0}\lambda_{1}\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[right= 0.5 cm of a,blob](b) {}; \vertex[above left=0.5 cm of
a,blob](f1) {}; \vertex[below left=0.5cm of a,blob](f2) {}; \vertex[above
left=0.7 cm of f1](ff1); \vertex[below left=0.7 cm of f2](ff2);
\vertex[right=0.7 cm of b](o1); \diagram* { (ff1) -- (f1) -- (a) -- (f2) --
(ff2), (a) -- (b) -- (o1) }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+2N_{0}\lambda_{2}\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[right= 0.5 cm of a,blob](b) {}; \vertex[above left=0.5 cm of
a,blob](f1) {}; \vertex[below left=0.5cm of a,blob](f2) {}; \vertex[above
left=0.7 cm of f1](ff1); \vertex[below left=0.7 cm of f2](ff2);
\vertex[right=0.7 cm of b](o1); \diagram* { (ff1) -- (f1) -- (a) -- (f2) --
(ff2), (a) -- [scalar] (b) -- (o1) }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+2N_{0}\lambda_{3}\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[right= 0.5 cm of a,blob](b) {}; \vertex[above left=0.5 cm of
a,blob](f1) {}; \vertex[below left=0.5cm of a,blob](f2) {}; \vertex[above
left=0.7 cm of f1](ff1); \vertex[below left=0.7 cm of f2](ff2);
\vertex[right=0.7 cm of b](o1); \diagram* { (ff1) -- (f1) --[scalar] (a)
--[scalar] (f2) -- (ff2), (a) -- (b) -- (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}$ (42)
$\displaystyle+2N_{0}\lambda_{4}\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[right= 0.5 cm of a,blob](b) {}; \vertex[above left=0.5 cm of
a,blob](f1) {}; \vertex[below left=0.5cm of a,blob](f2) {}; \vertex[above
left=0.7 cm of f1](ff1); \vertex[below left=0.7 cm of f2](ff2);
\vertex[right=0.7 cm of b](o1); \diagram* { (ff1) -- (f1) --[scalar] (a)
--[scalar] (f2) -- (ff2), (a) --[scalar] (b) -- (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+2N_{0}\lambda_{5}\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[right= 0.5 cm of a,blob](b) {}; \vertex[above left=0.5 cm of
a,blob](f1) {}; \vertex[below left=0.5cm of a,blob](f2) {}; \vertex[above
left=0.7 cm of f1](ff1); \vertex[below left=0.7 cm of f2](ff2);
\vertex[right=0.7 cm of b](o1); \diagram* { (ff1) -- (f1) --[scalar] (a) --
(f2) -- (ff2), (a) -- (b) -- (o1) }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+2N_{0}\lambda_{6}\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[right= 0.5 cm of a,blob](b) {}; \vertex[above left=0.5 cm of
a,blob](f1) {}; \vertex[below left=0.5cm of a,blob](f2) {}; \vertex[above
left=0.7 cm of f1](ff1); \vertex[below left=0.7 cm of f2](ff2);
\vertex[right=0.7 cm of b](o1); \diagram* { (ff1) -- (f1) --[scalar] (a) --
(f2) -- (ff2), (a) -- [scalar] (b) -- (o1) }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}},$
$\displaystyle\langle M^{2}(t)\rangle_{c}$ $\displaystyle=$ $\displaystyle
2N_{0}\lambda_{1}\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[right= 0.5 cm of a,blob](b) {}; \vertex[above left=0.5 cm of
a,blob](f1) {}; \vertex[below left=0.5cm of a,blob](f2) {}; \vertex[above
left=0.7 cm of f1](ff1); \vertex[below left=0.7 cm of f2](ff2);
\vertex[right=0.7 cm of b](o1); \diagram* { (ff1) --[scalar] (f1) -- (a) --
(f2) --[scalar] (ff2), (a) -- (b) -- (o1) }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+2N_{0}\lambda_{2}\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[right= 0.5 cm of a,blob](b) {}; \vertex[above left=0.5 cm of
a,blob](f1) {}; \vertex[below left=0.5cm of a,blob](f2) {}; \vertex[above
left=0.7 cm of f1](ff1); \vertex[below left=0.7 cm of f2](ff2);
\vertex[right=0.7 cm of b](o1); \diagram* { (ff1) --[scalar] (f1) -- (a) --
(f2) --[scalar] (ff2), (a) -- [scalar] (b) -- (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+2N_{0}\lambda_{3}\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[right= 0.5 cm of a,blob](b) {}; \vertex[above left=0.5 cm of
a,blob](f1) {}; \vertex[below left=0.5cm of a,blob](f2) {}; \vertex[above
left=0.7 cm of f1](ff1); \vertex[below left=0.7 cm of f2](ff2);
\vertex[right=0.7 cm of b](o1); \diagram* { (ff1) --[scalar] (f1) --[scalar]
(a) --[scalar] (f2) --[scalar] (ff2), (a) -- (b) -- (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}$ (43)
$\displaystyle+2N_{0}\lambda_{4}\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[right= 0.5 cm of a,blob](b) {}; \vertex[above left=0.5 cm of
a,blob](f1) {}; \vertex[below left=0.5cm of a,blob](f2) {}; \vertex[above
left=0.7 cm of f1](ff1); \vertex[below left=0.7 cm of f2](ff2);
\vertex[right=0.7 cm of b](o1); \diagram* { (ff1) --[scalar] (f1) --[scalar]
(a) --[scalar] (f2) --[scalar] (ff2), (a) --[scalar] (b) -- (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+2N_{0}\lambda_{5}\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[right= 0.5 cm of a,blob](b) {}; \vertex[above left=0.5 cm of
a,blob](f1) {}; \vertex[below left=0.5cm of a,blob](f2) {}; \vertex[above
left=0.7 cm of f1](ff1); \vertex[below left=0.7 cm of f2](ff2);
\vertex[right=0.7 cm of b](o1); \diagram* { (ff1) --[scalar] (f1) --[scalar]
(a) -- (f2) --[scalar](ff2), (a) -- (b) -- (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+2N_{0}\lambda_{6}\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[right= 0.5 cm of a,blob](b) {}; \vertex[above left=0.5 cm of
a,blob](f1) {}; \vertex[below left=0.5cm of a,blob](f2) {}; \vertex[above
left=0.7 cm of f1](ff1); \vertex[below left=0.7 cm of f2](ff2);
\vertex[right=0.7 cm of b](o1); \diagram* { (ff1) --[scalar] (f1) --[scalar]
(a) -- (f2) --[scalar] (ff2), (a) -- [scalar] (b) -- (o1) };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys<EMAIL_ADDRESS>
Just as it was the case in the last section, both expressions depend on the
set of parameters $\lambda_{i}$, which are not specified by the theory. With a
total of 6 unknown parameters for only two observables, it may seem at first
sight difficult to give a precise meaning to each term. However, the structure
of the diagrammatic expansion and the results of the one species model derived
earlier allows to do so, at least partially. In light of the results obtained
in the previous section, it might indeed be inferred that the diagrams
associated to $\lambda_{1}$ appearing in both expressions represent the
variance associated to the neutron’s offspring distribution. There is
furthermore no reason to expect that the value for $\lambda_{1}$ obtained in
the case of the one species model should be any different in the present
situation, except for the replacement of the reactivity $\rho$ with the prompt
reactivity $r=\rho-\beta$.

On the same model, the diagrams associated with $\lambda_{3}$ are most
certainly to be identified with the variance of the precursor emitted by
fission.

The diagrams associated to $\lambda_{2}$ and $\lambda_{4}$ share the same
output structure as for $\lambda_{1}$ and $\lambda_{3}$ but possess a
precursor at the origin of the branching term. This suggest that these
diagrams encode for a source of variance intrinsic to the precursor
population. The only process responsible for such kind of fluctuation in the
precursor population is their decay process which, being a Poisson process,
has a variance $\lambda_{D}$.

At this point, one can thus postulate,

$\begin{cases}\lambda_{1}=(\bar{\nu}_{n}^{(2)}-\rho+\beta)/2,\\\
\lambda_{3}=(\bar{\nu}_{m}^{(2)}+\beta)/2,\\\
\lambda_{2}=\lambda_{4}=\lambda_{D}/2,\end{cases}$ (44)

so that only two unknown parameters, $\lambda_{5}$ and $\lambda_{6}$, remain
to be determined for the two variances. Comparing equations Eq. 42-43 with a
microscopic balance of the neutron and precursor populations [6] for different
values of $\rho$, $\beta$ and $\lambda_{D}$, one gets (the offspring
distributions of neutrons and precursors were supposed independent),

$\begin{cases}\lambda_{5}=\overline{\nu_{m}(\nu_{n}-1)}=\bar{\nu}_{n}^{(2)}\frac{\beta}{1-\beta},\\\
\lambda_{6}=-\lambda_{D}.\end{cases}$ (45)

A comparison of the diagrammatic expansion with the results of the literature
is presented in Figure 1. The temporal evolution of the variance of both the
neutron and precursor populations are displayed, in all three regimes: sub-
critical ($\rho<0)$, critical $(\rho=0)$ and super-critical $(\rho>0)$. An
excellent agreement is found in all three cases.

Figure 1: Comparison of the variance of the neutron and precursor populations
extracted from the solution of the stochastic process (plain lines) and from
the calculation of Feynman diagrams (dotted lines), for $\rho=$ -0.2 (red), 0
(green), 0.2 (blue). The calculations are made assuming the following set of
parameters: $\beta=0.8$, $\lambda_{D}=0.03$, $N_{0}=1000$,
$\bar{\nu}_{n}^{(2)}=2$, $\bar{\nu}_{m}^{(2)}=1.2$. The two offspring
distributions were supposed independent so that
$\overline{\nu_{m}(\nu_{n}-1)}=\frac{\beta}{1-\beta}\bar{\nu}_{n}^{(2)}$.

##  V Stochastic diffusion equation

Now equipped with the main phenomena of 0 dimensional neutronics, our field
theoretic model can transition to proper neutron transport theory. A
convenient and easy way to do so, yet realistic 222the diffusion equation
accurately reproduces the transport equation for isotropic media and far from
the boundaries, see [42] for more details (it is used in the operational
context of reactor core monitoring), is to lean on the Brownian motion and its
associated diffusion equation. The diffusion equation, which can also formally
be built as a coarse grained description of a microscopic, stochastic process.
When supplemented by an external random noise term, the (one energy group)
equation can be written as

$\frac{\partial}{\partial
t}N(\vec{x},t)-D(\Delta+r)N(\vec{x},t)+\eta(\vec{x},t)=0,$ (46)

where the rescaled parameter $r=\rho/D$ as been defined for convenience and
where now $N(\vec{x},t)$ should be interpreted as the neutron density field.
The transcription of this $d+1$ diffusion equation into a response functional
follows, here also, exactly the same line as in section II. One thus have (in
$d$ spatial dimensions),

$\mathcal{S}=\int d^{d}\vec{x}\int
dt\widetilde{N}(\vec{x},t)\left(\partial_{t}-D(\Delta+r)\right)N(\vec{x},t)-\widetilde{N}\Gamma\widetilde{N}.$
(47)

The propagator of the theory is most conveniently written in the Fourier
transformed spatial domain ($\vec{p}$,$t$),

$\mathcal{G}(\vec{p},t)=e^{-D(\vec{p}^{\,2}-r)t}.$ (48)

The structure of the interaction term of the response functional is not
affected by the dimensionality of the problem. As a consequence, the
diagrammatic expansion of the theory is the same as the one obtained in
section III. It can for instance be immediately concluded that provided with
the initial condition $N(\vec{x},t)=N_{0}\,\delta(t)\delta^{(d)}(\vec{x})$,
one finds that

$\langle
N(\vec{x},t)\rangle=N_{0}\mathcal{G}(\vec{x},t)=N_{0}\,\sqrt{\frac{\pi}{4Dt}}\,e^{rDt-\frac{|\vec{x}|^{2}}{4Dt}}.$
(49)

A substantial amount of work has been done in the past, revolving around the
possibility that neutrons evolving in fissile matter might gather into
clusters and that the neutron density map should consequently present
patchiness. This neutron clustering phenomenon results from the very nature of
the underlying branching process they are subjected to [7, 8, 29, Agarwal2001,
9, 30, 15], and more specifically from the asymmetry between their uniform
capture in the medium and their production only at spatial sites correlated to
their parents Ṫhis phenomenon can be quantified through the equal-time spatial
correlation function $\langle N(\vec{x},t)N(\vec{y},t)\rangle$ which, in the
response functional formalism can be written as

$\displaystyle\langle
N(\vec{x},t)N(\vec{y},t)\rangle_{c}=2\lambda_{1}\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[above left=of a](f1); \vertex[below left=of a](f2); \vertex[right=of
a](o1); \diagram* { (f1) -- [anti fermion] (a) -- [fermion] (f2), (a) -- [anti
fermion] (o1) }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}=2\lambda_{1}\int
d^{d}\vec{w}\int_{0}^{t}dt_{1}\mathcal{G}(\vec{w},t_{1})\mathcal{G}(\vec{x}-\vec{w},t-t_{1})\mathcal{G}(\vec{y}-\vec{w},t^{\prime}-t_{1})$
$\displaystyle=2\lambda_{1}\int_{0}^{t}dt_{1}\iiint\frac{d^{d}\vec{p_{1}}}{(2\pi)^{d}}\frac{d^{d}\vec{p_{2}}}{(2\pi)^{d}}\frac{d^{d}\vec{p_{3}}}{(2\pi)^{d}}\mathcal{G}(\vec{p_{1}},t_{1})\mathcal{G}(\vec{p_{2}},t-t_{1})\mathcal{G}(\vec{p_{3}},t-t_{1})e^{-i\vec{p_{2}}\cdot\vec{x}}e^{-i\vec{p_{3}}\cdot\vec{y}}\int
d^{d}\vec{w}e^{-i\vec{w}\cdot(\vec{p_{1}}-\vec{p_{2}}-\vec{p_{3}})}$
$\displaystyle=2\lambda_{1}\int_{0}^{t}dt_{1}\iiint\frac{d^{d}\vec{p_{1}}}{(2\pi)^{d}}\frac{d^{d}\vec{p_{2}}}{(2\pi)^{d}}\frac{d^{d}\vec{p_{3}}}{(2\pi)^{d}}\mathcal{G}(\vec{p_{1}},t_{1})\mathcal{G}(\vec{p_{2}},t-t_{1})\mathcal{G}(\vec{p_{3}},t-t_{1})e^{-i\vec{p_{2}}\cdot\vec{x}}e^{-i\vec{p_{3}}\cdot\vec{y}}(2\pi)^{d}\delta^{(d)}(\vec{p_{1}}-\vec{p_{2}}-\vec{p_{3}})$
$\displaystyle=2\lambda_{1}\int_{0}^{t}dt_{1}\iint\frac{d^{d}\vec{p_{2}}}{(2\pi)^{d}}\frac{d^{d}\vec{p_{3}}}{(2\pi)^{d}}\mathcal{G}(\vec{p_{2}}+\vec{p_{3}},t_{1})\mathcal{G}(\vec{p_{2}},t-t_{1})\mathcal{G}(\vec{p_{3}},t-t_{1})e^{-i\vec{p_{2}}\cdot\vec{x}}e^{-i\vec{p_{3}}\cdot\vec{y}}.$
(50)

Note that because of the appearance of Dirac delta functions enforcing
momentum conservation at each vertex, it becomes useful to indicate on the
Feynman graphs the directions of the flowing momenta.

To proceed further with the calculation, one has to pay special attention to
the initial condition chosen for the neutron population. In the case of
neutrons evolving in an homogeneous and infinite system, and setting a uniform
initial neutron density $N_{0}$, one can use the translation invariance of the
problem to set $\vec{p_{1}}=\vec{0}$. Momentum conservation then implies that
Eq. V be transformed into

$\displaystyle\langle N(\vec{x},t)N(\vec{y},t)$
$\displaystyle\rangle_{c}=2\lambda_{1}N_{0}\int_{0}^{t}dt_{1}\int\frac{d^{d}\vec{p}}{(2\pi)^{d}}\mathcal{G}(\vec{0},t_{1})\times$
(51)
$\displaystyle\mathcal{G}(\vec{p},t-t_{1})\mathcal{G}(-\vec{p},t-t_{1})e^{-i\vec{p}\cdot(\vec{x}-\vec{y})}.$

Performing the integration over the momentum $\vec{p}$ in $d$-dimension, one
obtains,

$\langle
N(\vec{x},t)N(\vec{y},t)\rangle_{c}=2\lambda_{1}N_{0}e^{Drt}\int_{0}^{t}\frac{e^{Drt^{\prime}}}{(8\pi
Dt^{\prime})^{d/2}}\,e^{-\frac{|\vec{x}-\vec{y}|^{2}}{8Dt^{\prime}}}dt^{\prime},$
(52)

in agreement with the literature [31].

##  VI From neutronics to reactor physics - addition of a simple feedback
mechanism

Reactor physics is neutron transport theory supplemented by feedback
mechanisms, which happen whenever the reactor is operated at high power. These
feedback mechanisms are a central aspect of the stability of nuclear reactors
but their detailed modelling is generally difficult to implement in a neutron
balance, principally because they cannot be associated to a microscopic
stochastic process: for example, the Doppler effect is proportional to local
neutron densities, since neutrons heat up their surrounding media. In the
response functional formalism, these effect can however be directly
incorporated at the macroscopic level of description, i.e. in the mean field
equations.

In the following, this advantage is therefore exploited to explore a
simplified but realistic model of a reactor in the presence of counter-
reactions. It will be shown that the feedback mechanism imposed on the neutron
population has a strong impact on the structure of the fluctuations close to
the percolation transition. The field theoretic viewpoint adopted here will in
particular allow to establish that the neutron gas undergoes, whenever the
reactor is exactly critical, a phase transition whose universality class will
be identified.

###  VI.1 A simplified model with Doppler broadening

One of the most paradigmatic feedback effect affecting the neutron population
in a reactor core is the Doppler broadening of the capture resonances of 238U,
which couples the neutron population to the temperature of the fuel. As the
temperature increases, more and more neutrons are captured by the medium and
disappear: the overall neutron population is pulled back towards a reference
operating power. This feedback is quite difficult to include in a detailed
microscopic model because it involves the reactor power or the temperature of
the medium which are macroscopic quantities. It can however be
straightforwardly implemented under the form of a restoring force of the type
$\sigma N(N-N_{\text{ref.}})$ in the diffusion equation [4]. The behavior of
the neutron population evolving in a nuclear reactor operating at nominal
power (and thus subject to Doppler broadening) can consequently be described
by the equation

$\frac{\partial}{\partial t}N(\vec{x},t)-D(\Delta+r)N(\vec{x},t)+\sigma
N^{2}(\vec{x},t)=0,$ (53)

where the shifted parameter $r=(\rho+\sigma N_{\text{ref.}})/D$ has been
introduced, for convenience.

The response functional and the resulting theory to which this equation is
associated has been extensively studied. It corresponds to the functional of
the Directed Percolation (DP) universality class [24, 27]:

$\mathcal{S}=\int d^{d}\vec{x}\int
dt\widetilde{N}\left(\partial_{t}-D(\Delta+r)\right)N+\sigma\widetilde{N}N^{2}-\lambda\widetilde{N}^{2}N.$
(54)

The inclusion of the term $\sigma N^{2}$ has a strong impact on the structure
of the theory, for now the two following processes are allowed

$2\lambda\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[above left=of a](f1); \vertex[below left=of a](f2); \vertex[right=of
a](o1); \diagram* { (f1) --[anti fermion] (a) --[fermion] (f2), (a) --[anti
fermion] (o1) }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\hskip
56.9055pt2\sigma\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(a);
\vertex[above right=of a](f1); \vertex[below right=of a](f2); \vertex[left=of
a](o1); \diagram* { (f1) --[fermion] (a) --[anti fermion] (f2), (a)
--[fermion] (o1) }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}$

as well as their combination, i.e. the loop diagram

$4\sigma\lambda\,\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman[inline=(a)] \vertex(z);
\vertex[right=of z](a); \vertex[right=of a] (b); \vertex[above left=of a](f1);
\vertex[below left=of a](f2); \vertex[right=of b](o1); \diagram* { (z) --[anti
fermion] (a) \-- [fermion,half left] (b) \-- [fermion,half left] (a), (b)
--[anti fermion] (o1) }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\,,$ (55)

whose contribution leads to badly divergent behaviors of two different types
[27]:

  * •

a divergence is observed in dimension $d\geq 2$, because the momentum
associated to the loop is not constrained by any conservation rule. Its
integration thus runs over theoretically infinite momenta and include small
distances, short time, unbounded fluctuations of the neutron population. This
is a paradigmatic example of an Ultra-Violet (UV) divergence;

  * •

another and more serious divergence occurs in dimension $d\leq 4$. An Infrared
(IR) divergence is observed in the long time, long distance limit. The highest
dimension $d_{c}=4$ for which the divergence is observed is called the
critical dimension of the problem.

Both types of divergences can be eliminated through the well-known
renormalization procedure [32], which, in essence, consists in the application
of successive spatio-temporal scaling transformations on the problem, so that
the small distance and UV divergent fluctuations be truly integrated over. The
true macro-scale description of the problem is obtained at a fixed point of
the iterated scaling process which generally has the induced benefit of taming
the residual IR divergence of the problem.

Although the starting point of the model are the mean field equations of the
problem, the renormalized observables one ultimately extract from the theory
present asymptotic behaviors (i.e. critical exponents) that differ from the
predictions of the mean field equations, at least in the vicinity of the phase
transition. Close to the critical point, mean field theory breaks down and
only the renormalized field theoretic description of the system can correctly
describe the evolution of observables.

###  VI.2 The directed percolation transition

The largest deviation from the mean field predictions is obtained in the 0D
model. Indeed, setting $D=0$ in Eq. 53, the point-like DP model can be defined
as

$\frac{dN}{dt}=rN-\sigma N^{2}.$ (56)

As long as $r>0$, this equation possesses a positive stationary solution
$N_{\infty}\propto r^{\beta}$ ($\beta=1$ is the first critical exponent of
directed percolation). A detailed analysis with the help of the RG leads to a
strikingly different prediction [33, 34]: in zero dimension, whatever the
values chosen for $r$ and $\sigma$, the neutron population always
asymptotically vanish. This finding is in sharp contrast with the solution of
the mean field equations but has been tested and validated against Monte-Carlo
simulations. Translated into the domain of nuclear reactor physics, this means
that adding a feedback term of the form $\sigma N^{2}$ into the point kinetics
equation leads to a theory that do not possess a stable, non-trivial,
equilibrium solution: every point-like reactors subject to counter-reactions
are sub-critical !

The triviality of this point-like model can be easily circumvented though, by
considering instead the case of a tri-dimensional, infinite and homogeneous
reactor. The system then do admit $N_{\infty}\propto r^{\beta}$ as a solution,
but the RG predicts a value for the exponent $\beta$ close to 0.8 in 3$d$ [27,
24] (as compared to 1 for mean field). Although conceptually important, this
deviation in the scaling relation between the power level and the reactivity
of the reactor would be difficult to observe in practice.

Even though the feedback mechanism proposed in this paper might appear as a
simplistic approximation of the true processes at play in an operating
reactor, the asymptotic scaling laws and the overall behaviors predicted by
the RG analysis of the directed percolation process might still persist in a
more refined treatment. The DP conjecture [35, 24] indeed states that the DP
universality class is the most general universality class observed in
transitions between an active to an absorbing state and persist under very
general conditions. The modelling of critical nuclear reactors might thus fall
into the general class of the DP percolation transitions, even with a more
elaborate modelling of counter-reaction mechanisms.

Observable consequences of the DP transition and direct violation of the mean
field equation predictions might however be difficult to observe in practice.
Apart from the exponent $\beta$, two other critical exponents characterize the
DP universality class. They are associated to the correlation length $\xi$ and
the diffusion time scale $\tau$ that naturally appear in equation 53:

$\xi\sim|r|^{-\nu}\;\;\text{ and }\;\;\tau\sim\xi^{z}/D\sim|r|^{-z\nu}.$ (57)

The DP class is entirely characterized by these exponents and any observable
consequence of the DP transition can be traced back to a deviation of these
exponents from their mean field predictions. Deviations from the mean field
theory are expected to be larger for system in lower dimensions, in models
farther away from the critical dimension $d_{c}=4$. Table 1 list the values of
the exponents of directed percolation, extracted from numerical simulations,
for systems in 1, 2 or 3 dimensions.

Table 1: Critical exponents of the DP phase transition [36].  exponent | d = 1 | d = 2 | d = 3 | mean field
---|---|---|---|---
$\beta$ | 0.276486(8) | 0.5834(30) | 0.813(9) | 1
$\nu$ | 1.096854(4) | 0.7333(75) | 0.584(5) | 0.5
$z$ | 1.580745(10) | 1.7660(16) | 1.901(5) | 2

A detailed account of all of the possible experimental evidences of the DP
transition in nuclear reactors would by far exceed the scope of this paper. As
a matter of illustration, the following section will describe what observable
and concrete consequence the dynamic exponent $z$ might have in a real
reactor.

###  VI.3 A Consequence on reactor physics

Reactor physics aims at accurately calculating (with a precision of $\sim
1\%$) the spatial distribution of the neutron flux, from which diverse
measurable local quantities can be extracted, such as fission rates, neutron
counts, etc. Those quantities are measured either by in-core detection devices
(e.g. the aeroball measurement system) or by ex-core detection devices (using
source-range, intermediate-range or power-range neutron flux
instrumentations). Any significant deviation from the flux measured by these
detectors to ’expected’ values should trigger various safety warnings notified
by the automatic protection system, ultimately leading to the core tripping
procedure. Noticeably these expected values of the flux are set by an online
solving of a simplified -yet robust- neutron transport model.

During the intermediate and power range operation of the reactor, the thermal
feedback effects have to be taken into account as they modify the cross-
sections dictating the laws of propagation of the neutrons. As seen in the
previous section, this modification has another consequence: the macroscopic
properties of the medium emerge from an intimate coupling at the mesoscopic
scale between the random walk of the neutrons and the feedback reactions
depending on the neutrons local densities. Usual coupling schemes relying on
mean field equations are therefore not valid anymore since the mean field
equations themselves are wrong, and the properties of the system can only be
obtained through the critical exponents of the directed percolation defined in
the previous section. In the following, since we are interested in a
simplified model of reactor physics, we will restrict ourselves to $d=3$. In
this case, the value of $z$, the dynamic critical exponent, is about $1.9$.
This exponent describes how fast any local perturbation spreads [37] and can
be related to a modification of the diffusive ($z=2$) random walk which
becomes slightly super-diffusive for $d=3$ (indeed $z$ can be interpreted as
an anisotropy parameter). The scaling of the spatial correlation length versus
the time correlation length follows $\xi^{z}\sim D\tau$ (see Eq. 57, with
$z=1.901$ in the following). These Lévy flights (see for instance [38] for a
review on fractional Laplacians and their connection to Lévy flights) call for
a more quantitative modelling of the mean-field Eq. 53 which is hence modified
according to:

$D(-(-\Delta)^{z/2}+r)N(x,t)-\sigma N^{2}(x,t)=\frac{\partial}{\partial
t}N(x,t),$ (58)

where $-(-\Delta)^{z/2}$ is the fractional Laplacian defined through
$-(-\Delta)^{z/2}=\frac{\partial^{z}}{\partial|x|^{z}}$
($\frac{\partial^{z}}{\partial|x|^{z}}$ being the Riesz-Feller derivative of
fractional order $z$). Yet, whenever the reactor is in a stationary state, the
Doppler effect does hardly comes into action: it mainly stabilizes the core
upon local perturbations or control long-term statistical fluctuations [15].
Since the critical exponents obtained in the previous section did not
explicitly depend on $\sigma$, representing the strength of the coupling
between the neutron scalar field and the feedback mechanism, this equation can
further be simplified by considering the parametric regime where this coupling
is weak. The feedback term $\sigma N^{2}$ in the mean field equation can
consequently be safely neglected. In this case the stationary behavior of Eq.
58 boils down to

$D(-\Delta)^{z/2}N(x,t)=rN(x,t).$ (59)

This critical diffusion equation essentially tells that, in the frame of the
so-called one-group diffusion equation, moving from neutron transport theory
(without feedback) to reactor physics theory (i.e. neutronics with feedback)
mainly consists in replacing the standard Laplacian ($z=2$) by the fractional
Laplacian with $z=1.901$. The equilibrium spatial distribution of the neutron
field can therefore be obtained through the analysis of the first eigenmode of
this fractional Laplacian. Since the reactor has a finite spatial extent, it
is necessary to estimate this eigenmode in a bounded domain, keeping in mind
that the critical exponents obtained through the use of the renormalization
group are valid only for an infinite medium. To avoid the intricate analysis
of critical exponents of the directed percolation in finite domain, we will
anyway use this (infinite medium) dynamic critical exponent while keeping in
mind that the results of the analysis will be only valid within the bulk of
the medium (i.e. far from the boundaries). The precise topic of the properties
of the fractional Laplacian in bounded domains has been investigated in the
seminal paper of Zoia et al [39]. This paper proposes different way to
numerically tackle such equations. One possible -and elegant- way being to
perform a Monte Carlo simulation of the neutron transport where the standard
diffusion kernel is replaced by Lévy flights. Given the different hypotheses
of our reactor model (diffusion instead of transport, no energy dependence,
uniform medium, restriction to the bulk of the domain,…) we are only
interested in approximate solutions, and therefore another possibility is to
have recourse to approximate analytical expression as the one given by [40]
for example. Figure 2 reports the comparison between the exact solution of the
critical diffusion equation that models neutron transport and such an
approximate solution of our critical super-diffusion equation that models
reactor physics. This approximate solution has the form

$N(x)=A\cos\Bigl{(}\left(\frac{\pi}{2}-(2-z)\frac{\pi}{8}\right)\frac{x}{L}\Bigr{)}$
(60)

with A being a normalization constant, and is obtained for a $1d$ uniform
medium by imposing Dirichlet boundary conditions ($N(\pm L)=0$) at the
boundaries $L$ of the domain for $z=2$. For $z<2$, the neutron density at the
boundaries does not vanish: this can be interpreted as long-distance flight
bringing some neutrons to cross the boundary instead of being absorbed at it.
But in the present case we are not interested in the properties of the medium
close to the boundaries (where the dynamic critical coefficient is not valid
anymore), hence the proper boundaries conditions for $z=1.901$ are not
discussed in this paper. It is however worth noting that the topic of
boundaries conditions for Lévy flights and their physics interpretation is
currently under active investigation (see [41] for instance). Equation Eq. 60
can also be extended to $3d$ medium by multiplication of the modes along each
axis since $N(x,y,z)=N(x)N(y)N(z)$ (this separation of variables stands for
homogeneous partial differential equations in general [42]), making the study
of a $1d$ medium meaningful. Also, for $z$ values close to $2$ the
approximated solution comes close to being exact (the $L^{2}$ norm between the
exact and the approximate solution is controled by $(2-z)$ [40]). When $z=2$,
this approximate solution becomes exact since one retrieves the usual solution
of the critical diffusion equation
$N(x)\sim\cos\Bigl{(}\frac{\pi}{2}\frac{x}{L}\Bigr{)}$. The $1d$ comparison
between the spatial distributions of $z=2$ and $z=1.901$ (top plot of Figure
2) finally reveals that, as could be qualitatively guessed, the super-
diffusion slightly flatten the spatial distribution compared to diffusion.
This macroscopic point of view emerges from the fact that the far-distance
events characteristics of Lévy flights indeed tend to spatially redistribute
the neutrons. In the frame of all the aforementioned hypotheses, the
differences between the diffusion and the $z=1.901$ super-diffusion on local
neutron densities are of the order of $1-2\%$ (bottom plot of Figure 2), which
is precisely the accuracy targeted by high precision reactor physics codes and
methods, thus making these relative differences non-negligible.

Figure 2:  Top plot: comparison between the first normalized eigenfunction of
the linear diffusion equation ($z=2$ in plain line) with Dirichlet boundary
conditions in $-1$ and $1$, and an approximation of the first normalized
eigenfunction of the fractional Laplacian (as given by [40]) with $z=1.901$
(dashed line). The approximation is valid far from the boundaries, i.e. within
the bulk. Bottom plot: relative differences (in %) between both
eigenfunctions.

##  VII Conclusion

The response functional formalism of Janssen and De Dominicis has proven to be
a versatile and powerful companion for exploring the physics of branching
processes. Applied to the well known equations of neutron transport in fissile
media (neutronics), we have been able to recover with relative ease the main
recent results of the field : moments of the neutron population within the
reactor, effect of external neutron sources on the chain reaction, taking into
account of delayed neutrons, and characterization of spatial and temporal
correlations that may develop at low power. Along this path, the close
connection between the statistical field theory of Janssen and De Dominics and
the coherent path integral method of Doi and Peliti has been underlined. This
last approach is based on establishing a second quantized version of the
neutron’s branching process, and therefore relies on the solid-ground
microscopic description of branching spatial Brownian motion that has been
build in the past few decades. But the generality of the response functional
approach lend itself to the inclusion of very general ad-hoc models whose
microscopic origin would be difficult to implement. The restoring force that
has been used in section VI is an example of such a heuristic approach. It
allowed the study of a simplistic model of nuclear reactor, that couples
neutronics to a prototype feedback model (the Doppler effect), which is
usually not possible via the traditional master equations approach. Showing
that the phase transition of the neutron population at criticality falls into
the directed percolation universality class permitted to retrieve the critical
exponents of the phase transition obtained thanks to the well established
results of the renormalization group technics.

The main consequence of this work in the realm of nuclear reactor physics were
evaluated by focussing on the corrected mean field equation of our simplified
reactor physics model, using in particular the dynamic critical exponent $z$
previously obtained. In the same manner that this coefficient transforms the
diffusion law $z=2$ into a super-diffusive Lévy flights with $z=1.901$ for
$3d$ medium, the usual critical diffusion equation used to grasp the main
features of neutronics was shown to be replaced, in this prototype model of
reactor physics, by a critical super-diffusion equation. In the low feedback
regime, and far from the boundaries, the main eigenmode of this critical
super-diffusion equation was approximated using recent developments on the
spectrum of the fractional Laplacian in bounded domains. In this frame, the
solution of the standard mean field equation for the neutron local
distributions might show deviations up to $1-2\%$, which invalidates many
formal approaches and reactor physics analytical models (see [43, 44] among
many others) sometimes used to benchmark reactor physics codes, while
questioning the accuracy targeted by the current high precision reactor
physics codes themselves, which aim at calculating neutron fluxes with an
accuracy below this precision (see [45] for instance).

Finally, to further assess the accuracy of mean field approaches (and hence
reactor physics codes), it would be necessary to investigate the properties of
the directed percolation in finite medium, while studying the effects of
neutron energy (in the frame of a 2-group critical super-diffusion equation
for instance) and improving the modeling of the feedback mechanism. A
prominent role might also be played by the precursor of delayed neutrons.
Indeed, these precursors being not subject to any feedback mechanism, they
might inflect the system’s dynamic into a different universality class.
Ultimately, reactor physics experiments aiming at measuring and detecting
distortions of the neutron field could definitely sign this phenomenon of
neutron percolation in reactor cores.

## Acknowledgments

The authors are thankful to Léonie Canet and Andrea Zoia for insightful
discussions and comments.

##  Appendix A Analytical calculation of the generalized propagators
$\Sigma_{ab}$ ($a,b$ = $N,M$) in the two species model

Focusing on the evaluation of $\Sigma_{NN}(t,t^{\prime})$, and using the
shorthand notation $a=r+\lambda_{D}$, the diagram of Eq. 34 translates into
the infinite series,

$\Sigma_{NN}(t,t^{\prime})=e^{r(t-t^{\prime})}\left\\{1+\beta\lambda_{D}\int_{t^{\prime}}^{t}dt_{1}e^{-at_{1}}\int_{t^{\prime}}^{t_{1}}dt_{2}e^{at_{2}}+(\beta\lambda_{D})^{2}\int_{t^{\prime}}^{t}dt_{1}e^{-at_{1}}\int_{t^{\prime}}^{t_{1}}dt_{2}e^{at_{2}}\int_{t^{\prime}}^{t_{2}}dt_{3}e^{-at_{3}}\int_{t^{\prime}}^{t_{3}}dt_{4}e^{at_{4}}+...\right\\},$
(61)

or,

$\Sigma_{NN}(t,t^{\prime})=e^{r(t-t^{\prime})}\left(1+\sum_{n\geq
1}(\beta\lambda_{D})^{n}A_{n}(t,t^{\prime})\right),$ (62)

with,

$\displaystyle A_{n}(t,t^{\prime})=$
$\displaystyle\int_{t^{\prime}}^{t}\int_{t^{\prime}}^{t_{1}}...\int_{t^{\prime}}^{t_{2n-1}}dt_{1}...dt_{2n}e^{-a(t_{1}-t_{2}+...+t_{2n-1}-t_{2n})}$
(63)

One can first easily show that $A_{n}(t,t^{\prime})=A_{n}(t-t^{\prime},0)$, so
that only $A_{n}(t):=A_{n}(t,0)$ needs to be evaluated. Then, the form of the
exponential calls for a change of variable of the form
$t_{1},...,t_{2n}\rightarrow
t_{1}-t_{2},...,t_{2n-1}-t_{2n},t_{2n}=u_{1},...,u_{2n}$. Using this new set
of variables, $A_{n}(t)$ can be written,

$A_{n}(t)=\idotsint\limits_{\begin{subarray}{c}0<u_{1},\ldots,u_{2n}\\\
u_{1}+..+u_{2n}<t\end{subarray}}du_{1}\ldots
du_{2n}e^{-a(u_{1}+\ldots+u_{2n-1})},$ (64)

where the exponential now only depends on the odd terms $u_{2i+1}$. The
constraint $\sum_{i=1}^{2n}u_{i}<t$ can be enforced through the introduction
of an indicator function $\mathbbm{1}_{\sum_{i=1}^{2n}u_{i}<t}$, so that,

$\displaystyle A_{n}(t)=$
$\displaystyle\idotsint\limits_{0<u_{1},\ldots,u_{2n}}du_{1}\ldots
du_{2n}e^{-a(u_{1}+\ldots+u_{2n-1})}\mathbbm{1}_{\sum_{i=1}^{2n}u_{i}<t}$ (65)

The integration over even indices corresponds now to a volume integral over an
$n$-dimensional simplex,

$\idotsint\limits_{0<x_{1},\ldots,x_{p}}dx_{1}\ldots
dx_{p}\mathbbm{1}_{\sum_{i=1}^{p}x_{i}<\alpha}=\frac{\alpha^{p}}{p!}.$ (66)

Splitting Eq. A over even/odd indices, one can then show that,

$\begin{split}A_{n}(t)&=\idotsint\limits_{0<u_{1},\ldots,u_{2n-1}}du_{1}\ldots
du_{2n-1}\mathbbm{1}_{\sum_{i=0}^{n-1}u_{2i+1}<t}e^{-a(u_{1}+\ldots+u_{2n-1})}\idotsint\limits_{0<u_{2},\ldots,u_{2n}}du_{2}\ldots
du_{2n}\mathbbm{1}_{\sum_{i=1}^{n}u_{2i}<t-\sum_{i=0}^{n-1}u_{2i+1}}\\\
&=\idotsint\limits_{0<u_{1},\ldots,u_{2n-1}}du_{1}\ldots
du_{2n-1}\mathbbm{1}_{\sum_{i=0}^{n-1}u_{2i+1}<t}e^{-a(u_{1}+\ldots+u_{2n-1})}\frac{\left(t-u_{1}-\ldots
u_{2n-1}\right)^{n}}{n!}\end{split}$ (67)

The integrand in the last expression only depends on the sum
$u_{1}+\ldots+u_{2n-1}$. Thus, defining $x=\sum\limits_{i=0}^{n-1}u_{2i+1}$
and changing the integration variables accordingly, one ends up with,

$A_{n}(t)=\int\limits_{0}^{t}dxe^{-ax}\frac{\left(t-x\right)^{n}}{n!}\frac{x^{n-1}}{(n-1)!}.$
(68)

$\Sigma_{NN}$ can thus be written,

$\displaystyle\Sigma_{NN}(t,0)$ $\displaystyle=$ $\displaystyle
e^{rt}\Bigg{\\{}1+$ (69) $\displaystyle\left.\sum\limits_{n\geq
0}(\beta\lambda_{D})^{n+1}\int_{0}^{t}dx~{}e^{-ax}\frac{\left(t-x\right)^{n+1}x^{n}}{(n+1)!~{}n!}\right\\}.$

The latter expression involves the series development of the Modified Bessel
function of the first kind $I_{1}$. Reinstating the lower bound $t^{\prime}$
into the equation, one then has,

$\displaystyle\Sigma_{NN}(t,t^{\prime})$
$\displaystyle=e^{r(t-t^{\prime})}\Bigg{\\{}1+$
$\displaystyle\int_{0}^{t-t^{\prime}}d\tau
e^{-r\tau}\sqrt{\frac{\beta\lambda_{D}(t-t^{\prime}-\tau)}{\tau}}$
$\displaystyle\times$ $\displaystyle
I_{1}\left(2\sqrt{\beta\lambda_{D}\tau(t-t^{\prime}-\tau)}\right)\Bigg{\\}}.$
(70)

One the same model, one can show that,

$\displaystyle\Sigma_{MM}(t,t^{\prime})$
$\displaystyle=e^{-\lambda_{D}(t-t^{\prime})}\Bigg{\\{}1+$ (71)
$\displaystyle\int_{0}^{t-t^{\prime}}d\tau
e^{a\tau}\sqrt{\frac{\beta\lambda_{D}(t-t^{\prime}-\tau)}{\tau}}$
$\displaystyle\times
I_{1}\left(2\sqrt{\beta\lambda_{D}\tau(t-t^{\prime}-\tau)}\right)\Bigg{\\}},$
$\displaystyle\Sigma_{NM}(t,t^{\prime})$ $\displaystyle=\beta
e^{-\lambda_{d}(t-t^{\prime})}$ (72)
$\displaystyle\times\int_{0}^{t-t^{\prime}}d\tau~{}e^{a\tau}I_{0}\left(2\sqrt{\beta\lambda_{D}\tau(t-t^{\prime}-\tau)}\right),$

and

$\displaystyle\Sigma_{MN}(t,t^{\prime})$
$\displaystyle=\lambda_{D}~{}e^{r(t-t^{\prime})}$ (73)
$\displaystyle\times\int_{0}^{t-t^{\prime}}d\tau
e^{-a\tau}I_{0}\left(2\sqrt{\beta\lambda_{D}\tau(t-t^{\prime}-\tau)}\right).$

All these expressions contain integrals that can be further evaluated. After
successive change of variables, and setting
$\Omega=4\beta\lambda_{D}+\lambda_{D}^{2}+2\lambda_{D}r+r^{2}$, one can show
[46] that,

$\displaystyle\Sigma_{NN}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\frac{e^{\frac{1}{2}(r-\lambda_{D})(t-t^{\prime})}}{\sqrt{\Omega}}\left\\{\sqrt{\Omega}\cosh\left(\frac{1}{2}\sqrt{\Omega}(t-t^{\prime})\right)\right.$
(74)
$\displaystyle\left.+(r+\lambda_{D})\sinh\left(\frac{1}{2}\sqrt{\Omega}(t-t^{\prime})\right)\right\\},$

$\displaystyle\Sigma_{MM}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\frac{e^{\frac{1}{2}(r-\lambda_{D})(t-t^{\prime})}}{\sqrt{\Omega}}\left\\{\sqrt{\Omega}\cosh\left(\frac{1}{2}\sqrt{\Omega}(t-t^{\prime})\right)\right.$
(75)
$\displaystyle\left.-(r+\lambda_{D})\sinh\left(\frac{1}{2}\sqrt{\Omega}(t-t^{\prime})\right)\right\\},$

$\Sigma_{NM}(t,t^{\prime})=2\beta\frac{e^{\frac{1}{2}(r-\lambda_{D})(t-t^{\prime})}}{\sqrt{\Omega}}\sinh\left(\frac{1}{2}\sqrt{\Omega}(t-t^{\prime})\right),$
(76)

$\Sigma_{MN}(t,t^{\prime})=2\lambda_{D}\frac{e^{\frac{1}{2}(r-\lambda_{D})(t-t^{\prime})}}{\sqrt{\Omega}}\sinh\left(\frac{1}{2}\sqrt{\Omega}(t-t^{\prime})\right).$
(77)

## References

  * Spanier and Gelbard [1969] J. Spanier and E. M. Gelbard,  _Monte Carlo Principles and Neutron Transport Problems_ (Dover Publications, 1969).
  * Lux and Koblinger [1991] I. Lux and L. Koblinger,  _Monte Carlo particle transport methods : neutron and photon calculations_ (CRC Press, 1991).
  * Note [1] Neutron transport theory, also known as neutronics, aims at characterizing the properties of the neutron gas in fissile media.
  * Bell and Glasstone [1970] G. I. Bell and S. Glasstone,  _Nuclear Reactor Theory_ (Krieger Publishing Company, 1970).
  * Williams [1974] M. Williams,  _Random Processes in Nuclear Reactors_ (Pergamon, 1974).
  * Pázsit and Pál [2007] I. Pázsit and L. Pál,  _Neutron Fluctuations: A Treatise on the Physics of Branching Processes_ (Elsevier Science Ltd, 2007).
  * Dumonteil  _et al._ [2014] E. Dumonteil, F. Malvagi, A. Zoia, A. Mazzolo, D. Artusio, C. Dieudonné, and C. De Mulatier, Particle clustering in Monte Carlo criticality simulations, Annals of Nuclear Energy 63, 612 (2014).
  * Zoia  _et al._ [2014] A. Zoia, E. Dumonteil, A. Mazzolo, C. De Mulatier, and A. Rosso, Clustering of branching Brownian motions in confined geometries, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 90, 042118 (2014).
  * Sutton and Mittal [2017] T. M. Sutton and A. Mittal, Neutron clustering in Monte Carlo iterated-source calculations, Nuclear Engineering and Technology 49, 1211 (2017).
  * Harris [1963] T. E. Harris,  _The theory of branching processes_ (Springer Berlin, Heidelberg, 1963).
  * Zoia  _et al._ [2012] A. Zoia, E. Dumonteil, A. Mazzolo, and S. Mohamed, Branching exponential flights: Travelled lengths and collision statistics, J. Phys. A: Math. Theor. 45, 425002 (2012).
  * Young  _et al._ [2001] W. R. Young, A. J. Roberts, and G. Stuhne, Reproductive pair correlations and the clustering of organisms, Nature 2001 412:6844 412, 328 (2001).
  * Houchmandzadeh [2002] B. Houchmandzadeh, Clustering of diffusing organisms, Physical Review E 66, 052902 (2002).
  * Houchmandzadeh [2008] B. Houchmandzadeh, Neutral clustering in a simple experimental ecological community, Physical Review Letters 101, 078103 (2008).
  * Dumonteil  _et al._ [2021] E. Dumonteil, R. Bahran, T. Cutler, B. Dechenaux, T. Grove, J. Hutchinson, G. McKenzie, A. McSpaden, W. Monange, M. Nelson, N. Thompson, and A. Zoia, Patchy nuclear chain reactions, Communications Physics 2021 4:1 4, 1 (2021).
  * Janssen [1976] H. Janssen, On a lagrangean for classical field dynamics and renormalization group calculations of dynamical critical properties, Z Physik B 23, 378 (1976).
  * Dominicis [1976] C. D. Dominicis, Technics of field renormalization and dynamics of critical phenomena, J. Phys. (Paris) 1 (1976).
  * Doi [1976] M. Doi, Second quantization representation for classical many-particle system, Journal of Physics A: Mathematical and General 9, 1465 (1976).
  * Peliti [1985] L. Peliti, Path integral approach to birth-death processes on a lattice, J. Phys. France 46, 1469 (1985).
  * Cardy and Täuber [1996] J. Cardy and U. C. Täuber, Theory of branching and annihilating random walks, Phys. Rev. Lett. 77, 4780 (1996).
  * Martin  _et al._ [1973] P. C. Martin, E. D. Siggia, and H. A. Rose, Statistical dynamics of classical systems, Phys. Rev. A 8, 423 (1973).
  * Greiner [1996] W. Greiner,  _Field Quantization_ (Springer-Verlag Berlin and Heidelberg GmbH & Co. K, 1996).
  * Zinn-Justin [2004] J. Zinn-Justin,  _Path integrals in Quantum Mechanics_ (Oxford University Press, 2004).
  * Janssen, H. K. and Taüber, U. C. [2005] Janssen, H. K. and Taüber, U. C., The field theory approach to percolation processes, Annals Phys. 315, 147 (2005), arXiv:cond-mat/0409670 .
  * Garcia-Millan  _et al._ [2018] R. Garcia-Millan, J. Pausch, B. Walter, and G. Pruessner, Field-theoretic approach to the universality of branching processes, Phys. Rev. E 98, 062107 (2018).
  * Baez and Biamonte [2012] J. C. Baez and J. Biamonte, Quantum Techniques for Stochastic Mechanics, arXiv e-prints , arXiv:1209.3632 (2012), arXiv:1209.3632 [quant-ph] .
  * Taüber  _et al._ [2005] U. C. Taüber, M. Howard, and B. P. Vollmayr-Lee, Applications of field-theoretic renormalization group methods to reaction–diffusion problems, Journal of Physics A: Mathematical and General 38, R79 (2005).
  * Note [2] The diffusion equation accurately reproduces the transport equation for isotropic media and far from the boundaries, see [42] for more details.
  * De Mulatier  _et al._ [2015] C. De Mulatier, E. Dumonteil, A. Rosso, and A. Zoia, The critical catastrophe revisited, Journal of Statistical Mechanics: Theory and Experiment 2015, 10.1088/1742-5468/2015/08/P08021 (2015).
  * Dumonteil  _et al._ [2017] E. Dumonteil, G. Bruna, F. Malvagi, A. Onillon, and Y. Richet, Clustering and traveling waves in the Monte Carlo criticality simulation of decoupled and confined media, Nuclear Engineering and Technology 49, 1157 (2017).
  * de Mulatier [2015] C. de Mulatier, A random walk approach to stochastic neutron transport (2015).
  * Delamotte [2004] B. Delamotte, A hint of renormalization, American Journal of Physics 72, 170 (2004).
  * Canet and Hilhorst [2006] L. Canet and H. Hilhorst, Single-site approximation for reaction-diffusion processes., J Stat Phys 125, 517 (2006).
  * Benitez and Wschebor [2013] F. Benitez and N. Wschebor, Branching and annihilating random walks: Exact results at low branching rate, Phys. Rev. E 87, 052132 (2013).
  * Hinrichsen [2000] H. Hinrichsen, Non-equilibrium critical phenomena and phase transitions into absorbing states, Advances in Physics 49, 815 (2000).
  * Henkel  _et al._ [2008] M. Henkel, H. Hinrichsen, and S. Lübeck,  _Non-equilibrium Phase Transitions, Vol. 1: Absorbing Phase Transitions_ (Springer, 2008).
  * Täuber [2012] U. C. Täuber,  _Critical Dynamics: A Field Theory Approach to Equilibrium and Non-Equilibrium Scaling Behavior_ (Cambridge University Press, 2012) pp. 1–511.
  * Lischke  _et al._ [2020] A. Lischke, G. Pang, M. Gulian, F. Song, C. Glusa, X. Zheng, Z. Mao, W. Cai, M. M. Meerschaert, M. Ainsworth, and G. E. Karniadakis, What is the fractional Laplacian? A comparative review with new results, Journal of Computational Physics 404, 109009 (2020).
  * Zoia  _et al._ [2007] A. Zoia, A. Rosso, and M. Kardar, Fractional Laplacian in bounded domains, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 76, 10.1103/PHYSREVE.76.021116 (2007).
  * Kwaśnicki [2012] M. Kwaśnicki, Eigenvalues of the fractional Laplace operator in the interval, Journal of Functional Analysis 262, 2379 (2012).
  * Constantin and Ignatova [2016] P. Constantin and M. Ignatova, Remarks on the Fractional Laplacian with Dirichlet Boundary Conditions and Applications, International Mathematics Research Notices 2017, 1653 (2016), https://academic.oup.com/imrn/article-pdf/2017/6/1653/11066188/rnw098.pdf .
  * Duderstadt and Martin [1979] J. J. Duderstadt and W. R. Martin,  _Transport theory_ (John Wiley & Sons, New York, 1979) p. 613.
  * Dellacherie  _et al._ [2016] S. Dellacherie, E. Jamelot, and O. Lafitte, A simple monodimensional model coupling an enthalpy transport equation and a neutron diffusion equation, Applied Mathematics Letters 62, 35 (2016).
  * Tomatis  _et al._ [2021] D. Tomatis, R. Gross, and E. Gilad, On the ronen method in simple 1-d geometries for neutron transport theory solutions, Journal of Computational and Theoretical Transport 50, 134 (2021), https://doi.org/10.1080/23324309.2020.1843496 .
  * Demazière, Christophe  _et al._ [2020] Demazière, Christophe, Sanchez-Espinoza, Victor Hugo, and Chanaron, Bruno, Advanced numerical simulation and modelling for reactor safety - contributions from the cortex, hpmc, mcsafe and nuresafe projects, EPJ Nuclear Sci. Technol. 6, 42 (2020).
  * Gradshteyn  _et al._ [2014] I. S. Gradshteyn, I. M. Ryzhik, D. Zwillinger, and V. Moll,  _Table of integrals, series, and products; 8th ed._ (Academic Press, Amsterdam, 2014).