Modeling of neutron fluctuations with temperature
feedback using stochastic differential equations

B. Dechenaux

ªAutorité de Sûreté Nucléaire et de Radioprotection
(ASNR), PSN-RES/SNC/LN, Fontenay-aux-Roses, F-92260, France

Abstract

Neutron fluctuations, also coined zero power noise, is concerned with the study of the random fluctuations observed in the neutron population evolving in a fissile medium, due to the fundamentally intrinsic stochastic nature of the neutrons’ interactions in matter.

The canonical theoretical modeling of such fluctuations rely on the estab-
lishment of a microscopic probability balance equation, from which one can
extract, with the highest precision, all of the information about the stochastic
behavior of the neutron population.

An approximate yet practical scheme has been gaining in popularity in
the recent years. It relies on a coarse grained approach of the problem, where
the evolution of the neutron population is modeled through (Itô prepoint)
Stochastic Differential Equations (SDE). In this simplified scheme, the crux
is to find the SDE that best approximates the fluctuations of the underlying
microscopic stochastic process.

In the present work, an alternative approach to the scheme usually found
in the neutronics literature is proposed. A family of SDE that best approx-
imate the fine grain behavior of the neutron population is derived from an
explicit truncation of the translation operators found in the microscopic mas-
ter equation of the problem. The approach is developed over a concrete and
relatively original application case, namely that of the point reactor model
including mono-energetic neutrons, one group of precursors, and a tempera-
ture feedback process.

*Keywords:* Neutron fluctuations, zero power noise, stochastic differential equations, temperature feedback
# 1. Introduction

The characterization of the noise induced by the stochastic fluctuations associated with the random nature of neutron interactions in fissile media is an important concern in many different areas of neutronics, being for instance an important topic for the monitoring and measurement of the kinetic parameters in a nuclear reactor [1, 2]. The accurate modeling of these random fluctuations, observed on the macroscopic neutron population, relies fundamentally on a microscopic balance of the probabilistic events to which the neutrons are subject. The result of such a probability balance is referred to as master equation. It concretely comes in the form of a partial differential equation for the probability function of finding *N* neutrons (and *C* precursors, etc.) in the system at time *t*. These master equations however often prove a little unwieldy, and it is most of the time impossible to explicitly solve them to obtain analytical results. In practice, it is therefore common to sacrifice a little precision in the physical description of the problem at the microscopic scale, in order to obtain approximate equations that are easier to handle.

In the field of neutron fluctuations, an approximate modeling has been proposed in [3], which figures a recipe for turning the celebrated point-kinetic equations into a set of so-called Stochastic Differential Equations (SDE), that must be satisfied by the neutrons and precursors populations. From these SDE, numerical solutions for the mean and variance of the populations can easily be extracted, via a simple and straightforward numerical integration of the system of differential equations. This approximate treatment thus appears much easier to handle than the heavy machinery of the master equation technique. It has therefore seen growing interest in the neutronics literature. The approach has for instance been used to re-derive classical results about neutron detection statistics in [4]. More recently, SDE were used to study the stochastic behavior of the neutron population in a model figuring a phenomenological feedback mechanism [5]. Another interesting and fertile direction has been the use of the SDE method to model and study the behavior of dynamic Monte-Carlo simulations [6]. Finally, some more exotic use of SDE in the context of neutronics have also been recently published [7, 8, 9]. In these, SDE reveal an invaluable tool to propose an original, field theoretical, approach to neutron noise analysis. This original approach bridges the gap between neutronics and out-of-equilibrium statistical mechanics. It allowed for instance to interpret criticality as a second order phase transition [8].
The approximation scheme that is proposed in the present work, although providing equivalent results than the concurrent approach described in [3], follows a seemingly more tortuous, yet more robust, route. It consists in a direct and straightforward application of the textbook treatment of the matter, as can for instance be found in [10, 11, 12]. The alternative method that will be proposed has the advantage of highlighting more clearly the approximation that is being made in resorting to SDE to model the problem. This is so in particular because it forces one to explicitly implement a truncation in the master equation of the problem and because, after this is done, one is left with an approximate probability balance equation, the so-called Fokker-Planck equation, that can still be compared to the initial problem. In this setting, the SDE describing the approximate system is only obtained indirectly, as a byproduct of the Fokker-Planck equation: it can indeed be rigorously demonstrated that, associated to any Fokker-Planck equation, there corresponds a (family of) SDE [10].

To recap, the approximation method that is being proposed consists in the following successive steps: starting from a detailed microscopic model of the neutrons’ interactions, that is, from a given (forward) master equation, a unique Fokker-Planck equation is derived. The equation is obtained by explicitly performing a truncated expansion of the operators appearing in the master equation, hoping this does not alter too much the statistical properties of the solution. From this, a (family of) SDE can be readily extracted. These SDE can be quite easily numerically integrated and the solutions compared to the detailed microscopic solution of the problem.

Note that beyond providing a simple and convenient approximation scheme, there is a more fundamental reason for wanting to model neutron fluctuations using SDE. This is to make contact with the other known source of neutron noise, namely power neutron noise. There indeed exist another source of noise impacting the neutron population, it is the noise imprinted on the neutron population, originating from the external thermo-mechanical vibrations and constraints occurring in a nuclear reactor operating at high power.

One does not generally dispose, for the modeling of these complex random fluctuations, of a detailed microscopic model, and thus, their coupling to the neutron population must usually rely on a coarse grained approach (with notable exceptions, e.g. [13, 14]). One is naturally led to resort to a phenomenological, 'mesoscopic', approach, where the noise term resulting from the external vibrations of the reactor is modeled via SDE [15], which reveal an invaluable tool in situations where a microscopic model of the noise
source and of its coupling to the neutron population cannot be easily estab-
lished. By bridging the gap between the microscopic stochastic modeling of
zero power noise and the phenomenological description of power noise, SDE
seem to provide a unique and convenient framework in which both source of
noise could ultimately be incorporated and studied in a unified manner. Be-
cause it is generally (and rightly) assumed that power noise is much bigger
than zero power neutron noise, the application cases where such a unified
treatment would be necessary are scarce. It would in particular appear that
neutron fluctuations would only apply in situations where temperature feed-
back mechanisms, such as the Doppler or density effects, can be neglected.
There exists nevertheless various arguments in favor for a modeling including
feedback mechanisms, and the matter has already been extensively discussed
in the literature [5, 8, 9, 14, 16].

The present work revives a model first introduced in [17], which figures
a description of neutron fluctuations in a simplified, 0-dimensional, model in
which the neutron population is coupled to the temperature of the ambient
medium via a thermal feedback mechanisms. The stochastic microscopic
model is introduced in section 2. The forward master equation of the problem
is established, as well as the equations satisfied by the mean and variance
of the populations. Section 3 then deals with presenting the approximation
scheme that will be used to derive a system of SDE representative of the
overall problem. This approximation scheme is finally applied, in section 4,
to the stochastic model with thermal feedback.

## 2. A Stochastic Model with Temperature Feedback

A relatively simple reference model will be used throughout, to illustrate
the approximation method on a concrete application case. The latter is
straightforwardly adapted from a model found in the literature [17], with
only minor adjustments (namely, the addition of a source term and of a state
feedback mechanism on the heat extraction rate). Let one thus consider a
population of neutrons *N* (considered mono-energetic with constant velocity
*v*) and precursors *C* (only one group is considered, for simplicity), that evolve
in an infinite and homogeneous medium of temperature *T*. It will be consid-
ered that the population {*N*, *C*, *T*} evolve with time, under the influence of
the following random events:

*Neutron Capture* which occurs at a rate $\lambda_c = \Sigma_c v$ (per second and per neutron) and induces transitions in the population of the type $\{N, C, T\} \rightarrow$
$\{N-1, C, T + \tau_c\}$. For the sake of simplicity, it is assumed that each capture releases the exact same amount of energy in the ambient medium. Thus, each capture contribute to a fixed increase in temperature $\tau_c$. To account for the temperature feedback on the system's cross-sections, one can introduce a temperature dependency in the capture rate (e.g. [5]). A simple linear model with $\lambda_c = \lambda_c[T] = \lambda_{c,0} + \alpha T$ will be considered in the following.

*Precursor Decay* occurs at a rate $\lambda_d$ (per second and per precursor) and induces transitions in the population of the type $\{N, C, T\} \rightarrow \{N+1, C-1, T\}$.

*Neutron Source and/or Spontaneous Fission* which occurs at a rate $S$ (per second) and induces transitions in the population of the type $\{N, C, T\} \rightarrow \{N+i, C, T\}$ ($i \in \mathbb{N}$). The probability that $i$ neutrons are emitted by the source is denoted $\mu_i$. The mean and second factorial moment of this probability distribution will be denoted $\bar{\mu}$ and $\bar{\mu}(\bar{\mu}-1)$ respectively.

*Neutron Induced Fission* occurs at a rate $\lambda_f = \Sigma_f v$ (per second and per neutron) and induces transitions in the population of the type $\{N, C, T\} \rightarrow \{N-1+i, C+j, T+\tau_f\}$. The joint probability that $i$ neutrons and $j$ precursors are emitted by fission is denoted $\nu_{ij}$. It will be considered, for simplicity, that the marginal probabilities associated to neutrons and precursors are independent. One thus has $\nu_{ij} = \nu_{p,i} \times \nu_{d,j}$, where the index $p$ stands for the prompt neutrons and the index $d$ for the precursors (i.e. the delayed neutrons). The moments of $\nu_p$ and $\nu_d$ are written with the same nomenclature as for the neutron source.

*Temperature extraction* which occurs at a rate $\lambda_x[T] = \lambda_{x,0}T$ (per second) and induces transitions of the type $\{N, C, T\} \rightarrow \{N, C, T - \tau_x\}$.

Note that with the above convention, the reference temperature of the medium, that is, the temperature at which the feedback mechanisms vanish, is $T=0$. This reference temperature could be set to any desired value by a simple shift $T \rightarrow T - T_0$ in the functions $\lambda_c[T]$ and $\lambda_x[T]$.

In any sufficiently small interval $\Delta t$, the probability $P(N, C, T, t+\Delta t)$ to observe $N$ neutrons, $C$ precursors and a temperature $T$ at time $t+\Delta t$ can be deduced from the same data at time $t$ ([2, 1]). This observation is embodied
in the forward master equation,

$$
\begin{equation}
\begin{aligned}
\frac{d}{dt} P(N, C, T, t) ={}& \lambda_d \{(C+1)P(N-1, C+1, T, t) - CP(N, C, T, t)\} \\
& + \lambda_c [T - \tau_c](N+1)P(N+1, C, T - \tau_c, t) - \lambda_c [T] NP(N, C, T, t) \\
& + S \left\{ \sum_i \mu_i P(N-i, C, T, t) - P(N, C, T, t) \right\} \\
& + \lambda_f \left\{ \sum_{ij} \nu_{ij}(N+1-i)P(N+1-i, C-j, T - \tau_f, t) - NP(N, C, T, t) \right\} \\
& + \lambda_x [T + \tau_x] P(N, C, T + \tau_x, t) - \lambda_x [T] P(N, C, T, t).
\end{aligned}
\tag{1}
\end{equation}
$$

The equation appears quite formidable, and it would seem that there is
absolutely no hope of solving it directly. One can nevertheless make use
of this equation to derive a system of ordinary differential equations that
the mean populations must obey. Making use of the fact that, for instance,
$\langle N \rangle = \sum_{n,c,\theta} n P(n,c,\theta,t)$, one can readily show that,

$$
\frac{d}{dt} \langle N \rangle = \frac{\rho - \beta}{\Lambda} \langle N \rangle + \lambda_d \langle C \rangle + S \bar{\mu} - \alpha \langle NT \rangle, \quad (2)
$$

$$
\frac{d}{dt} \langle C \rangle = \frac{\beta}{\Lambda} \langle N \rangle - \lambda_d \langle C \rangle, \quad (3)
$$

$$
\frac{d}{dt} \langle T \rangle = \mathcal{Q} \langle N \rangle + \alpha \tau_c \langle NT \rangle - \lambda_{x,0} \tau_x \langle T \rangle. \quad (4)
$$

The following quantities have been introduced:

• $(\rho - \beta)/\Lambda = \lambda_f(\overline{\nu}_p - 1) - \lambda_{c,0}$. The parameter $\rho$ is the total reactivity of the system;

• $\beta = \overline{\nu}_d / (\overline{\nu}_p + \overline{\nu}_d)$ is the fraction of delayed neutrons;

• $\Lambda = 1/\lambda_f \times 1/(\overline{\nu}_p + \overline{\nu}_d)$ is the neutron generation time;

• $\mathcal{Q} = \lambda_f \tau_f + \lambda_{c,0} \tau_c$ is the total temperature increase rate induced by the neutrons' reactions.

One recognizes, in the system Eq. 2-4, the celebrated point kinetics equa-
tions, coupled to the ambient temperature of the medium via a thermal
feedback mechanism. One should note the crucial appearance of the moment $\langle NT \rangle$ in the equations. The system is not closed and must be supplemented with the equations for the second order moments. Those in turn involve moments of higher order, etc. There is a priori no reason to believe that this intrusion of higher moments stops at some point, that is, there is little hope that the system of equations closes at some order $n$ of the moments of the problem. To progress further, one has to impose some additional constraints on the problem.

## 2.1. The non-closure problem

One simplifying hypothesis would be to simply replace the covariance terms $\langle NT \rangle$ by the product $\langle N \rangle \times \langle T \rangle$ (note in particular that something of the sort must happen with SDE) but it is usually much better to apply a simplification at the level of the covariances of the problem. Introducing the standard notation $\mu_{IJ} = \langle IJ \rangle - \langle I \rangle \langle J \rangle - \langle I \rangle \delta_{IJ}$ [2, 1], one can show that the equations for the covariances of the populations read,

$$ \frac{d}{dt}\mu_{NN} = 2\frac{\rho - \beta}{\Lambda}\mu_{NN} + 2\lambda_d\mu_{NC} + \lambda_f\overline{\nu_p(\nu_p - 1)}\langle N \rangle + \overline{S\mu(\mu - 1)} \\ - \alpha[\langle N^2T \rangle - \langle N \rangle \langle NT \rangle - \langle NT \rangle], \quad (5) $$

$$ \frac{d}{dt}\mu_{CC} = -2\lambda_d\mu_{CC} + 2\frac{\beta}{\Lambda}\mu_{NC} + \overline{\nu_d(\nu_d - 1)}\langle N \rangle, \quad (6) $$

$$ \frac{d}{dt}\mu_{NC} = \left(\frac{\rho - \beta}{\Lambda} - \lambda_d\right)\mu_{NC} + \frac{\beta}{\Lambda}\mu_{NN} + \lambda_d\mu_{CC} + \lambda_f\overline{\nu_p\nu_d}\langle N \rangle \\ - \alpha[\langle NCT \rangle - \langle NT \rangle \langle C \rangle], \quad (7) $$

$$ \frac{d}{dt}\mu_{CT} = Q\mu_{NC} + \frac{\beta}{\Lambda}\mu_{NT} - \lambda_d\mu_{CT} - \lambda_{x,0}\tau_x\mu_{CT} + \frac{\beta}{\Lambda}\tau_f\langle N \rangle \\ + \alpha\tau_c[\langle NCT \rangle - \langle NT \rangle \langle C \rangle], \quad (8) $$

$$ \frac{d}{dt}\mu_{NT} = Q\mu_{NN} + \frac{\rho - \beta}{\Lambda}\mu_{NT} + \lambda_d\mu_{CT} - \lambda_{x,0}\tau_x\mu_{NT} + \tau_f\left\{\frac{\rho - \beta}{\Lambda} + \lambda_{c,0} + \lambda_f\right\}\langle N \rangle \\ + \alpha\tau_c[\langle N^2T \rangle - \langle NT \rangle (\langle N \rangle + 1)] - \alpha[\langle NT^2 \rangle - \langle NT \rangle \langle T \rangle], \quad (9) $$
$$
\begin{align*}
\frac{d}{dt}\mu_{TT} &= 2Q\mu_{NT} + [\lambda_{c,0}\tau_c(\tau_c - 1) + \lambda_f\tau_f(\tau_f - 1)]\langle N \rangle \\
&\quad - 2\lambda_{x,0}\mu_{TT} + \lambda_{x,0}\tau_x(\tau_x - 1)\langle T \rangle \\
&\quad - \alpha\tau_c[2\langle NT^2 \rangle - \langle NT \rangle \langle T \rangle - \langle NT \rangle] + \alpha\tau_c^2(\mu_{NT} + \langle N \rangle \langle T \rangle) \quad (10)
\end{align*}
$$

As anticipated, the system of equations depends on moments of higher order ($\langle NCT \rangle$, etc.). To progress any further, constraints must be imposed on the moments of order $\ge 3$. Probably the most elegant approach allowing to easily tackle the problem consists in postulating that the probability function $P(N, C, T, t)$ is normal distributed. In this case, all of the moments of the problem can be written as a sum of products of mean and covariances [2]. One can then use standard results in probability theory to treat the terms of high order and close the system of equations. The simplified system of equations obtained with this assumption is provided in Appendix A.

This 'Gaussian model' will be used, in the following, as the reference theoretical prediction, to which all of the simulation results and other approximation schemes will be compared.

### 3. The Fokker-Planck equation and Stochastic Differential Equations

It was already stated in the introduction that the whole subject of neutron fluctuations can somehow be viewed as the matter of characterizing the imprint of the microscopic random nature of the neutrons' interactions on the macroscopic averaged observables that one is interested in. In this interpretation, SDE are but a description of the problem at an intermediate scale, lying in between the microscopic and macroscopic scales. It indeed concretely consists in a separation between the microscopic and random degrees of freedom of the system, along with the provision of a continuous function that describes the evolution of the quantity of interest at a scale that is really neither truly macroscopic, nor microscopic. In the present case, this is reflected in the nature of the populations {$N, C, T$}, which vary at the microscopic level in discrete increments. It is evident that one cannot directly use these quantities to establish a set of differential equations, which require continuous functions, at a minimum. Therefore, the first step in deriving a set of SDE is to introduce coarse-grained continuous functions for the system. Various methods from the literature can be adapted for this purpose [11, 18, 19].
These methods are crucial in controlling the accuracy of the SDE treatment, but their detailed examination in the context of neutron fluctuations is left for future work.

On the model of the earlier work on the subject, [3, 4, 5, 6], it will simply be considered that the overall populations sizes are, at any time, much larger than the typical random increments of the underlying stochastic process. In other words, it will be considered that the increments of the stochastic process are real and that the populations {$N(t), C(t), T(t)$} are continuous functions. In this most simple scheme, the Fokker-Planck approximation that is introduced in the present section might appear somehow unjustified and unnatural. One should thus keep in mind that more elaborate schemes could be employed to better justify the approximation.

### 3.1. The Fokker-Planck equation

The master equation Eq. 1 can be written more compactly by formally introducing translation operators $\hat{E}_m^I (m \in \mathbb{R} \text{ and } I \in \{N, C, T\})$, such that, for example,

$$ \hat{E}_m^N f(N, C, T, t) = f(N + m, C, T, t), \quad (11) $$

for any function $f(N, C, T, t)$. Equipped with these operators, one can show that the master equation Eq. 1 takes the form,

$$ \begin{aligned} \frac{d}{dt} P(N, C, T, t) = & \left( \hat{E}_1^N \hat{E}_{-\tau_c}^T - 1 \right) \lambda_c[T] NP(N, C, T, t) \\ & + \lambda_d \left( \hat{E}_1^C \hat{E}_{-1}^N - 1 \right) CP(N, C, T, t) \\ & + S \left( \sum_i \mu_i \hat{E}_{-i}^N - 1 \right) P(N, C, T, t) \\ & + \left( \hat{E}_{\tau_x}^T - 1 \right) \lambda_x[T] P(N, C, T, t) \\ & + \lambda_f \left( \sum_{ij} \nu_{ij} \hat{E}_{1-i}^N \hat{E}_{-j}^C \hat{E}_{-\tau_f}^T - 1 \right) NP(N, C, T, t). \end{aligned} \quad (12) $$

At the coarse grained level at which the system is observed, the translation operators can be given an explicit form [12, 20],

$$ \hat{E}_m^I = \exp(m\partial_I) = 1 + m \frac{\partial}{\partial I} + \frac{m^2}{2} \frac{\partial^2}{\partial I^2} + \dots \quad (13) $$
It is tempting to actually stop the series expansion of the exponential at, say, the second order, and see what happens. Assume concretely, that one can exactly write,

$$ \hat{E}_m^I = 1 + m \frac{\partial}{\partial I} + \frac{m^2}{2} \frac{\partial^2}{\partial I^2}. \qquad (14) $$

Injecting this expression into the master equation Eq. 12, one obtains the celebrated Fokker-Planck equation,

$$ \frac{d}{dt} P(N, C, T, t) = - \sum_{I=N,C,T} \partial_I [A_I P(N, C, T, t)] + \frac{1}{2} \sum_{IJ} \{[BB^t]_{IJ} P(N, C, T, t)\}, \qquad (15) $$

with, in the present case,

$$ \vec{A} = \begin{pmatrix} \frac{\rho-\beta}{\Lambda} N + \lambda_d C + S \bar{\mu} - \alpha NT \\ \frac{\beta}{\Lambda} N - \lambda_d C \\ QN - \lambda_{x,0} \tau_x T + \alpha \tau_c NT \end{pmatrix} \qquad (16) $$

and,

$$ \mathbf{BB}^t = \begin{pmatrix} [BB^t]_{NN} & [BB^t]_{NC} & [BB^t]_{NT} \\ [BB^t]_{CN} & [BB^t]_{CC} & [BB^t]_{CT} \\ [BB^t]_{TN} & [BB^t]_{TC} & [BB^t]_{TT} \end{pmatrix}, \qquad (17) $$

with,

* $[BB^t]_{NN} = \left[ \frac{1}{\Lambda} \frac{\nu_p(\nu_p-1)}{\nu_p+\bar{\nu}_d} - \frac{\rho-\beta}{\Lambda} \right] N + S \left[ \bar{\mu}(\mu-1) + \bar{\mu} \right] + \lambda_d C + \alpha NT$;

* $[BB^t]_{CC} = \left[ \frac{1}{\Lambda} \frac{\nu_d(\nu_d-1)}{\nu_p+\bar{\nu}_d} + \frac{\beta}{\Lambda} \right] N + \lambda_d C$;

* $[BB^t]_{TT} = (\lambda_f \tau_f^2 + \lambda_{c,0} \tau_c^2 + \alpha \tau_c^2 T) N + \lambda_{x,0} \tau_x^2 T$;

* $[BB^t]_{NC} = [BB^t]_{CN} = \frac{\beta}{\Lambda}(\bar{\nu}_p - 1)N - \lambda_d C$;

* $[BB^t]_{NT} = [BB^t]_{TN} = -QN + \tau_f(\frac{\rho-\beta}{\Lambda} + \lambda_f + \lambda_{c,0})N - \alpha \tau_c NT$;

* $[BB^t]_{CT} = [BB^t]_{TC} = \frac{1}{2}\frac{\beta}{\Lambda}\tau_f N$.
### 3.2. Stochastic Differential Equations

The derivation of a system of SDE for the problem is now but a textbook result [10]. It can indeed be shown that to the above Fokker-Planck equation Eq. 15, there correspond a system of SDE that takes the form (with $\vec{X} = (N, C, T)$),

$$d\vec{X} = \vec{A}(\vec{X}, t)dt + \mathbf{B} d\vec{W}, \quad (18)$$

where the components of the vector $d\vec{W}$ are independent Wiener processes.

There thus would appear to be little extra work to establish an approximate system of SDE for the problem. One only needs to find the "square root" of the symmetric matrix $\mathbf{BB}^t$. Note that the matrix $\mathbf{B}$ appearing in the SDE is only defined up to multiplication by any orthogonal matrix $\mathbf{O}$ (because $(\mathbf{BO}) \times (\mathbf{BO})^t = BB^t$). The matrix $\mathbf{B}$ itself, and therefore, the whole system of SDE representing the system is not unique.

## 4. Application to the stochastic model

In this last section, the approximation scheme associated to Eq. 15 is applied to the microscopic point kinetics model with temperature feedback that was introduced in section 2. The approach will be incremental, and the different populations are all introduced individually, one after the other. This staged approach will reveal useful to introduce, in particular, the modeling strategies and the parameters selections that will be used throughout.

### 4.1. Neutrons

As a first application case, consider a simplified version of the stochastic model, where one imposes $C = T = 0$ at any time in the master equation Eq. 12. It can then easily be shown that the Fokker-Planck approximation introduced in the previous section yields the SDE,

$$dN(t) = \left[ \frac{\rho}{\Lambda} N + S \bar{\mu} \right] dt + B(N) dW(t), \quad (19)$$

with $B(N) = \sqrt{\left[ \frac{1}{\Lambda} \frac{\overline{\nu_p(\nu_p-1)}}{\overline{\nu_p}} - \frac{\rho}{\Lambda} \right]} N + S \left[ \overline{\mu} + \overline{\mu(\mu-1)} \right]$. This result is in agreement with the latest literature [4, 5, 8], and also agrees with the concurrent approach documented in [3], provided one make the further approximation $\overline{\nu_p^2} = \overline{\nu_p^2}$. This slight difference will be further discussed in the next section.
The Figure 1 presents the result of the numerical integration of the SDE, along with a comparison with the mean and variance derived from the master equation simulation. The simulations were performed using the following set of parameters (all of the parameters used in the simulations can be directly deduced from those listed below):

* $\overline{\nu}_p = 2.5$ and $\overline{\mu} = 5$. For simplicity, it was assumed that the neutron emission probability distributions both follow Poissonian profiles. In this case, it can be shown that $\overline{\nu}_p(\overline{\nu}_p - 1) = \overline{\nu}_p^2$ and $\overline{\mu}(\overline{\mu} - 1) = \overline{\mu}^2$.

* $\Lambda = 0.1$ s. Note that this value ultimately sets the dynamical scale of the problem (along with the source strength $S$). As such, it can always be adjusted to whatever value by a global rescaling of the time coordinate. For the numerical integration of the SDE, increments of size $dt = 10^{-3}$ s were used.

* $\rho = -500$ pcm and $S = 10$ s$^{-1}$. The source term is chosen so that it exactly compensates the negative reactivity. The neutron population is thus stationary and one imposes $N(t = 0) = N_{\text{stat.}} = S\overline{\mu}/|\rho| = 1000$.

The SDE Eq. 19 is integrated 100 times using the Euler-Maruyama method [21]. This process generates the hundred individual trajectories, displayed as light blue plain curves on the upper plot of the Figure 1. The mean behavior of the population is obtained by performing an ensemble average over these curves at each time. The mean and variance resulting from this operation are shown as plain, deep blue curves on both the upper (mean) and lower (variance) plots of the figure.

The results of the approximate SDE treatment can be compared to a direct simulation of the microscopic stochastic process, performed with the well known and well documented Gillespie algorithm [22]. The approach is essentially the same as above: 100 individual trajectories are simulated and the mean and variance are estimated by taking the ensemble average over this population. The mean and variance obtained from the simulation of the microscopic process are shown as plain red curves on the Figure 1.

Finally, both results are compared to the theoretical predictions for the mean and variance of the neutron populations, that follows from solving the moment equations established section 2. The equations for the mean and
<figure><img src="image_1.png" /><figcaption>Figure 1: Comparison of the numerical integration of the SDE Eq. 19 with the simulation of the associated master equation. In both cases, the mean (upper plot) and variance (lower plot) are compared to the solution of the theoretical prediction, extracted from the system Eq. 20-21.</figcaption></figure>

rescaled variance $\mu_{NN}$ read, in this simplified setting,

$$ \frac{d}{dt}\langle N \rangle = \frac{\rho}{\Lambda}\langle N \rangle + S\bar{\mu}, \qquad (20) $$

$$ \frac{d}{dt}\mu_{NN} = 2\frac{\rho}{\Lambda}\mu_{NN} + \lambda_f\overline{\nu_p(\nu_p - 1)}\langle N \rangle + S\overline{\mu(\mu - 1)}. \qquad (21) $$

Both equations can be shown to possess analytical solutions [1, 2]. The latter are shown, for reference, as dotted green curves on both plots in Figure 1.

A very good agreement is obtained for both the mean and the variance of the neutron population. It would thus appear that the SDE treatment of the problem, and the underlying Fokker-Planck approximation, manage to reproduce quite well the main features of the detailed model (with more than 10 times faster computation times than the Gillespie algorithm).

A more detailed inspection of the problem would however still be needed to completely characterize the fidelity with which one is able to reproduce,
with the Fokker-Planck approximation, the underlying complexity of the microscopic stochastic process. In this particularly simple regime, where analytical solutions (at least partially and with additional simplifications on the source term [23]) of both the master equation and the Fokker-Planck equation could be obtained, it would be interesting to completely characterize the agreement of both probability distributions. This careful inspection of the problem is left for future work.

## 4.2. Inclusion of the precursors

The simplest next step is the inclusion of the precursor population into the problem. Setting $T = 0$ in Eq. 12, one ends up with the system of coupled SDE,

$$ \left( \frac{dN}{dC} \right) = \left( \frac{\rho - \beta}{\Lambda} N + \lambda_d C + S \bar{\mu} \right) dt + \mathbf{B} \left( \frac{dW_1}{dW_2} \right), \quad (22) $$

where the matrix $\mathbf{B}$ is whatever matrix satisfying,

$$ \mathbf{BB}^t = \begin{pmatrix} \left[ \frac{1}{\Lambda} \frac{\overline{\nu_p(\nu_p-1)}}{\overline{\nu_p} + \overline{\nu_d}} - \frac{\rho - \beta}{\Lambda} \right] N + S \left[ \overline{\mu} + \overline{\mu}(\mu - 1) \right] + \lambda_d C & \frac{\beta}{\Lambda} (\overline{\nu_p} - 1) N - \lambda_d C \\ \frac{\beta}{\Lambda} (\overline{\nu_p} - 1) N - \lambda_d C & \left[ \frac{\beta}{\Lambda} + \frac{1}{\Lambda} \frac{\overline{\nu_d}(\overline{\nu_d} - 1)}{\overline{\nu_p} + \overline{\nu_d}} \right] N + \lambda_d C \end{pmatrix}. \quad (23) $$

The SDE system that is obtained agrees with the results of the concurrent method [3] only if one further assumes that $\overline{\nu}_i^2 \approx \overline{\nu}_i^2$ ($i = p, d$). Because the present approach is more firmly based on the underlying microscopic stochastic process, this further hypothesis isn't necessary and one is thus able to better reproduce the microscopic model of the problem. Apart from this difference, both models can be shown to be in total correspondence, and one would thus expect to find similar results with both approaches. A further investigation of this model thus doesn't appear really necessary, and the interested reader is referred to the original literature [3] for more information and concrete application of the approximate equations.

## 4.3. The full model with temperature feedback

The full model introduced in section 2 is finally considered. To this end, much more data must be provided in order to completely characterize the problem. The data used earlier are thus supplemented with the following set of parameters:
* The fraction of delayed neutrons $\beta$ is set to one percent. The typical decay constant for the precursors is set so that $\lambda_d = \ln 2 \text{ s}^{-1}$. The initial number of precursors is set to zero, i.e. $C(0) = 0$.

* The initial temperature is set to $T(0) = 0$.

* The temperature increase rate is set to $Q = 0.01 \text{ s}^{-1}$. On the basis that fission releases around 200 MeV of energy, while radiative capture releases around 10 MeV, one sets $\tau_f = 20 \tau_c$.

* The increase rate of the capture cross-section with temperature is set to $\alpha = 0.005$ (per second and per degree).

* The temperature extraction mechanism is fully specified by both the extraction rate $\lambda_{x,0} = 4.6 \text{ s}^{-1}$ and the associated temperature $\tau_x = 50 \tau_f$. These parameters were manually adjusted, so that the system reaches stationarity (or at least appears stationary) in a sufficiently short time interval.

Following the Fokker-Planck approximation of section 3, the problem can be described as the system of SDE,

$$ \left( \begin{matrix} dN \\ dC \\ dT \end{matrix} \right) = \left( \begin{matrix} \frac{\rho - \beta}{\Lambda} N + \lambda_d C + S \bar{\mu} - \alpha N T \\ \frac{\beta}{\Lambda} N - \lambda_d C \\ Q N - \lambda_{x,0} \tau_x T + \alpha \tau_c N T \end{matrix} \right) dt + \mathbf{B} \left( \begin{matrix} dW_1 \\ dW_2 \\ dW_3 \end{matrix} \right). \quad (24) $$

The last ingredient that must be specified is the matrix $\mathbf{B}$, whose value must be deducted from the product $\mathbf{B}\mathbf{B}^t$ of Eq. 17. The most direct approach along this way consists in diagonalizing the matrix product, and deduce $\mathbf{B}$ from the square root of the diagonal matrix. Because of the multiple changes of basis, the resulting matrix looks quite involved, but it can be quite readily used to numerically integrate the SDE.

As in the previous case, the system of SDE is numerically integrated 100 times, using the Euler-Maruyama scheme. The results are compared to the theoretical predictions extracted from the Gaussian model established in section 2, but also to a direct simulation of the stochastic process, using, here also, the Gillespie algorithm.

The Figure 2 shows the temporal evolution of the neutron population as a result of the numerical integration of the SDE system and the direct simulation of the underlying stochastic process. Both results agree with
<figure><img src="image_2.png" /><figcaption>Figure 2: Comparison of the neutron population extracted from the numerical integration of the SDE Eq. 24 with the simulation of the microscopic stochastic process. The mean (upper plot) and variance (lower plot) of the neutron population are compared to the Gaussian model.</figcaption></figure>

the reference Gaussian model, whose solution is derived from a numerical integration of the associated system of moment equations shown in Appendix A.

Note the similarity between the behavior observed in Figure 2 and those obtained previously for the case where only neutrons were considered. In both cases the neutron population appears stationary (except for a very small and short transient at the beginning of the simulations, and invisible on the figure, in the current model). As for the variance of the neutron population, one observe in both cases a transient increase, followed by a stabilization of the variance at a relatively high value. One can observe that the introduction of a feedback mechanism contributes to considerably shorten the transient evolution of the variance. The latter stabilizes sooner, at a lower value (almost half the original variance observed for neutrons alone). The feedback thus act as a stabilizing agent, decreasing the overall fluctuations in the neutron population.

Very similar results can be observed for the temperature of the ambient
medium, which is displayed on the Figure 3. After an initial increase of the temperature, during which the neutron population and the temperature smoothly settle to their equilibrium positions, the temperature stabilizes to the equilibrium value predicted by the theoretical model.

<figure><img src="image_3.png" /><figcaption>Figure 3: Comparison of the temperature (expressed in arbitrary units) extracted from the numerical integration of the SDE Eq. 24 with the simulation of the microscopic stochastic process. The mean (upper plot) and variance (lower plot) are compared to the Gaussian model.</figcaption></figure>

For the sake of completeness, the results for the precursor population are shown in Figure 4. Here also, one observe that the precursor population stabilizes relatively quickly, following the general equilibrium resulting from the competition between neutron-induced reactions and temperature.

From the good agreement observed on the above figures, it appears that the approximate SDE treatment reveals very accurate and is able to reproduce precisely the overall features of this quite complex microscopic model. As a final note, one could wonder if the agreement extends beyond the means and variances of the populations, and ask, for instance, if the overall probability distributions $P(N, C, T, t)$, for any fixed time $t$, agree. A proper and complete investigation of the matter would be quite involved, and would far
<figure><img src="image_4.png" /><figcaption>Figure 4: Comparison of the precursor population extracted from the numerical integration of the SDE Eq. 19 with the simulation of the associated master equation. The mean (upper plot) and variance (lower plot) are compared to the Gaussian closure model, whose equations are provided in Appendix A.</figcaption></figure>

exceed the scope of the present introduction. To get a sense of how well the SDE is capable of reproducing the overall probability distribution, the Figure 5 shows the distribution of the number of neutrons, precursors and the temperature of the medium at $t = 40$ s for each individual trajectories simulated, either from the SDE integration (in orange) or from the simulation of the stochastic process (in blue). In all three cases, the distributions obtained from the two simulation methods appear quite similar, at least by visual inspection. This fact was however already reflected in the good agreement that was found for the means and variances of the populations. It would be difficult, in the present case to really go any further, because this would necessitate to explore rarer events, corresponding to the tails of the probability distributions. This in turn would necessitate more statistics to properly characterize the overall agreement of the distributions.

It would for instance be interesting to evaluate the predictions of the SDE model associated to the extinction or survival probabilities. This would
<figure><img src="image_5.png" /><figcaption>Figure 5: Distribution of the number of neutrons, number of precursors and of the temperature of the medium at time $t = 40s$ extracted from the individual independent trajectories simulated with the SDE model (orange) or the Gillespie algorithm (blue).</figcaption></figure>

however probably require the introduction of the more refined and more powerful approximation schemes cited at the beginning of the section. This further inspection of the properties of the SDE is thus left for future work.

## 5. Conclusion

Stochastic Differential Equations are a versatile and powerful tool for simplifying the modeling of systems subjected to a source of noise. They allow for the modeling of a problem on a mesoscopic scale, where the microscopic degrees of freedom are clearly isolated and treated separately. As was clearly established in the present work, this approximate treatment corresponds to a truncation of the stochastic microscopic model. In this setting, the underlying probability function associated with the solution of the SDE is described by a Fokker-Planck equation, which is supposed to correctly capture the main characteristics of the problem.

Because of its versatility and ease of manipulation, the SDE formalism is gaining popularity for modeling neutron noise. However, its relationship
to the more traditional treatment, i.e., the master equation approach, has never been clearly established, and the underlying hypotheses and required approximations have never been clearly stated, at least in the neutronics literature. The goal of the present work was therefore to fill this gap and clearly expose the thread of approximations that connects the two approaches to the problem. One of the main advantages of the Fokker-Planck approach proposed here is that it clearly identifies the weak points in the approximation scheme, thus revealing potential areas for improvement. Among the possible improvements of the approximation scheme, the implementation of the refined system size expansion method [11], would for instance appear as promising candidate to better justify the truncation of the master equation. Its application in the context of neutron fluctuations, would thus appear a natural extension of this work.

The SDE formalism has been applied in the relatively original context of a so-called zero power noise model including feedback mechanisms. However weak they may be, the thermal feedback mechanisms can be shown to play a very important role to correctly describe the statistical properties of the neutron population. For example, these mechanisms have been shown to play an important role in the appearance of possible *noise-induced transitions* [24], which would also seem to be a natural lead for further exploration.
Appendix A. Moment equations after Gaussian closure

Assuming that the probability function P(N, C, T, t) is Gaussian, it can be readily shown that the system of equations established in section 2.1, reads,

$$ \frac{d}{dt}\langle N \rangle = \frac{\rho - \beta}{\Lambda}\langle N \rangle + \lambda_d\langle C \rangle + S\bar{\mu} - \alpha(\mu_{NT} + \langle N \rangle\langle T \rangle), $$

$$ \frac{d}{dt}\langle C \rangle = \frac{\beta}{\Lambda}\langle N \rangle - \lambda_d\langle C \rangle, $$

$$ \frac{d}{dt}\langle T \rangle = Q\langle N \rangle - \lambda_{x,0}\tau_x\langle T \rangle + \alpha\tau_c(\mu_{NT} + \langle N \rangle\langle T \rangle), $$

$$ \frac{d}{dt}\mu_{NN} = 2\frac{\rho - \beta}{\Lambda}\mu_{NN} + 2\lambda_d\mu_{NC} + \lambda_f\overline{\nu_p(\nu_p - 1)}\langle N \rangle + S\bar{\mu}(\mu - 1) \\ - 2\alpha(\mu_{NT} + \langle N \rangle\langle T \rangle)(\langle N \rangle + 1) $$

$$ \frac{d}{dt}\mu_{CC} = -2\lambda_d\mu_{CC} + 2\frac{\beta}{\Lambda}\mu_{NC} + \overline{\nu_d(\nu_d - 1)}\langle N \rangle $$

$$ \frac{d}{dt}\mu_{NC} = \left(\frac{\rho - \beta}{\Lambda} - \lambda_d\right)\mu_{NC} + \frac{\beta}{\Lambda}\mu_{NN} + \lambda_d\mu_{CC} + \lambda_f\overline{\nu_p\nu_d}\langle N \rangle \\ + \alpha(\mu_{NT} + \langle N \rangle\langle T \rangle)\langle C \rangle $$

$$ \frac{d}{dt}\mu_{CT} = Q\mu_{NC} + \frac{\beta}{\Lambda}\mu_{NT} - \lambda_d\mu_{CT} - \lambda_{x,0}\tau_x\mu_{CT} + \frac{\beta}{\Lambda}\tau_f\langle N \rangle \\ - \alpha\tau_c(\mu_{NT} + \langle N \rangle\langle T \rangle)\langle C \rangle $$

$$ \frac{d}{dt}\mu_{NT} = Q\mu_{NN} + \frac{\rho - \beta}{\Lambda}\mu_{NT} + \lambda_d\mu_{CT} - \lambda_{x,0}\tau_x\mu_{NT} + \tau_f\left\{\frac{\rho - \beta}{\Lambda} + \lambda_{c,0} + \lambda_f\right\}\langle N \rangle \\ - \alpha\tau_c(\mu_{NT} + \langle N \rangle\langle T \rangle)(\langle N \rangle + 1) + \alpha\tau_c(\mu_{NT} + \langle N \rangle\langle T \rangle) $$

$$ \frac{d}{dt}\mu_{TT} = 2Q\mu_{NT} + [\lambda_{c,0}\tau_c(\tau_c - 1) + \lambda_f\tau_f(\tau_f - 1)]\langle N \rangle - 2\lambda_{x,0}\mu_{TT} + \lambda_{x,0}\tau_x(\tau_x - 1)\langle T \rangle \\ + \alpha\tau_c[2\mu_{TT}\langle N \rangle + \mu_{NT}(2\langle T \rangle - 1)] + \alpha\tau_c^2(\mu_{NT} + \langle N \rangle\langle T \rangle) $$
References

[1] I. Pázsit, L. Pál, Neutron fluctuations: a treatise on the physics on branching processes, 1st ed ed., Elsevier, Amsterdam ; Boston, 2008.

[2] M.M.R. Williams, Random Processes in Nuclear Reactors, pergamon ed., 1974.

[3] J. G. Hayes, E. J. Allen, Stochastic point-kinetics equations in nu-
clear reactor dynamics, Annals of Nuclear Energy 32 (2005) 572–587.
doi:10.1016/j.anucene.2004.11.009.

[4] C. Dubi, R. Atar, Modeling neutron count distribution in a subcritical core by stochastic differential equations, Annals of Nuclear Energy 111 (2018) 608–615. doi:10.1016/j.anucene.2017.09.040.

[5] C. Dubi, A. K. Prinja, A Diffusion Scale Approximation for Stochastic Fission Chains with State Feedback, Nuclear Science and Engineering (2024) 1–15. doi:10.1080/00295639.2024.2324532.

[6] T. Z. Hajas, G. Tolnai, M. Margoczi, D. Legrady, Noise term modeling of dynamic monte carlo using stochastic differential equations 195 (2024) 110061. doi:10.1016/j.anucene.2023.110061.

[7] R. Garcia-Millan, J. Pausch, B. Walter, G. Pruessner, Field-theoretic approach to the universality of branching processes, Physical Review E 98 (2018) 062107. doi:10.1103/PhysRevE.98.062107, arXiv:1808.08418 [cond-mat, physics:math-ph].

[8] B. Dechenaux, T. Delcambre, E. Dumonteil, Percolation properties of the neutron population in nuclear reactors, Physical Review E 106 (2022) 064126. doi:10.1103/PhysRevE.106.064126.

[9] B. Dechenaux, Field-theoretic approach to neutron noise in nuclear reactors, Physical Review E 109 (2024) 044145. doi:10.1103/PhysRevE.109.044145.

[10] C. W. Gardiner, Stochastic methods: a handbook for the natural and social sciences, number 13 in Springer series in synergetics, 4th ed ed., Springer, Berlin Heidelberg, 2009.
[11] N. Van Kampen, Stochastic Processes in Physics and Chemistry, 3rd edition ed., North Holland, 1981.

[12] A. Altland, B. D. Simons, Condensed Matter Field Theory, 2 ed., Cambridge University Press, 2010. doi:10.1017/CBO9780511789984.

[13] L. Pál, I. Pázsit, Theory of Neutron Noise in a Temporally Fluctuating Multiplying Medium, Nuclear Science and Engineering 155 (2007) 425–440. doi:10.13182/NSE07-A2674.

[14] I. Pázsit, Z. Kuang, A. Prinja, A unified theory of zero power and power reactor noise via backward master equations, Annals of Nuclear Energy 29 (2002) 169–192. doi:10.1016/S0306-4549(01)00029-9.

[15] Handbook of nuclear engineering, 2010. doi:10.1007/978-0-387-98149-9.

[16] C. Sunde, et al., Research and Development Program in Reactor Diagnostics and Monitoring with Neutron Noise Methods: Stage 11 and 12. Final Report, Technical Report, Report 2006:34, SKI (Swedish Nuclear Power Inspectorate), 2006.

[17] D. R. Harris, V. Prescop, Stability and stationarity of a reactor as a stochastic process with feedback, Nuclear Science and Engineering 37 (1969) 171–179. doi:10.13182/NSE69-A20675.

[18] D. T. Gillespie, The chemical Langevin equation, The Journal of Chemical Physics 113 (2000) 297–306. doi:10.1063/1.481811.

[19] F. Benitez, C. Duclut, H. Chaté, B. Delamotte, I. Dornic, M. A. Muñoz, Langevin equations for reaction-diffusion processes, Physical Review Letters 117 (2016) 100601. doi:10.1103/PhysRevLett.117.100601, arXiv:1604.06740 [cond-mat].

[20] J. J. Sakurai, Advanced quantum mechanics, Addison-Wesley, Reading, Mass. Bonn, 1967.

[21] M. S. Miguel, R. Toral, Stochastic Effects in Physical Systems (1997). doi:10.48550/arXiv.cond-mat/9707147, arXiv:cond-mat/9707147 version: 1.
[22] D. T. Gillespie, Stochastic Simulation of Chemical Kinetics, Annual Review of Physical Chemistry 58 (2007) 35–55.
doi:10.1146/annurev.physchem.58.032806.104637.

[23] G. I. Bell, Probability distribution of neutrons and precursors in a multiplying assembly, Annals of Physics 21 (1963) 243–283.
doi:https://doi.org/10.1016/0003-4916(63)90108-8.

[24] W. Horsthemke, R. Lefever, Noise-Induced Transitions, volume 15 of Springer Series in Synergetics, Springer Berlin Heidelberg, 2006.
doi:10.1007/3-540-36852-3.