© A. Fauvel et al., Published by EDP Sciences, 2025

1. Introduction

Power reactor neutron noise refers to the small variations of the neutron flux around its steady state [1]. Noise is typically induced by various kinds of perturbations resulting from fluid-structure interactions, encompassing the vibration of fuel pins and control rods, fuel assemblies or the full vessel.

The perturbations are conveyed by the neutron flux throughout the core, so that the pertaining information can in principle be retrieved even by ex-core neutron detectors located far away from the noise source. In this respect, noise monitoring provides a useful diagnostic tool to infer the occurrence, nature and location of core perturbations by measuring the deviations and correlations of multiple detectors collecting noise spectra [1]. The potential of neutron noise analysis to detect and localize core anomalies in their early stage is of interest for industrial applications, and thus calls for the development of accurate, efficient and robust simulation tools capable of solving the equations describing the propagation of the noise field. This has fostered the creation or the improvement of a wide class of deterministic and Monte Carlo neutron noise solvers either in the time or frequency domain[2–10], especially under the auspices of the CORTEX H2020 European project (2017-2021) [11].

In this work, we will focus on Monte Carlo simulation for noise problems, bearing in mind that the use of Monte Carlo methods allows for reference calculations with almost no approximations, as opposed to faster but approximate deterministic methods [12, 13]. Almost all Monte Carlo methods for noise analysis aim at solving the noise equations in the (Fourier-transformed) frequency domain, which takes the form of a fixed-source transport problem with complex-valued source and operators, and whose solution is also complex-valued. Following an idea introduced by Yamamoto [14], the complex terms occurring in the noise equation are handled by allowing particles to carry complex-valued weights. Two main methods have been proposed, differing on the interpretation of the Fourier-transformed time-derivative occurring in the noise equation. Yamamoto’s method is flight-based: the time-derivative is taken into account during flights, rotating the particle weight in the complex plane [14]. Rouchon et al.’s method is collision-based: the time-derivative is interpreted as a special copy collision event, with real cross-section and complex yield [15].

Both methods, while effective within the so-called plateau region of the frequency domain, have shown poor convergence properties outside the plateau [16]. A peculiar behavior that may occur concerns the symmetric divergence of the particle populations having opposite signs: the average contribution to the sought response is collectively compensated, but the variance may grow dramatically, while the simulation is prevented from ending [16]. To forestall these effects, population control methods have been introduced: weight cancellation, which consists in annihilating particles of opposite sign, has a positive impact on the figure of merit [14, 17, 18] and on convergence [17]. Besides, other sampling strategies, encompassing implicit capture, forced fission, Russian roulette [16] or particle splitting into real and imaginary parts [19], also non-trivially affect the convergence.

In a previous work, we proposed a preliminary formulation of a general theoretical framework for the analysis of the convergence of Monte Carlo algorithms for neutron noise [19]. We have shown in particular that convergence can be assessed in terms of a family of eigenvalues. The main limitation of this investigation was that it focused on infinite homogeneous media, without spatial effects. In this paper, we set out to considerably extend the perimeter of validity of our convergence analysis, in view of establishing a comprehensive theoretical framework supported by numerical evidence in more realistic systems.

This work is organized as follows. In Section 2 we briefly recall the neutron noise equation. The main features of the Monte Carlo algorithms for noise problems are described in Section 3. Then, in Section 4, we propose a formalism for Monte Carlo games associated to a class of complex-valued integral equations. Relying on these findings, we show that the convergence criterion of neutron noise algorithms can be mapped onto two families of eigenvalues and we discuss the effects of weight cancellation. In order to substantiate these results, in Section 5 we investigate a few simple benchmark configurations where analytical reference solutions can be derived. Comparison to previous literature findings and recommendations on implementation are further examined in Section 6. Conclusions are finally drawn in Section 7.

2. The noise equations

In order for this paper to be self-contained, we begin by briefly recalling the derivation of the neutron noise equation. For more details, the readers are referred to e.g. Ref. [1]. For conciseness, the position $\overrightarrow{r}$, direction $\overrightarrow{\mathrm{\Omega }}$ and energy E variables are implicit, whenever possible. The time-dependent Boltzmann equation is written symbolically

$$\begin{array}{c}\hfill \mathcal{B}\phi (t)=0\end{array}$$(1)

with φ the neutron flux and ℬ the Boltzmann transport operator, reading

$$\begin{array}{cc}\hfill \mathcal{B}& =\frac{1}{v}\frac{\partial }{\partial t}+\overrightarrow{\mathrm{\Omega }}·\overrightarrow{\mathrm{\nabla }}+{\mathrm{\Sigma }}_t-\int \int {\mathrm{\Sigma }}_s{f}_{s}d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }\hfill \\ \hfill & -\frac{1}{4\pi k_{\text{eff}}}\int \int {\chi }_p{\nu }_p{\mathrm{\Sigma }}_{f}d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }\hfill \\ \hfill & -\frac{1}{4\pi k_{\text{eff}}}\sum _{j=1}^J{\int }^t{\lambda }_j{e}^{-{\lambda }_j(t-t^{\prime })}\int \int {\chi }_{d,j}{\nu }_{d,j}{\mathrm{\Sigma }}_{f}d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }d{t}^{\prime }\hfill \end{array}$$(2)

Here v is the velocity, Σ_t, Σ_s and Σ_f are the the total, scattering and fission cross sections, respectively, f_s is the scattering kernel, k_eff is the effective multiplication factor, χ is the fission energy spectrum, ν the fission multiplicity, and λ the precursors decay constant. The indices p and {d, j} refer to prompt fissions and delayed fissions precursors family j, respectively, considering J families. In the following, we will also use the notation of the total multiplicity

$$\begin{array}{c}\hfill \nu ={\nu }_p+\sum _{j=1}^J{\nu }_{d,j}\end{array}$$(3)

the delayed neutron fractions

$$\begin{array}{c}\hfill {\beta }_j=\frac{{\nu }_{d,j}}{\nu },\beta =\sum _{j=1}^J{\beta }_j\end{array}$$(4)

and the total spectrum

$$\begin{array}{c}\hfill \chi =(1-\beta ){\chi }_p+\sum _{j=1}^J{\beta }_j{\chi }_{d,j}\end{array}$$(5)

Core perturbation are introduced as periodic changes δΣ of cross-section Σ around their time-averages ⟨Σ⟩, namely

$$\begin{array}{c}\hfill \mathrm{\Sigma }(t)=⟨\mathrm{\Sigma }⟩+\delta \mathrm{\Sigma }(t)\end{array}$$(6)

Reaction yields and outgoing particles distributions (scattering laws and fission spectra) can be decomposed accordingly [17], so that for a generic kernel K we have the compact notation

$$\begin{array}{c}\hfill K(t)=⟨K⟩+\delta \left[K\right]\end{array}$$(7)

Notation δ[K] emphasizes that the perturbations are applied to the whole kernel, which typically involves a product of cross-section, multiplicity and outgoing particles distribution, rather than one of its individual factors. Thus, the Boltzmann transport operator of the perturbed system ℬ is split into a time-averaged transport operator ℬ₀ and a perturbation operator δℬ, i.e.

$$\begin{array}{c}\hfill \mathcal{B}={\mathcal{B}}_0+\delta \mathcal{B}\end{array}$$(8)

with

$$\begin{array}{cc}\hfill {\mathcal{B}}_0& =\frac{1}{v}\frac{\partial }{\partial t}+\overrightarrow{\mathrm{\Omega }}·\overrightarrow{\mathrm{\nabla }}+⟨{\mathrm{\Sigma }}_t⟩-\int \int ⟨{\mathrm{\Sigma }}_s{f}_s⟩d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }\hfill \\ \hfill & -\frac{1}{4\pi k_{\text{eff}}}\int \int ⟨{\chi }_p{\nu }_p{\mathrm{\Sigma }}_f⟩d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }\hfill \\ \hfill & -\frac{1}{4\pi k_{\text{eff}}}\sum _{j=1}^J{\int }^t{\lambda }_j{e}^{-{\lambda }_j(t-t^{\prime })}\int \int ⟨{\chi }_{d,j}{\nu }_{d,j}{\mathrm{\Sigma }}_f⟩d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }d{t}^{\prime }\hfill \end{array}$$(9)

and

$$\begin{array}{cc}\hfill \delta \mathcal{B}& =\delta {\mathrm{\Sigma }}_t-\int \int \delta \left[{\mathrm{\Sigma }}_s{f}_s\right]d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }\hfill \\ \hfill & -\frac{1}{4\pi k_{\text{eff}}}\int \int \delta \left[{\chi }_p{\nu }_p{\mathrm{\Sigma }}_f\right]d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }\hfill \\ \hfill & -\frac{1}{4\pi k_{\text{eff}}}\sum _{j=1}^J{\int }^t{\lambda }_j{e}^{-{\lambda }_j(t-t^{\prime })}\int \int \delta \left[{\chi }_{d,j}{\nu }_{d,j}{\mathrm{\Sigma }}_f\right]d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }d{t}^{\prime }\hfill \end{array}$$(10)

The restriction of operator ℬ₀ to the linear subspace of time-independent functions is the usual time-independent Boltzmann operator

$$\begin{array}{cc}\hfill {\mathcal{B}}_c& ={\mathcal{B}}_0{|}_{\frac{\partial }{\partial t}=0}=\overrightarrow{\mathrm{\Omega }}·\overrightarrow{\mathrm{\nabla }}+⟨{\mathrm{\Sigma }}_t⟩-\int \int ⟨{\mathrm{\Sigma }}_s{f}_s⟩d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }\hfill \\ \hfill & -\frac{1}{4\pi k_{\text{eff}}}\int \int ⟨\chi \nu {\mathrm{\Sigma }}_f⟩d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }\hfill \end{array}$$(11)

Solving the steady-state equation

$$\begin{array}{c}\hfill {\mathcal{B}}_c{\phi }_c=0\end{array}$$(12)

yields the critical flux φ_c and the effective multiplication factor k_eff, which is henceforth fixed. All fission sources in the time-averaged and perturbation operators are rescaled by k_eff to ensure criticality, which is a surrogate compensation taking into account the imperfect knowledge of technological specifications, material compositions and nuclear data. Criticality might be alternatively ensured by considering more realistic approaches, such as explicitly modeling stabilizing feedback from Doppler effect and thermal-hydraulics, or configuration parameters like control rod insertion or boron concentration in the moderator, but this possibility will not be considered in the present work.

The neutron flux φ is expressed as the sum of the (time-independent) critical flux φ_c and the noise δφ, i.e.

$$\begin{array}{c}\hfill \phi (t)={\phi }_c+\delta \phi (t)\end{array}$$(13)

Using this decomposition, the Boltzmann transport equation of the perturbed system in equation (1) is rewritten as

$$\begin{array}{c}\hfill {\mathcal{B}}_0{\phi }_c+{\mathcal{B}}_0\delta \phi (t)+\delta \mathcal{B}{\phi }_c(t)+\delta \mathcal{B}\delta \phi (t)=0\end{array}$$(14)

The term ℬ₀φ_c vanishes according to equation (12). Within the linear noise theory, one customarily assumes the cross-section perturbations to be small compared to their mean values, and correspondingly the noise to be small compared to the critical flux, which leads to the orthodox linearization: the term δℬδφ is considered to be a second-order contribution and is therefore neglected [1]. Finally, the Fourier transform 𝔉 is applied to obtain the linear noise equation in the frequency domain, namely,

$$\begin{array}{c}\hfill \widehat{{\mathcal{B}}_0}\widehat{\delta \phi }(\omega )=-\widehat{\delta \mathcal{B}}{\phi }_c(\omega )\end{array}$$(15)

whose solution yields the frequency-domain noise $\widehat{\delta \phi }=\mathfrak{F}(\delta \phi )$ for a given angular frequency ω. The operators occurring in equation (15) are respectively defined as

$$\begin{array}{cc}\hfill \widehat{{\mathcal{B}}_0}& =i\frac{\omega }{v}+\overrightarrow{\mathrm{\Omega }}·\overrightarrow{\mathrm{\nabla }}+⟨{\mathrm{\Sigma }}_t⟩-\int \int ⟨{\mathrm{\Sigma }}_s{f}_s⟩d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }\hfill \\ \hfill & -\frac{1}{4\pi k_{\text{eff}}}\int \int ⟨{\chi }_p{\nu }_p{\mathrm{\Sigma }}_f⟩d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }\hfill \\ \hfill & -\frac{1}{4\pi k_{\text{eff}}}\sum _{j=1}^J\frac{{\lambda }_j}{{\lambda }_j+i\omega }\int \int ⟨{\chi }_{d,j}{\nu }_{d,j}{\mathrm{\Sigma }}_f⟩d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }\hfill \end{array}$$(16)

and

$$\begin{array}{cc}\hfill \widehat{\delta \mathcal{B}}& =\mathfrak{F}(\delta {\mathrm{\Sigma }}_t)-\int \int \mathfrak{F}(\delta \left[{\mathrm{\Sigma }}_s{f}_s\right])d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }\hfill \\ \hfill & -\frac{1}{4\pi k_{\text{eff}}}\int \int \mathfrak{F}(\delta \left[{\chi }_p{\nu }_p{\mathrm{\Sigma }}_f\right])d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }\hfill \\ \hfill & -\frac{1}{4\pi k_{\text{eff}}}\sum _{j=1}^J\frac{{\lambda }_j}{{\lambda }_j+i\omega }\int \int \mathfrak{F}(\delta \left[{\chi }_{d,j}{\nu }_{d,j}{\mathrm{\Sigma }}_f\right])d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }\hfill \end{array}$$(17)

The noise equation (15) has the form of a fixed-source problem, with a noise source on the right-hand side, and a transport operator $\widehat{{\mathcal{B}}_0}$ on the left-hand side, acting on the unknown $\widehat{\delta \phi }$; all these quantities are complex-valued. The angular frequency ω is a parameter rather than a variable, and equation (15) can be thus solved independently for each value of ω. The operator $\widehat{{\mathcal{B}}_0}$ is similar to the time-independent transport operator $\widehat{{\mathcal{B}}_0}{|}_{\omega =0}={\mathcal{B}}_c$, the sole differences appearing in the complex factors iω/v, corresponding to the Fourier-transformed time-derivative, and λ_j/(λ_j + iω), corresponding to the precursor delay. In particular, for the null angular frequency ω = 0, we recover the critical operator ℬ_c.

Useful insight on the behavior of equation (15) is gained from point-kinetics analysis, i.e., by integrating out all phase space variables except ω [1]. Analytical expressions can be derived for the so called fundamental transfer function between the noise source and the noise field: for illustration, the resulting shape is displayed in Figure 1. Two cut-off frequencies appear, namely,

$$\begin{array}{c}\hfill {\omega }_{\text{LF}}={\lambda }_{\text{eff}},{\omega }_{\text{HF}}={\lambda }_{\text{eff}}+{\beta }_{\text{eff}}/{\mathrm{\Lambda }}_{\text{eff}}\approx {\beta }_{\text{eff}}/{\mathrm{\Lambda }}_{\text{eff}}\end{array}$$(18)

with Λ_eff the effective mean neutron generation time (or equivalently effective mean neutron lifetime, since criticality is enforced), β_eff the effective delayed fraction and λ_eff the effective mean decay constant [1, 16]. We indicate with f_LF = ω_LF/2π and f_HF = ω_HF/2π the associated frequencies. The cut-off frequencies define three regions: in the low-frequency regime (ω ≤ ω_LF), the noise source tends to be amplified; in the plateau region (ω_LF ≤ ω ≤ ω_HF) the noise tends to propagate and the amplification factor depends very little on the angular frequency; in the high-frequency regime (ω_HF ≤ ω), the noise tends to be damped. In most cases of interest in reactor physics applications, perturbations occur at frequencies lying in the plateau region.

3. Monte Carlo solvers for neutron noise

Over the last few years, a strategy has been proposed and progressively refined in order to solve the linear noise equation (15) by Monte Carlo methods, with two main variants depending on how the Fourier-transformed time-derivative is interpreted [14, 15]. In both cases, the algorithm starts by solving the usual steady-state equation (12) by power iteration, in order to estimate φ_c, which is required to determine the noise source $\widehat{\delta \mathcal{B}}{\phi }_c(\omega )$. Observe that the critical equation includes time-averaged quantities, from the original unperturbed system as well as from any non-zero-average perturbation: this therefore corresponds to a new steady-state [16], whose k_eff may be different from the one of the unperturbed system; in practice, most noise solvers neglect this difference and consider exclusively the unperturbed system. Once the power iteration has achieved convergence, some additional cycles are computed to tally the complex noise source $\widehat{\delta \mathcal{B}}{\phi }_c$, which depends on the specific noise problem. The noise source is then injected into a fixed-source computation to solve equation (15).

Complex quantities are handled in Monte Carlo by allowing each particle to carry complex weights, instead of the usual statistical weights, i.e. real positive weights [14]. The complex factors occurring in the transport operator $\widehat{{\mathcal{B}}_0}$ are used as weight multipliers during the transport routines. The tallies associated to the sought response functions are multiplied by the complex weights and are thus complex.

3.1. Flight and collision-based methods

Two methods are currently available to interpret the Fourier-transformed time-derivative iω/v. The flight-based method [14] associates this term to the total cross-section ⟨Σ_t⟩ by introducing a modified complex-valued total cross-section ${\tilde{\mathrm{\Sigma }}}_t^{\text{fb}}=⟨{\mathrm{\Sigma }}_t⟩+i\omega /v$, which is given practical meaning as follows. The flight length s is sampled using only the total cross-section ⟨Σ_t⟩. If the cross-section is piece-wise constant, as customary, this corresponds to sample from the density f(s)=⟨Σ_t⟩exp(−⟨Σ_t⟩s), using e.g. surface tracking or delta tracking. Then, the particle weight is multiplied by exp(−iωs/v) to take into account the modified total cross-section.

The collision-based method [15] introduces an additional copy collision event, with a real and positive cross-section ${\mathrm{\Sigma }}_{\eta }=\eta \frac{\omega }{v}$, where the positive factor η can be chosen freely, and associated complex multiplicity (η − i)/η. The total cross-section is increased accordingly, i.e., ${\tilde{\mathrm{\Sigma }}}_t^{\text{cb}}=⟨{\mathrm{\Sigma }}_t⟩+\eta \frac{\omega }{v}$, and is still real-valued. The phase-space coordinates of the particles emitted from copy collision events are identical to those of the incident particles, and their weight is multiplied by (η − i)/η to take into account multiplicity. Taking η = 1 generally yields a good figure of merit for source frequencies in the plateau region. Previous works have shown that increasing η by several order of magnitudes might improve convergence properties at high frequencies.

In both methods, the weights of the particles emitted by delayed fission event associated to precursor family j are multiplied by λ_j/(λ_j + iω) to take into account the Fourier-transformed delayed fission operator.

3.2. Russian roulette

To ensure termination of the computation, in particular when using forced fission and implicit capture, a Russian roulette protocol must be devised for particles carrying complex weights w ∈ ℂ. Most algorithms use the two-part roulette introduced by Yamamoto [14], which is performed separately on the real ℜw and the imaginary ℑw parts of the weight. In detail, the weight w′∈ℂ of the particle after the roulette is selected as follows. Let w₀ ≤ 1 be the roulette threshold. If |ℜw| ≤ w₀, the real part of w′ is sg(ℜw) with probability |ℜw|, or zero with complementary probability; otherwise, the real part of w′ is ℜw. Similarly, if |ℑw| ≤ w₀, the imaginary part of w′ is sg(ℑw) with probability |ℑw|, or zero with complementary probability; otherwise, the imaginary part of w′ is ℑw. The particle is only suppressed if neither its real or imaginary part survive roulette, i.e. if w′=0.

In a previous work, we have introduced an alternative roulette procedure based on the norm of the complex weight [19]. If |w| ≤ w₀, the particle survives with probability |w| with a new weight w′=w/|w|, or is suppressed with complementary probability. This ensures that, in the case of survival, the weight argument is conserved. This idea was initially proposed to simplify the theoretical analysis of the Monte Carlo algorithms. We will show in the following that this choice also improves the convergence of the noise simulations.

3.3. Weight cancellation

In most applications of noise calculations using Monte Carlo methods, weight cancellation has proven crucial [14, 17, 18]. This population control technique relies on multiple particles weights being summed together, resulting in the cancellation of opposite weight contributions and fewer particles. This in turn requires the use of rendez-vous points for the particle: in practice, this is achieved by restructuring the transport routine and decomposing particle histories over successive generations; weight cancellation is simultaneously applied to all particles reaching the end of each generation.

The implementation proposed by Yamamoto for noise problems is based on a binning procedure [14]: a fictitious mesh is superposed to the model geometry and used to sum the weights of the particles falling within each cell of the mesh. Particles are then re-sampled from those in the cell, and inherit a fraction of the total weight of the cell. All particles from a given bin emerge with the same weight. This procedure is approximate, in that it introduces discretizations errors due to the bin-averaging procedure. The bias can be reduced by reducing the bin size, at the expense of a loss of efficiency due to fewer particles being available for cancellation in the same bin. Other cancellation methods have been proposed; for an overview, we refer the reader to e.g. Ref. [16]. Recently, an exact regional cancellation method for noise problems has been also proposed [17], inspired by a technique originally introduced by Booth [20] and refined by Belanger et al. [18]. However, the current computational cost of these methods overcomes their precision gain with respect to the simpler (albeit biased) mesh-based procedures [17].

The definition of the generation used as a rendez-vous for particle histories is in principle arbitrary. As mentioned, weight cancellation nonetheless needs to aggregate sufficiently large groups of particles to be efficient. A strong argument can be made in favor of using fission generations: particles emerging from fission events are sent to the next generation, while any other secondary particles are recycled into the current generation. Indeed, the emission of particles from fission events being typically isotropic and independent of the incident neutron energy, weight cancellation techniques can merely discriminate particles by their positions alone. For this reason, if a collision-based algorithm is used for the simulation of particle histories, the copy reaction is required to be recycled within the same generation, instead of being added to the particles banked for the next generation, since the scattering kernel of the copy reaction is singular¹.

3.4. Review of known limitations of the current algorithms

Monte Carlo algorithms for neutron noise have been implemented in several codes and verified against each other or deterministic solvers [4, 15, 16, 21]. Furthermore, preliminary validation results have been obtained using neutron noise experiments performed in research reactors [22]. In most cases, convergence issues in the noise solvers have been raised in the literature, ranging from mild to severe. In this section, convergence will be used in a loose sense, meaning that simulations were successfully completed and provided coherent results; a more rigorous definition will be provided later on. An overview of such problems has been provided by Rouchon [16].

Rouchon’s investigations were conducted for collision and flight-based strategies, with pseudo-analog ² and non-analog games: either sampling one collision channel at a time without Russian roulette, or applying implicit capture, forced fission and copy where relevant, combined with Russian roulette. For Russian roulette, the two-part implementation was used. Weight cancellation was tested exclusively for the flight-based method. Noise problems at frequencies in the range from about the lower cut-off up to well beyond the upper cut-off were examined.

Both collision and flight-based methods were found to converge for noise source frequencies in the plateau region and slightly above the upper cut-off frequency, without the need to enforce weight cancellation. Outside this region, the flight-based method required weight cancellation to converge. Alternatively, the collision-based method could be made to converge, provided that a pseudo-analog game was used. Besides, increasing the factor η for collision-based methods was observed to be necessary, in addition to using a pseudo-analog game, to obtain satisfactory results for higher frequencies.

The robustness and reliability of Monte Carlo methods for noise problems is key to their implementation in production codes for real-world applications. In this respect, precisely determining convergence properties of these algorithms and firmly establishing how they behave under the combined effects of different population-control techniques is mandatory. Inspired by these considerations, in the following sections we propose a broad theoretical framework that makes it possible to formulate these questions. Although some forms of the noise source $-\widehat{\delta \mathcal{B}}{\phi }_c$, such as mechanical vibrations, are known to introduce additional complications for the noise solvers [16, 17], the convergence properties of the Monte Carlo algorithms are independent of the source and this term will not be addressed in the following.

4. Theory of complex-valued Monte Carlo

Before considering complex-valued Monte Carlo games, we begin by reminding the basic formulation of standard fixed-source Monte Carlo games, by building on the general principles formulated in Refs. [12] and [13]. Let us consider a phase-space 𝒟, whose points are denoted by x, the space L¹(𝒟) of summable functions on 𝒟, and a positive integral operator 𝒦 on L¹(𝒟) with non-negative kernel K:

$$\begin{array}{c}\hfill \mathcal{K}\psi (x)={\int }_{\mathcal{D}}K(x,x^{\prime })\psi (x^{\prime })d{x}^{\prime }\end{array}$$(19)

Consider a non-negative source q and a bounded detector response function h. We want to determine the response

$$\begin{array}{c}\hfill R={\int }_{\mathcal{D}}h(x)\psi (x)dx\end{array}$$(20)

associated to an unknown function ψ, which is the solution of the Fredholm integral equation of the second kind

$$\begin{array}{c}\hfill \psi =\mathcal{K}\psi +q\end{array}$$(21)

The existence of solutions to equation (21) is characterized by the spectrum of operator 𝒦. In particular, conditions may be formulated on the spectral radius, which is the supremum of the norm of the elements of the spectrum of the operator. Intuitively, the spectral radius can be assimilated with the modulus of the dominant eigenvalues of this operator; although this equality is admitted to hold in reactor physics problems, to the best of our knowledge it has not been demonstrated in all generality.

Assuming that the spectral radius of 𝒦 is striclty smaller than 1, the Neumann series of the operator 𝒦, namely

$$\begin{array}{c}\hfill \sum _{n=0}^{+\infty }{\mathcal{K}}^n\end{array}$$(22)

converges in the operator norm topology, and the function

$$\begin{array}{c}\hfill \psi =\sum _{n=0}^{+\infty }{\mathcal{K}}^{n}q\end{array}$$(23)

is a non-negative solution of equation (21). Moreover, the response is given by

$$\begin{array}{c}\hfill R=\sum _{n=0}^{+\infty }{\int }_{\mathcal{D}}h(x){\mathcal{K}}^{n}q(x)dx\end{array}$$(24)

In order to estimate the response R by Monte Carlo games, a crucial remark is that, if a random variable x taking values in 𝒟 is sampled from a probability density f, if K(⋅, x) is a probability density function, and y is sampled from K(⋅, x), then the probability density of y is 𝒦f. This fact stems from the definition of a conditional probability density. The main idea of fixed-source Monte Carlo games is to sample a particle starting position x₀ according to density q, then sequentially position x_n + 1 from density K(⋅, x_n). By tallying the sum of terms h(x_n), we build the Wasow random variable [13]

$$\begin{array}{c}\hfill \widehat{R}=\sum _{n=0}^{+\infty }h(x_n)\end{array}$$(25)

which can be shown to be an unbiased estimator of the sought response R.

In practice, K(⋅, x) is rarely a probability density function. Notably, if such was the case for any x, then 𝒦 would be an isometry on the cone of positive functions and therefore its spectral radius would be 1. The lack of normalization can be compensated in several ways. The first method is to use multiplicity: from a particle at x, sample ν particles, where ν is a random variable taking integer values and expectation A(x)=∫_𝒟K(x′,x)dx′, each from the probability density K(y, x)/A(x). The second method is to attribute each particle a statistical weight: from a particle at x with weight w, sample a particle from the probability density K(y, x)/A(x), adjust the particle weight to w′=A(x)w, and tally the sum of terms wh(x). More generally, a combination of the above may be used.

The normalization of the source q can be dealt with using either of the above methods. Alternatively, the starting position can be sampled from the normalized source, and the tally estimator multiplied by the source norm at the end of the computation, thanks to linearity.

These procedures can be immediately extended to the case of q(y)dy and K(y, x)dy being positive finite measure, i.e. probability measures up to normalization. In particular, this allows handling Dirac delta functions in the kernels.

4.1. Complex-weighted Monte Carlo games

To account for the idea of complex particle weights, we treat the complex weight as an additional phase-space coordinate. For the sake of simplicity, we start by considering analog complex-weighted games, which we define as games where weights w are restricted to be of modulus 1. Consequently, a complex weight w can be described by its phase³θ, which belongs to 𝒮₁ = ℝ/2πℤ, the real numbers modulo 2π. Let us consider the extended phase-space $\tilde{\mathcal{D}}=\mathcal{D}\times {\mathcal{S}}_1$ whose points will be denoted $\tilde{x}=(x,\theta )$.

In this space, multiplication by a unitary complex number e^iμ amounts to a shift in phase, represented by the kernel

$$\begin{array}{c}\hfill Z_{\mu }(\theta ,{\theta }^{\prime })=\delta (\theta -{\theta }^{\prime }-\mu )\end{array}$$(26)

where δ is the Dirac delta. This choice is not unique, but it is optimal in a sense, as discussed in Appendix A. Let us consider the kernel K to be a complex-weighted sum of non-negative kernels, i.e.

$$\begin{array}{c}\hfill K(x,x^{\prime })=\sum _j{e}^{i{\mu }_j(x,x^{\prime })}{K}_j(x,x^{\prime })\end{array}$$(27)

where μ_j(x, x′) are real phases, that may depend on both in- and outgoing coordinates. We introduce the extended kernel

$$\begin{array}{c}\hfill \tilde{K}(\tilde{x},{\tilde{x}}^{\prime })=\sum _j{Z}_{{\mu }_j(x,x^{\prime })}(\theta ,{\theta }^{\prime })K_j(x,x^{\prime })\end{array}$$(28)

the extended detector function

$$\begin{array}{c}\hfill \tilde{h}(\tilde{x})=e^{i\theta }h(x)\end{array}$$(29)

and the extended source $\tilde{q}$

$$\begin{array}{c}\hfill \tilde{q}(\tilde{x})=|q(x)|\delta (\theta -arg(q(x)))\end{array}$$(30)

We can now write the extended Fredholm equation

$$\begin{array}{c}\hfill \tilde{\psi }=\tilde{\mathcal{K}}\tilde{\psi }+\tilde{q}\end{array}$$(31)

where the integral operator $\tilde{\mathcal{K}}$ has the kernel $\tilde{K}$, and the extended response

$$\begin{array}{c}\hfill \tilde{R}={\int }_{\tilde{\mathcal{D}}}\tilde{h}(\tilde{x})\tilde{\psi }(\tilde{x})d\tilde{x}\end{array}$$(32)

If the spectral radius of $\tilde{\mathcal{K}}$ is stricly smaller than 1, a solution of equation (31) exists and is given by

$$\begin{array}{c}\hfill \tilde{\psi }=\sum _{n=0}^{+\infty }{\tilde{\mathcal{K}}}^n\tilde{q}\end{array}$$(33)

Since equation (31) only involves non-negative quantities, it is suitable for the standard Monte Carlo procedure. Moreover, under the previous assumptions, it is straightforward to show that $\tilde{R}=R$, so that this Monte Carlo game can be used to estimate the sought response of the original complex-valued equation. In the following, we will denote all the quantities pertaining to equation (21) and operator 𝒦 as physical, and all the quantities pertaining to equation (31) and operator $\tilde{\mathcal{K}}$ as extended.

It is worth noting that this extension process can also be applied to take into account statistical weights, i.e. the usual, real positive weights, but yields no new information compared to previous formalizations. This can also be generalized to non-analog complex weighted games, in which the norm of complex weights is another additional dimension that is equivalent to statistical weights. This approach is described in Appendix B. Negative-weighted games are also a special case of the above, where e^{iμ_j(x, x′)} = ±1; in other words, μ_j(x, x′) are either 0 or π.

4.2. On the notion of convergence

In order to assess the convergence pitfalls of complex-weighted Monte Carlo games, it is necessary to agree on the meaning of convergence. We define a Monte Carlo game to be convergent if the tally estimator series given in equation (25), as a series of random variables, is convergent in mean and is an unbiased estimator of the desired tally in equation (20), regardless of the summable source q. According to general theorems, this implies convergence in probability.

A tightly related notion is that of sub-criticality [12], which is defined as the spectral radius of operator 𝒦 being strictly lower than 1. As exposed above, when sampling successively from the (potentially not-normalized) density K(⋅, x), sub-criticality ensures both that the tally estimator series converges in mean and is unbiased, for any bounded response function h. It is therefore a sufficient condition for convergence. Little can be said about convergence when sub-criticality is not ensured without additional assumptions⁴; counter-examples of critical or supercritical operators that produce convergent games for specific tallies may be designed⁵.

The notion of convergence is distinct from both that of feasibility, that Lux and Koblinger define as the Monte Carlo game ending with probability 1 [12], and that of statistical convergence with an increasing number of particles, which would be ensured by the central limit theorem provided finite second moment of the tally estimator.

4.3. Twofold convergence criterion

Fixed-source problems are typically formulated using an integro-differential form, which must be cast into an integral form to fit equation (21). Consider the fixed-source problem

$$\begin{array}{c}\hfill \mathcal{B}\phi =Q\end{array}$$(34)

where the Boltzmann operator ℬ could be either $\widehat{{\mathcal{B}}_0}$ from equation (16) or ℬ_c from equation (11). We introduce the collision density ψ = Σφ, where Σ is the total macroscopic cross-section of the problem under consideration, and we decompose the operator ℬ into a differential leakage operator $\mathcal{L}={\mathrm{\Sigma }}^{-1}\overrightarrow{\mathrm{\Omega }}·\overrightarrow{\mathrm{\nabla }}+1$, an integral scattering operator 𝒮, and an integral fission operator ℱ, each factorized by the total cross-section, i.e.

$$\begin{array}{c}\hfill \mathcal{B}=(\mathcal{L}-\mathcal{S}-\mathcal{F})\mathrm{\Sigma }\end{array}$$(35)

Equation (34) can be integrated along straight lines with direction $\overrightarrow{\mathrm{\Omega }}$, which are the characteristic curves of the leakage operator ℒ. This is formally equivalent to inverting ℒ, whose inverse 𝒯 = ℒ⁻¹ is the integral flight operator [12]. We can then write

$$\begin{array}{c}\hfill \psi =\mathcal{T}(\mathcal{S}+\mathcal{F})\psi +\mathcal{T}Q\end{array}$$(36)

which is a Fredholm integral equation of the second kind, with first collision source 𝒯Q and operator 𝒯(𝒮 + ℱ).

Generations can be formally accounted for by factorizing the operator ensuring transition between generations, here ℱ for fission generations. The missing step is to move 𝒯𝒮ψ to the left-hand side, and invert the operator (ℐ − 𝒯𝒮), where ℐ is the identity operator. This leads to

$$\begin{array}{c}\hfill \psi ={(\mathcal{I}-\mathcal{T}\mathcal{S})}^{-1}\mathcal{T}\mathcal{F}\psi +{(\mathcal{I}-\mathcal{T}\mathcal{S})}^{-1}\mathcal{T}Q\end{array}$$(37)

Equation (37) can be expressed as a Fredholm integral equation of the second kind like equation (21), with operator

$$\begin{array}{c}\hfill \mathcal{K}={(\mathcal{I}-\mathcal{T}\mathcal{S})}^{-1}\mathcal{T}\mathcal{F}\end{array}$$(38)

and first-absorption source

$$\begin{array}{c}\hfill q={(\mathcal{I}-\mathcal{T}\mathcal{S})}^{-1}\mathcal{T}Q\end{array}$$(39)

To assess whether the operator (ℐ − 𝒯𝒮) can indeed be inverted, we consider the following c-eigenvalue problem:

$$\begin{array}{c}\hfill c{\psi }_c=\mathcal{T}\mathcal{S}{\psi }_c\end{array}$$(40)

If the spectral radius of operator 𝒯𝒮 is strictly smaller than 1, then (ℐ − 𝒯𝒮) is invertible and the Neumann series

$$\begin{array}{c}\hfill {(\mathcal{I}-\mathcal{T}\mathcal{S})}^{-1}=\sum _{j=0}^{+\infty }{(\mathcal{T}\mathcal{S})}^j\end{array}$$(41)

holds. The operator 𝒦 can then be written as

$$\begin{array}{c}\hfill \mathcal{K}=\sum _{j=0}^{+\infty }{(\mathcal{T}\mathcal{S})}^j\mathcal{T}\mathcal{F}\end{array}$$(42)

which provides a formulation as a positive integral operator. A solution of equation (37) exists and is given by equation (23), provided that the spectral radius of 𝒦 is also strictly smaller than 1, which is expressed through the eigenvalue equation

$$\begin{array}{c}\hfill k{\psi }_k=\mathcal{K}{\psi }_k\end{array}$$(43)

Equation (43) is equivalent to

$$\begin{array}{c}\hfill {\psi }_k=\mathcal{T}(\mathcal{S}+\frac{1}{k}\mathcal{F}){\psi }_k\end{array}$$(44)

where we recover the usual k-eigenvalue problem.

Let us denote by 𝔎 and 𝔠 the spectral radii of operators 𝒦 and 𝒯𝒮, respectively. Considering the above, if the following criterion hold,

$$\begin{array}{c}\hfill \mathfrak{c}<1,\mathfrak{K}<1\end{array}$$(45)

then equation (23) gives a solution of the Fredholm integral equation in equation (37), equation (24) holds and application of 𝒦 can be computed by repeated application of 𝒯𝒮, according to equation (42). In particular, if the source Q and the kernels of operators 𝒯, 𝒮 and ℱ are real non-negative, a Monte Carlo game derived from operator 𝒦 is convergent.

Criterion (45) can be interpreted from a computational point of view. The procedure of inverting the operator (ℐ − 𝒯𝒮) as the Neumann series in equation (41) is the computation of one generation, composed of successive flights and scattering events until the absorption of the particle. Therefore, the condition 𝔠 <  1 corresponds to sub-criticality inside a generation. This condition is always satisfied in usual reactor physics problems. The condition 𝔎 <  1 has its usual interpretation of overall sub-criticality [12], i.e. configurations where multiplication is dominated by captures and leakage.

4.4. Noise eigenvalues

The previous considerations can be applied to Monte Carlo computations of the linear noise equation (15) in the frequency domain, using either the physical form of the operator $\widehat{{\mathcal{B}}_0}$ given in equation (16), or its extended form as described in Section 4.1. The former case yields the physical eigenvalues k and c and radii 𝔎 and 𝔠. The latter yields their extended counterpart, noted with tildes, e.g. $\tilde{\mathfrak{K}}$ and $\tilde{\mathfrak{c}}$. Considering a physical eigenvalue associated to an eigenvector ψ(x), one can show that $\tilde{\psi }(\tilde{x})=\psi (x)e^{-i\theta }$ is an eigenvector of the corresponding extended operator for the same eigenvalue. Consequently, the physical eigenvalues are a subset of the extended eigenvalues. Any other extended eigenvalue is called non-physical, resulting from the extension process and bearing no physically meaningful information. More generally, one can show that the physical spectrum is a subset of the extended spectrum, but we omit the proof for the sake of conciseness. A consequence of this inclusion property is the relation between the spectral radii, namely

$$\begin{array}{c}\hfill \mathfrak{c}\le \tilde{\mathfrak{c}},\mathfrak{K}\le \tilde{\mathfrak{K}}\end{array}$$(46)

Since the Monte Carlo game is performed on the extended noise equation, the convergence criterion (45) is replaced by

$$\begin{array}{c}\hfill \tilde{\mathfrak{c}}<1,\tilde{\mathfrak{K}}<1\end{array}$$(47)

Note also that both 𝔎 and $\tilde{\mathfrak{K}}$ of the noise equation differ from k_eff, the dominant eigenvalue of the critical equation, which is already taken into account in the operator expression equation (16). For completeness, the expressions of operators 𝒯, 𝒮 and ℱ for neutron noise are provided in Appendix C.

4.5. Weight cancellation

A complete description of weight cancellation is currently missing in the literature. We therefore resort to Belanger et al.’s linear cancellation model, which may provide useful insights despite the drastic simplifications that it relies on [18].

Concisely, the linear cancellation model assumes that weight cancellation can be represented as a linear operator 𝒜. A perfect weight cancellation operator is a projector on the linear span of the physical eigenvectors, whose kernel is the linear span of the non-physical eigenvectors. All physical eigenvalues are unmodified by the weight cancellation, while non-physical eigenvalues are cast to 0. Therefore, the non-null eigenvalues of the extended operator with perfect weight cancellation $\mathcal{A}\tilde{\mathcal{K}}$ are exactly the eigenvalues of the physical operator 𝒦. Weight cancellation is only applied at the end of a generation: the set of $\tilde{c}$-eigenvalues remains unchanged. Therefore, with weight cancellation between fission generations, the convergence criterion in equation (47) becomes

$$\begin{array}{c}\hfill \tilde{\mathfrak{c}}<1,\mathfrak{K}<1\end{array}$$(48)

Belanger et al.’s linear cancellation model also provides a description of imperfect weight cancellation. Such an effect arises in particular because of the finite number of particle in the simulation. Given a perfect weight cancellation operator 𝒜, the corresponding imperfect weight cancellation operator 𝒜^(α) with efficiency α ∈ [0, 1] is

$$\begin{array}{c}\hfill {\mathcal{A}}^{(\alpha )}=\alpha \mathcal{A}+(1-\alpha )\mathcal{I}\end{array}$$(49)

The spectral radius of the extended operator with imperfect weight cancellation ${\mathcal{A}}^{(\alpha )}\tilde{\mathcal{K}}$ is

$$\begin{array}{c}\hfill {\tilde{\mathfrak{K}}}^{(\alpha )}=max(\mathfrak{K},(1-\alpha )\tilde{\mathfrak{K}})\end{array}$$(50)

and the corresponding criterion is

$$\begin{array}{c}\hfill \tilde{\mathfrak{c}}<1,{\tilde{\mathfrak{K}}}^{(\alpha )}<1\end{array}$$(51)

Criterion (47) and (48) are retrieved for α = 0 and α = 1 respectively. Since ${\tilde{\mathfrak{K}}}^{(\alpha )}$ is a decreasing function of α, this criterion may only be satisfied for a sufficient weight cancellation efficiency α. This may for instance be achieved by simulating a sufficient number of particles.

One might also consider weight cancellation inside a generation. However, this does not seems practically feasible: it would need to differentiate the particle on all the dimensions, instead of the spatial position for the fission source alone, which would require a tremendous amount of particles for cancellation to be effective. Besides, performing weight cancellation at each collision would result in an unreasonable computational cost. Alternatively, one may consider cancellation every n collision, which reduce the computational cost but potentially limit the effects in case of imperfect cancellation, since the relevant spectral radius would then be $max({\mathfrak{c}}^n,(1-\alpha ){\tilde{\mathfrak{c}}}^n)$ and therefore require even more particles to achieve sufficient cancellation. These strategies can easily be included in our formalism and provide a corresponding criterion but are not investigated in the following.

We emphasize that all the given criteria are relevant for a Fredholm equation of the second kind, i.e. a fixed-source problem. The dominant solution to a Fredholm equation of the first kind, i.e. a critical problem, is usually computed by Monte Carlo in a related manner, relying on power iteration. The convergence criterion for such games is quite different from those given above. The condition $\tilde{\mathfrak{c}}<1$ remains. However, since the physical eigenvector is sought, one must ensure that its eigenvalue is dominant, which is typically not the case without weight cancellation. Weight cancellation is therefore mandatory, along with the condition $(1-\alpha )\tilde{\mathfrak{K}}<\mathfrak{K}$. Beside, the usual condition 𝔡 <  1 on the dominance ratio 𝔡 of operator ${\mathcal{A}}^{(\alpha )}\tilde{\mathcal{K}}$ applies. Therefore, the criterion for power iteration is

$$\begin{array}{c}\hfill \tilde{\mathfrak{c}}<1,(1-\alpha )\tilde{\mathfrak{K}}<\mathfrak{K},\mathfrak{d}<1\end{array}$$(52)

4.6. Convergence regimes

Based on the previous considerations and definitions of Section 4.2, we can infer three possible regimes for the fixed-source, complex-weighted Monte Carlo algorithm used to solve the linear neutron noise equation. Let again $\mathfrak{c},\tilde{\mathfrak{c}},\mathfrak{K},\tilde{\mathfrak{K}}$ be the physical and extended eigenvalues of the operators 𝒯𝒮 and 𝒦 = (ℐ − 𝒯𝒮)⁻¹𝒯ℱ (without weight cancellation), and remember that the detector response function h is assumed to be bounded. The three convergence regime are the following:

1. Convergent. If $\mathfrak{K}\le \tilde{\mathfrak{K}}<1$ and $\mathfrak{c}\le \tilde{\mathfrak{c}}<1$, both computation of generations and the overall computation are sub-critical, therefore the Monte Carlo game is convergent.

2. Conditionally convergent. If $\mathfrak{K}<1\le \tilde{\mathfrak{K}}$ and $\mathfrak{c}\le \tilde{\mathfrak{c}}<1$, then the computation of each generation is sub-critical. However, the overall computation is not, unless weight cancellation is enforced, with a sufficiently high efficiency. Without sufficient weight cancellation, this regime is peculiar: the series in equation (33) is divergent in $L^1(\tilde{\mathcal{D}})$. Empirical evidence suggests that the second moment of the tally is not finite, and even the first moment or the random variable altogether may not be defined. However, we are unable to provide a general characterization of the convergence properties in this regime.

3. Divergent. If the previous conditions are not met, the computation is not subcritical, even with weight cancellation. This can arise for two reasons. If $1<\mathfrak{K}\le \tilde{\mathfrak{K}}$ and $\mathfrak{c}\le \tilde{\mathfrak{c}}<1$, generations are sub-critical, but the overall simulation is not, even with weight cancellation. If $1<\tilde{\mathfrak{c}}$, regardless of the other eigenvalues, generations are super-critical. In general, the Monte Carlo game does not converge in this regime.

5. Application to simple benchmark problems

We probe now the theoretical framework developed in the previous section on two benchmark problems where exact reference solutions can be established and used for comparison.

The spectral radii are given in the spaces L^p (𝒟) for physical values and $L^p(\tilde{\mathcal{D}})$ for extended values, for any 1 <  p <  ∞. The relevant operators can be fully diagonalized using Fourier series, i.e. we obtain a countable basis of the functional space with eigenvectors. This ensures that the whole point spectrum has been found, and that the spectrum is exactly the closure of the point spectrum, so that the spectral radius is indeed the supremum of the modulus of eigenvalues. Moreover, this value is independent of the choice of p. The case L¹, although the natural space to formulate Monte Carlo, is particularly arduous, because of the classical issue of the trigonometric series not necessarily converging in the L¹ norm. Therefore, the expressions provided are only lower bounds of the spectral radii. This is assumed to be a technicality for reactor physics applications, and is not investigated in the following.

5.1. Problem 1: Infinite homogeneous medium

We start by considering the simple case of an infinite homogeneous medium, with a single group of energy and a single family of precursors (i.e. we drop the index j = 1). For clarity, since we are only concerned with operator $\widehat{{\mathcal{B}}_0}$, we also drop the brackets in the time-averaged quantities, e.g. writing νΣ_f rather than ⟨νΣ_f⟩. Using the notation

$$\begin{array}{c}\hfill \mathrm{{\rm Y}}=1-\beta +\frac{\beta \lambda }{\lambda +i\omega }\end{array}$$(53)

for the factor acting on the fission operator, the operator $\widehat{{\mathcal{B}}_0}$ reads

$$\begin{array}{c}\hfill \widehat{{\mathcal{B}}_0}=\frac{i\omega }{v}+{\mathrm{\Sigma }}_t-{\mathrm{\Sigma }}_s-\frac{\mathrm{{\rm Y}}}{k_{\text{eff}}}\nu {\mathrm{\Sigma }}_f\end{array}$$(54)

As customary, we set Σ_a = Σ_t − Σ_s.

5.1.1. Criticality

Let us first consider the critical equation, i.e. ω = 0. It is convenient to introduce the parameters

$$\begin{array}{c}\hfill k_{\infty }=\frac{\nu {\mathrm{\Sigma }}_f}{{\mathrm{\Sigma }}_a},c_{\infty }=\frac{{\mathrm{\Sigma }}_s}{{\mathrm{\Sigma }}_t}\end{array}$$(55)

The effective multiplication factor k_eff is the dominant k-eigenvalue of this unperturbed configuration. We define accordingly c_eff the dominant c-eigenvalue of the same configuration. In the absence of leakage, we have

$$\begin{array}{c}\hfill k_{\text{eff}}=k_{\infty },c_{\text{eff}}=c_{\infty }\end{array}$$(56)

Criticality is enforced by dividing fission sources by k_eff in equation (54).

5.1.2. Noise eigenvalues

We introduce the extended delay factor

$$\begin{array}{c}\hfill \widetilde{\mathrm{{\rm Y}}}=1-\beta +\frac{\beta \lambda }{|\lambda +i\omega |}\end{array}$$(57)

which satisfies $|\mathrm{{\rm Y}}|\le \widetilde{\mathrm{{\rm Y}}}$. For the flight-based method, we obtain the dominant eigenvalues

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathfrak{c}}_1^{\text{fb}}=\left|\frac{c_{\infty }}{1+i\omega \tau }\right|,& {\mathfrak{K}}_1^{\text{fb}}=\left|\frac{\mathrm{{\rm Y}}}{1+i\omega \mathrm{\Lambda }}\right|\hfill \\ \hfill {\tilde{\mathfrak{c}}}_1^{\text{fb}}=c_{\infty },& {\tilde{\mathfrak{K}}}_1^{\text{fb}}=\widetilde{\mathrm{{\rm Y}}}\hfill \end{array}\end{array}$$(58)

where we have defined the mean neutron generation time (or equivalently mean neutron lifetime) of the critical system

$$\begin{array}{c}\hfill \mathrm{\Lambda }=\frac{1}{v\nu {\mathrm{\Sigma }}_f/k_{\text{eff}}}=\frac{1}{v{\mathrm{\Sigma }}_a}\end{array}$$(59)

and the mean time between collisions

$$\begin{array}{c}\hfill \tau =\frac{1}{v{\mathrm{\Sigma }}_t}\end{array}$$(60)

For the collision-based method, we obtain the dominant eigenvalues

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathfrak{c}}_1^{\text{cb}}=\left|\frac{c_{\infty }+(\eta -i)\omega \tau }{1+\eta \omega \tau }\right|,& {\mathfrak{K}}_1^{\text{cb}}=\left|\frac{\mathrm{{\rm Y}}}{1+i\omega \mathrm{\Lambda }}\right|\hfill \\ \hfill {\tilde{\mathfrak{c}}}_1^{\text{cb}}=\frac{c_{\infty }+\sqrt{1+{\eta }^2}\omega \tau }{1+\eta \omega \tau },& {\tilde{\mathfrak{K}}}_1^{\text{cb}}=\frac{\widetilde{\mathrm{{\rm Y}}}}{1-\kappa \omega \mathrm{\Lambda }}\hfill \end{array}\end{array}$$(61)

where for the sake of conciseness we set

$$\begin{array}{c}\hfill \kappa =\sqrt{1+{\eta }^2}-\eta \end{array}$$(62)

The flight-based and collision-based physical k-eigenvalue equations correspond to two formulations of the same problem, so that we can set

$$\begin{array}{c}\hfill {\mathfrak{K}}_1={\mathfrak{K}}_1^{\text{fb}}={\mathfrak{K}}_1^{\text{cb}}\end{array}$$(63)

5.1.3. Analysis of the eigenvalues

The theoretical spectral radii given in equations (58) and (61) are illustrated graphically in Figure 2, using the numerical values provided in Table 1.

Using the flight-based method, for all frequencies we have $\tilde{\mathfrak{c}}<1$ and $\tilde{\mathfrak{K}}<1$: there is no issue of convergence inside a given generation and the algorithm should overall always be convergent. Note however that $\tilde{\mathfrak{K}}\approx 1$ at low frequency, which corresponds to a fixed-source problem in nearly-critical conditions, which is known to be very hard to converge. Therefore, we expect the algorithm to perform fairly bad at low frequency, unless special precautions are adopted⁶.

Using the collision-based method, different convergence regimes are observed. Although not visible in the figures, at frequencies lower than a frequency $f_{\mathrm{VLF}}^{\text{cb}}={\omega }_{\mathrm{VLF}}^{\text{cb}/2}\pi$, the collision-based method should be conditionally convergent. However, this only concerns frequencies several orders of magnitude lower than the low frequency cut-off f_LF, e.g. here ${\omega }_{\mathrm{VLF}}^{\text{cb}}=2.5\times {10}^{-6}\mathrm{rad}{\mathrm{s}}^{-1}$ which correspond to a period about a month, and the amplification factor of the fundamental transfer function is about 4.5 × 10⁶ - for these frequencies, a model based on the linear noise equation without feedback is not physically relevant. At other low frequencies ${\omega }_{\mathrm{VLF}}^{\text{cb}}<\omega \le {\omega }_{\text{LF}}$, the algorithm should be convergent, but the problem is again nearly-critical and we expect the same behavior as for the flight-based method. In the plateau region and slightly above the upper cut-off frequency, until a frequency $f_{\text{CC}}^{\text{cb}}={\omega }_{\text{CC}}^{\text{cb}/2}\pi$, we have $\tilde{\mathfrak{c}}<1$ and $\tilde{\mathfrak{K}}<1$, and the algorithm should be convergent. However, for higher frequency, the algorithm should only be conditionally convergent and weight cancellation would be necessary. Finally, at some frequency higher than $f_{D^{\text{cb}}}={\omega }_{\text{D}}^{\text{cb}/2}\pi$, we have $\tilde{\mathfrak{c}}\ge 1$. Recalling that weight cancellation inside generations is not considered here, the algorithm should therefore be divergent for these very high frequencies. Weight cancellation inside generations would in any case yield small gains.

The cut-off frequencies can be computed exactly. The cut-off frequency ${\omega }_{\text{D}}^{\text{cb}}$ is defined by the equation $\tilde{\mathfrak{c}}=1$, which yields

$$\begin{array}{c}\hfill {\omega }_{\mathrm{D}}^{\text{cb}}=1/\kappa \mathrm{\Lambda }\end{array}$$(64)

The cut-off frequencies ${\omega }_{\text{VLF}}^{\text{cb}}$ and ${\omega }_{\mathrm{CC}}^{\text{cb}}$ are defined by the equation $\tilde{\mathfrak{K}}=1$. Solving it reveals a third-degree equation whose roots read⁷

$$\begin{array}{c}\hfill {\omega }_{\mathrm{VLF}}^{\text{cb}}=\frac{2}{3}\frac{\beta }{\kappa \mathrm{\Lambda }}\times (1-\sqrt{1-3s^2}cos(\frac{1}{3}arccos\frac{1-18s^2}{{(1-3s^2)}^{\frac{3}{2}}}))\end{array}$$(65)

$$\begin{array}{c}\hfill {\omega }_{\mathrm{CC}}^{\text{cb}}=\frac{2}{3}\frac{\beta }{\kappa \mathrm{\Lambda }}\times (1+\sqrt{1-3s^2}cos(\frac{1}{3}arccos\frac{1-18s^2}{{(1-3s^2)}^{\frac{3}{2}}}+\frac{\pi }{3}))\end{array}$$(66)

where

$$\begin{array}{c}\hfill s=\kappa \frac{\lambda \mathrm{\Lambda }}{\beta }\end{array}$$(67)

More handy approximate expressions, under the usual condition λ ≪ β/Λ, read

$$\begin{array}{c}\hfill {\omega }_{\mathrm{VLF}}^{\text{cb}}\approx 2\frac{\kappa {\lambda }^2\mathrm{\Lambda }}{\beta },{\omega }_{\mathrm{CC}}^{\text{cb}}\approx \beta /\kappa \mathrm{\Lambda }\end{array}$$(68)

5.1.4. Numerical validation - settings

The previous theory has been tested numerically, with the numerical data provided in Table 1. This proceeds in two steps. First, simulations of the $\mathfrak{c},\tilde{\mathfrak{c}},\mathfrak{K},\tilde{\mathfrak{K}}$

radii were performed to confirm their analytical formulas. Second, fixed-source noise simulations were performed to check their convergence behavior compared to the one predicted by their spectral radii. All simulations were performed using the Monte Carlo mini-app MGMC [25], taking advantage of the flexibility of this code to implement all the required new features.

Simulations were performed using both an analog and a non-analog game. The analog game samples a single reaction per collision, with relevant multiplicity, only applying weight multipliers of modulus 1. The non-analog game uses implicit capture, forced fission and Russian roulette, while complex weights evolve accordingly. Russian roulette is taken as norm-based. When weight cancellation is applied, we use a simple binning procedure. However, contrary to Yamamoto’s implementation [14], all binned particles are re-used, instead of re-sampling some of the particles only, similarly to Zhang’s binning procedure for signed weights [26]. Various frequencies were tested, using either the collision-based and flight-based transport methods, with 10,000 particles.

The power iteration algorithm to compute $\tilde{\mathfrak{K}}$ is structured as the usual power iteration on the extended counterpart of equation (44), using the noise transport technique described in Section 3. The eigenvalue $\tilde{\mathfrak{K}}$ is estimated as the ratio of the total statistical weight of the population, that is the sum of the absolute values of the complex weights, between two consecutive generations.Simulations for $\tilde{\mathfrak{c}}$ are performed in a similar way, with power iteration on the extended counterpart of equation (40), which translates in the following. The fission operator ℱ does not appear in this equation and therefore no fission descendant is sampled. Any particle emerging from scattering (or copy reaction) is sent to the next generation. Consequently, the generation ends after a single flight and collision of each particle - or after the first real collision of the particle, when using delta-tracking or variants. Note that the total cross-section Σ, that appears in all operators, is the same in $\tilde{\mathfrak{K}}$ and $\tilde{\mathfrak{c}}$ simulation: it is considered to be a parameter of the problem, rather than the sum of all partial cross sections. In particular, capture does not appear in any of the operators and is not handled explicitly.

Simulations for 𝔎 or 𝔠 are performed by including the application of weight cancellation, assuming that the efficiency α is sufficient for the physical eigenvalues to be dominant. The ratio of population at each consecutive generation must therefore be computed after full application of operators ${\mathcal{A}}^{(\alpha )}\tilde{\mathcal{K}}$ and ${\mathcal{A}}^{(\alpha )}\tilde{\mathcal{T}}\tilde{\mathcal{S}}$, or, in other words, after weight cancellation. However, when using the collision-based method, no $\tilde{\mathfrak{K}}$ nor 𝔎 simulation were run for frequencies for which $\tilde{\mathfrak{c}}>1$, since these would be unable to complete a single generation. For each case, 1000 inactive and 1000 active generations were used, and the standard deviation was estimated using 20 independent replicas.

The fixed-source noise simulation starts with a uniform source of particles with weight w = 1. Simulations that are expected to diverge were handled very basically. In cases where $\tilde{\mathfrak{c}}<1$, simulations were halted at a predetermined generation n_max, and the population evolution was monitored to observe the divergence. The number n_max was chosen so that the expected population size would be 1000 times the source size, i.e. $n_{\max }=3/log\tilde{\mathfrak{K}}$, when weight cancellation was not enforced; when weight cancellation was enforced, we chose n_max = 3/log𝔎. In cases where $\tilde{\mathfrak{c}}\ge 1$, we simply monitored that the algorithm got stuck after a few generations, which usually happens at the first one.

5.1.5. Numerical validation - analysis

Results of the eigenvalues power-iteration computations are presented in Figure 3, for non-analog simulations.Analog simulations give similar results and are not presented here.

In most configurations, a very good agreement is obtained between the theoretical findings and the Monte Carlo simulations, with relative differences typically lower than 3 times the standard error of the mean. In particular, the eigenvalue computed with weight cancellation is coherent with the physical spectral radius, as predicted by the linear cancellation model presented in Section 4.5.

A statistically significant discrepancy is observed in the simulation of 𝔎 at frequency 10⁴ Hz with the collision-based method: Student’s t-test yields a p-value of 1.6 × 10⁻⁴. This is the point where weight cancellation is the most dramatic, as it is expected to reduce the eigenvalue from $\tilde{\mathfrak{K}}=7.49586$ to 𝔎 = 0.42785, according to the theoretical formulas equation (61), so that an efficiency α >  0.94 is needed to compute the physical eigenvalue. In this case, because of the finite number of particles in the simulation, cancellation is largely imperfect. To substantiate this statement, a simulation was run with 50 times more particles, to enable a greater degree of weight cancellation, which resulted in a relative difference with respect to analytical formulas of 17 ± 24 pcm, statistically compatible with 0.

Results of the fixed-source computation are also satisfactory. Overall, for all tested cases we observe a coherent behavior of Monte Carlo simulations with respect to the theoretical framework. In particular, algorithms do not converge in their conditionally convergent regime, the computation does not end, any estimate of the tally seems to diverge - unless weight cancellation is applied.

5.2. Problem 2: the rod model

We now address a richer, yet exactly solvable model, namely the rod model [27]. The viable space is a homogeneous, finite, one-dimensional rod of half-length ℓ. Boundary conditions correspond to null incoming flux.Particles are bound to move in one of the two directions Ω ∈ { − 1, +1} along the rod. Scattering is assumed to be isotropic. We consider here two energy groups, with no up-scattering. The energy groups are indexed by g = 1 for the fast and g = 2 for the thermal group. The fission spectrum χ_g is in the fast group only. This amounts to setting the conditions

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{s,2\to 1}=0,{\chi }_2=0\end{array}$$(69)

We take J families of precursors indexed by 1 ≤ j ≤ J. We introduce the delay factor

$$\begin{array}{c}\hfill \mathrm{{\rm Y}}=1-\beta +\sum _{j=1}^J\frac{{\beta }_j{\lambda }_j}{{\lambda }_j+i\omega }\end{array}$$(70)

Defining ${\phi }_{\mathrm{\Omega }}^{g(x)}$ the flux at position x, direction Ω and group g, the operator $\widehat{{\mathcal{B}}_0}$ reads

$$\begin{array}{cc}\hfill \widehat{{\mathcal{B}}_0}{\phi }_{\mathrm{\Omega }}^g& =(\frac{i\omega }{v_g}+\mathrm{\Omega }\frac{d}{\mathit{dx}}+{\mathrm{\Sigma }}_{t,g}){\phi }_{\mathrm{\Omega }}^g\hfill \\ \hfill & -\frac{1}{2}\sum _{{\mathrm{\Omega }}^{\prime },g^{\prime }}{\mathrm{\Sigma }}_{s,g^{\prime }\to g}{\phi }_{{\mathrm{\Omega }}^{\prime }}^{g^{\prime }}-\frac{\mathrm{{\rm Y}}}{k_{\text{eff}}}\frac{{\chi }_g}{2}\sum _{{\mathrm{\Omega }}^{\prime },g^{\prime }}\nu {\mathrm{\Sigma }}_{f,g^{\prime }}{\phi }_{{\mathrm{\Omega }}^{\prime }}^{g^{\prime }}\hfill \end{array}$$(71)

The boundary conditions impose

$$\begin{array}{c}\hfill \forall g,\forall \mathrm{\Omega },{\phi }_{\mathrm{\Omega }}^g(-\mathrm{\Omega }\ell )=0\end{array}$$(72)

Solving this problem for the $c,\tilde{c},k$ and $\tilde{k}$ eigenvalues formulation is a straightforward but tedious task. Intermediate steps are omitted for the sake of conciseness.

5.2.1. Solutions for the critical problem

Let us start by considering the critical problem, without noise, to understand the structure of the solutions. Since there is no up-scattering, the c-eigenvalue problem is triangular: the eigenvalues are those of the two single-speed problems, for the fast and thermal group. Their coupling through Σ_{s, 1 → 2} plays no role. Using the notation of the corresponding infinite medium eigenvalues for each group, i.e.

$$\begin{array}{c}\hfill c_{\infty ,g}=\frac{{\mathrm{\Sigma }}_{s,g\to g}}{{\mathrm{\Sigma }}_{t,g}}\end{array}$$(73)

we obtain

$$\begin{array}{c}\hfill c_{\text{eff}}=\underset{1\le g\le 2}{max}\frac{c_{\infty ,g}}{1+{cot}^2{\gamma }_g^0}\end{array}$$(74)

where γ_g⁰ is a geometric factor, corresponding to the smallest positive solution of

$$\begin{array}{c}\hfill \gamma tan\gamma =\ell {\mathrm{\Sigma }}_{t,g}\end{array}$$(75)

Anticipating the results given in Section 5.2.2, we will have to consider geometric factors γ that are solutions of the equation

$$\begin{array}{c}\hfill \gamma (tan\gamma +ucot\gamma )=a(1-u)\end{array}$$(76)

In this equation, a is a dimensionless scale factor, and u is an anisotropy parameter. Then, the geometric factor γ_g⁰ is the smallest positive solution for parameters

$$\begin{array}{cc}\hfill {\gamma }_{g^0}:& \{\begin{array}{c}a=\ell {\mathrm{\Sigma }}_{t,g}\hfill \\ u=0\hfill \end{array}\hfill \end{array}$$(77)

In order to interpret the results for the k-eigenvalue problem, we rely on a variant of the six-factor formula from reactor physics, which reads

$$\begin{array}{c}\hfill k_{\text{eff}}=ϵp\xi P_1{P}_2\end{array}$$(78)

where ξ is the product of the thermal fission factor and thermal utilization factor, p is the resonance escape probability, ϵ is the fast-fission factor, P₁ and P₂ are the fast and thermal non-leakage probability, respectively. Their expressions are

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \xi & =\frac{\nu {\mathrm{\Sigma }}_{f,2}}{({\mathrm{\Sigma }}_{t,2}-{\mathrm{\Sigma }}_{s,2\to 2})},\hfill \\ \hfill p& =\frac{{\mathrm{\Sigma }}_{s,1\to 2}}{({\mathrm{\Sigma }}_{t,1}-{\mathrm{\Sigma }}_{s,1\to 1})},\hfill \\ \hfill ϵ& =1+\frac{\nu {\mathrm{\Sigma }}_{f,1}}{({\mathrm{\Sigma }}_{t,1}-{\mathrm{\Sigma }}_{s,1\to 1})}\frac{1}{P_{2}p\xi },\hfill \\ \hfill P_g& =\frac{1}{1+{\alpha }_0^{{2M}_{g^2}}/{\ell }^2}\hfill \end{array}\end{array}$$(79)

where we introduced the migration area in each group,

$$\begin{array}{c}\hfill M_g^2=\frac{1}{{\mathrm{\Sigma }}_{t,g}({\mathrm{\Sigma }}_{t,g}-{\mathrm{\Sigma }}_{s,g\to g})}\end{array}$$(80)

The geometric parameter α₀ occurring in equation (79) is the smallest positive solution of the following system of equations for (α, 𝜚):

$$\begin{array}{cc}& \{\begin{array}{c}{\varrho }^2-{\alpha }^2={\ell }^2(\frac{1}{M_1^2}+\frac{1}{M_2^2}-\frac{\nu {\mathrm{\Sigma }}_{f,1}{\mathrm{\Sigma }}_{t,1}}{k_{\text{eff}}})\hfill \\ \frac{1-\frac{\alpha }{\ell {\mathrm{\Sigma }}_{t,1}}tan\alpha }{1+\frac{\varrho }{\ell {\mathrm{\Sigma }}_{t,1}}tanh\varrho }=\frac{{\ell }^2-M_2^2{\varrho }^2}{{\ell }^2+M_2^2{\alpha }^2}\frac{1-\frac{\alpha }{\ell {\mathrm{\Sigma }}_{t,2}}tan\alpha }{1+\frac{\varrho }{\ell {\mathrm{\Sigma }}_{t,2}}tanh\varrho }\hfill \end{array}\hfill \end{array}$$(81)

The k_eff eigenvalue can be replaced by a function of α₀ using the previous expressions. The parameter 𝜚₀ is conjugate to α₀, in that the same numerical value of k_eff is obtained by swapping α₀ with −𝜚₀ in the expressions of P_g. In fact, ℓ/α₀ is the characteristic length of the sinusoidal part of the solution, while ℓ/𝜚₀ is the characteristic length of the hyperbolic part of the solution, whence the tangent and hyperbolic tangent functions in the previous expressions.

More generally, we will have to consider solutions α of the following system of equations on (α, 𝜚):

$$\begin{array}{cc}& \{\begin{array}{c}{\varrho }^2-{\alpha }^2=m_1+m_2-g(\alpha ),\hfill \\ \frac{1-\frac{\alpha }{s_1}tan\alpha }{1+\frac{\varrho }{s_1}tanh\varrho }=\frac{m_2-{\varrho }^2}{m_2+{\alpha }^2}\frac{1-\frac{\alpha }{s_2}tan\alpha }{1+\frac{\varrho }{s_2}tanh\varrho }\hfill \end{array}\hfill \end{array}$$(82)

This equation depends on dimensionless parameters m₁, m₂, s₁, s₂ and function g. Then, the geometric factor α₀ is the smallest positive solution of equation (82) for parameters

$$\begin{array}{cc}\hfill {\alpha }_0:& \{\begin{array}{c}{m}_g={\ell }^2/M_g^2\hfill \\ s_g=\ell {\mathrm{\Sigma }}_{t,g}\hfill \\ g(\alpha )={\ell }^2\nu {\mathrm{\Sigma }}_{f,1}{\mathrm{\Sigma }}_{t,1}/k_{\text{eff}}\hfill \end{array}\hfill \end{array}$$(83)

When referring to solutions of equations (76) or (82), we implicitly refer to the solution that yields the smallest spectral radius. It is usually, yet not systematically, the solution with the smallest real part.

5.2.2. Noise eigenvalues

Concerning the c-eigenvalue problem, the spectral radii are the maximum of the spectral radius of each group. For the flight-based method, we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathfrak{c}}_2^{\text{fb}}& =\underset{1\le g\le 2}{max}\left|\frac{c_{\infty ,g}}{(1+i\omega {\tau }_g)(1+{cot}^2{\gamma }_{g^{\text{fb}}})}\right|,\hfill \\ \hfill {\tilde{\mathfrak{c}}}_{2^{\text{fb}}}& =c_{\text{eff}}\hfill \end{array}\end{array}$$(84)

where ${\gamma }_g^{\text{fb}}$ is solution of equation (76) for parameters

$$\begin{array}{cc}\hfill {\gamma }_g^{\text{fb}}:& \{\begin{array}{c}a=\ell {\mathrm{\Sigma }}_{t,g}(1+i\omega {\tau }_g),\hfill \\ u=0\hfill \end{array}\hfill \end{array}$$(85)

We have introduced the characteristic times

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_g=\frac{1}{v_g({\mathrm{\Sigma }}_{t,g}-{\mathrm{\Sigma }}_{s,g\to g})},{\tau }_g=\frac{1}{v_g{\mathrm{\Sigma }}_{t,g}}\end{array}$$(86)

which can be interpreted as the mean neutron lifetime and the time to collision in each group, neglecting leakage.

For the collision-based method, we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathfrak{c}}_2^{\text{cb}}& =\underset{1\le g\le 2}{max}|\frac{c_{\infty ,g}+(\eta -i)\omega {\tau }_g}{1+\eta \omega {\tau }_g}+\frac{{\gamma }_{g^{\text{cb}}}cot{\gamma }_{g^{\text{cb}}}}{\ell {\mathrm{\Sigma }}_{t,g}}|,\hfill \\ \hfill {\tilde{\mathfrak{c}}}_2^{\text{cb}}& =\underset{1\le g\le 2}{max}\frac{c_{\infty ,g}+\sqrt{1+{\eta }^2}\omega {\tau }_g}{1+\eta \omega {\tau }_g+\frac{\stackrel{~}{{\gamma }_{g^{\text{cb}}}}cot\widetilde{{\gamma }_g^{\text{cb}}}}{\ell {\mathrm{\Sigma }}_{t,g}}}\hfill \end{array}\end{array}$$(87)

where the geometric factors are solutions of solution of equation (76) respectively for parameters

$$\begin{array}{cc}\hfill {\gamma }_g^{\text{cb}}:& \{\begin{array}{c}a=\ell {\mathrm{\Sigma }}_{t,g}(1+\eta \omega {\tau }_g),\hfill \\ u=\frac{(\eta -i)\omega {\tau }_g}{(\eta -i)\omega {\tau }_g+c_{\infty ,g}}\hfill \end{array}\hfill \end{array}$$(88)

and

$$\begin{array}{cc}\hfill {\tilde{\gamma }}_{g^{\text{cb}}}:& \{\begin{array}{c}a=\ell {\mathrm{\Sigma }}_{t,g}(1+\eta \omega {\tau }_g)\hfill \\ u=\frac{\sqrt{1+{\eta }^2}\omega {\tau }_g}{\sqrt{1+{\eta }^2}\omega {\tau }_g+c_{\infty ,g}}\hfill \end{array}\hfill \end{array}$$(89)

For the 𝔎 and $\tilde{\mathfrak{K}}$ eigenvalues problems, we can express the results using the six-factor formula. Regardless of the method, the physical radius is

$$\begin{array}{c}\hfill {\mathfrak{K}}_2={\mathfrak{K}}_2^{\text{fb}}={\mathfrak{K}}_2^{\text{cb}}=\left|k_2\right|\end{array}$$(90)

with k₂ the complex eigenvalue being

$$\begin{array}{c}\hfill k_2=\mathrm{{\rm Y}}\frac{P_1^{\prime }}{P_1}\frac{1}{ϵ}(ϵ-1+\frac{P_2^{\prime }}{P_2})\end{array}$$(91)

where P_g′ takes into account leakage and dampening, i.e.

$$\begin{array}{c}\hfill P_g^{\prime }=\frac{1}{1+i\omega {\mathrm{\Lambda }}_g+\frac{{({\alpha }^{\prime })}^2}{1+i\omega {\tau }_g}\frac{M_g^2}{{\ell }^2}}\end{array}$$(92)

and we have introduced the extended delay factor

$$\begin{array}{c}\hfill \widetilde{\mathrm{{\rm Y}}}=1-\beta +\sum _{j=1}^J\frac{{\beta }_j{\lambda }_j}{|{\lambda }_{j+i}\omega |}\end{array}$$(93)

In equation (92), the geometric factor α′ is obtained as a solution of equation (82) for parameters

$$\begin{array}{cc}\hfill {\alpha }^{\prime }:& \{\begin{array}{c}{m}_g={\ell }^2(1+i\omega {\tau }_g)(1+i\omega {\mathrm{\Lambda }}_g)/M_g^2\hfill \\ s_g=\ell {\mathrm{\Sigma }}_{t,g}(1+i\omega {\tau }_g)\hfill \\ g(\alpha )={\ell }^2\nu {\mathrm{\Sigma }}_{f,1}{\mathrm{\Sigma }}_{t,g}(1+i\omega {\tau }_1)/k_2\hfill \end{array}\hfill \end{array}$$(94)

The extended spectral radii are

$$\begin{array}{c}\hfill {\tilde{\mathfrak{K}}}_2^{\text{fb}}=\widetilde{\mathrm{{\rm Y}}},{\tilde{\mathfrak{K}}}_2^{\text{cb}}=\widetilde{\mathrm{{\rm Y}}}\frac{{\tilde{P}}_1^{\text{cb}}}{P_1}\frac{1}{ϵ}(ϵ-1+\frac{{\tilde{P}}_2^{\text{cb}}}{P_2})\end{array}$$(95)

where ${\tilde{P}}_{g^{\text{cb}}}$ takes into account leakage and dampening, i.e.

$$\begin{array}{c}\hfill {\tilde{P}}_g^{\text{cb}}=\frac{1}{1-\kappa \omega {\mathrm{\Lambda }}_g+\frac{{\tilde{\alpha }}^2}{1-\kappa \omega {\tau }_g}\frac{M_g^2}{{\ell }^2}}\end{array}$$(96)

The geometric factor $\tilde{\alpha }$ is obtained as a solution of equation (82) for parameters

$$\begin{array}{cc}\hfill \tilde{\alpha }:& \{\begin{array}{c}{m}_g={\ell }^2(1-\kappa \omega {\tau }_g)(1-\kappa \omega {\mathrm{\Lambda }}_g)/M_g^2\hfill \\ s_g=\ell {\mathrm{\Sigma }}_{t,g}(1-\kappa \omega {\tau }_g)\hfill \\ g(\alpha )={\ell }^2\nu {\mathrm{\Sigma }}_{f,1}{\mathrm{\Sigma }}_{t,g}(1-\kappa \omega {\tau }_1)/{\tilde{\mathfrak{K}}}_2^{\text{cb}}\hfill \end{array}\hfill \end{array}$$(97)

5.2.3. Analysis of the eigenvalues

The theoretical values given in equations (84), (87), (90), and (95) are displayed in Figure 4, using the numerical values given in Table 2, which were chosen to be comparable with those of Problem 1. In particular, the cross sections of the thermal group are identical, except for the fission cross-section, which in Problem 1 compounds the resonance escape probability. The one-precursor family for Problem 1 roughly corresponds to the one-family approximation of the six families from Problem 2, with comparable total delayed fraction and effective decay constants. The half-length ℓ of the rod was chosen to ensure criticality.

The dependence of the geometrical factors γ and α with respect to the angular frequency ω is weak and has no significant impact.

The two-group structure impacts the regime transition from plateau to high frequencies, in particular the value of ω_CC^cb. It is of small consequences in this example, but we may not exclude that a multiple group or continuous energy structure may affect more significantly this transition.

Having multiple precursors families is responsible for the behavior of the 𝔎 and $\tilde{\mathfrak{K}}$ spectral radii around the low-frequencies cut-off f_LF. Their decrease is spread from λ₁/2π to λ₆/2π, whose values lie from the low frequency cut-off to the mid-plateau frequencies, instead of being centered at the low frequency cut-off f_LF.

Nevertheless, the behavior of Problem 2 is qualitatively similar to the one of Problem 1. Both models have the same convergence regime, with only numerical corrections to the convergence speed and cut-off between regimes brought by the more complex structure of the equations associated to Problem 2.

5.2.4. Numerical validation

The previous theory has been tested numerically, as described for Problem 1 in Section 5.1.4, with the exception of 𝔠 simulations that were not performed. Contrary to Problem 1, where particle directions can be safely neglected for the purpose of weight cancellation, in Problem 2 (as well as in more general settings) the estimation of 𝔠 would have required a special implementation of weight cancellation taking into account directions and energy groups, which we did not investigate.

Results of the eigenvalues simulations for non-analog simulations are presented in Figure 5, and compared to the analytical formulas given in Section 5.2.2, used as a reference. Analog simulations give similar results and are not presented here.

Results are similar to those of Problem 1. Statistical agreement is found between theoretical findings and Monte Carlo simulations, for both the power iteration eigenvalue computation and the fixed-source convergence computation. A discrepancy is observed in the simulation of $\tilde{\mathfrak{c}}$ with the flight-based method; the amplitude of the discrepancy is however rather small, of the order of less than 1 pcm. The origin of this disagreement is currently under investigation, and is thought to be due to truncation errors.

6. Discussion

The theoretical framework that has been applied in this work to complex-valued Monte Carlo for neutron noise problems is in fact more general, and can be in principle applied to any Monte Carlo games involving non-positive-weighted particles. In reactor physics applications, this would include e.g. the class of heterogeneous buckling models [28], or negative-weighted delta-tracking methods for multi-physics problems involving spatially continuous cross sections [29, 30]. The same strategy might also be extended to complex-valued equations beyond the domain of reactor physics, e.g. for fermion wave-function computation [31].

The two proposed benchmark models examined in Section 5, despite admittedly drastic simplifications, do provide hints and general trends that are expected to apply to less restrictive conditions. These models may provide useful estimates of the extended radii $\tilde{\mathfrak{K}}$ and $\tilde{\mathfrak{c}}$, as well as the cut-off frequencies, provided effective parameters have been computed. Their accuracy in realistic settings has not been probed. Potential difficulties may reside in the influence of the energy structure on the spectral radii at high frequencies. Besides, the estimation of the radius with weight cancellation may suffer from imperfect cancellation, the efficiency of weight cancellation being generally unknown.

Overall, the models reproduce the convergence pattern observed in the literature for other configurations and benchmarks [16]. In particular, they show that in the plateau region convergence is achieved for both collision-based and flight based methods without need of weight cancellation. For the collision-based method, convergence without weight cancellation carries on to higher frequencies, up to a previously unidentified limit that we were able to identify as the situation where $\tilde{\mathfrak{K}}=1$, which yields the cut-off frequency ω_CC^cb ≈ β/κΛ.

Nevertheless, some observations are new and contrast with Rouchon’s findings reported in Ref. [16], to which we dedicate the following discussion.

6.1. The effect of Russian roulette

Reference [16] reported poorer convergence properties without weight cancellation than what we observed in our investigations. In particular, it was stated that weight cancellation was necessary for convergence with the flight-based method, outside the plateau region. This discrepancy is attributed to the use of two-part Russian roulette, as opposed to the norm-based proposed in this work. In fact, although the two-part Russian roulette provides unbiased results, it displaces particles in the θ dimension. It can be shown that this process does indeed increase the spectral radii compared to simulations without Russian roulette, as discussed in Appendix A. To support this hypothesis, computation of $\tilde{\mathfrak{K}}$ and $\tilde{\mathfrak{c}}$ eigenvalues with the two-part roulette were performed, as illustrated in Figure 6, along with fixed source computations. We observe that the spectral radii with the two-part roulette are consistently higher than with the norm-based roulette. On the opposite, the norm-based roulette can be shown not to increase the spectral radius compared to simulations without roulette. This fact is supported by our observation of similar radii in the analog and non-analog games tested.

Consequently, the cut-off frequencies ω_CC^cb and ω_D^cb are lowered when using the two-part Russian roulette, compared to the norm-based roulette or no roulette. Therefore, we suggest to use the new implementation, since it displays better convergence properties than the two-part roulette at no significant additional computation cost.

6.2. Analog games

Previous works claimed that, using the collision-based method without weight cancellation, convergence was possible at higher frequencies if a pseudo-analog game, rather than a non-analog game, was performed [16]. Our theoretical framework, even using a more restrictive definition of analog games, does not support this claim, and our Monte Carlo simulations using either form of games did not show any substantial difference.

We believe that the observations reported in Ref. [16] are in fact also due to the two-part implementation of the Russian roulette, since the pseudo-analog game did not use any roulette, contrary to the non-analog game. Switching to the pseudo-analog game allow to get rid of the weakened convergence properties of the two-part Russian roulette, and in particular sets ${\omega }_{\mathrm{CC}}^{\text{cb}}$ at higher frequencies. The symmetric divergence of the positive and negative total weights previously observed with the non-analog game, and therefore with the two-part Russian roulette, is typical of the Conditionally Convergent regime.

We also issue a call for caution. The implementation of pseudo-analog games in Ref. [16], contrary to our current analog games, used in the copy reaction a weight multiplier (η − i)/η, whose modulus is greater than 1,without enforcing any splitting procedure. In the Conditionally Convergent regime, any particle dies almost surely yet the particle expected weight diverge with respect to the number of generation. This game is feasible, in the sense of Lux and Koblinger of all particle histories ending with probability 1 [12], yet not convergent. In this case, the expected value of the tally estimator seems undefined. Such game may be misleading : it provides a value for the tally estimator anytime it is run yet this value does not allow to estimate the desired tally. Only independent replica may hint at this statistical discrepancy, which are not systematically run by Monte Carlo users.

6.3. Adjustment of the η factor

Rouchon’s works have shown that, when using the collision-based method, increasing the η factor may overcome convergence issues at high frequencies [16]. This statement is indeed reproduced in the derivation discussed in equation (68): the cutoff frequencies ${\omega }_{\mathrm{CC}}^{\text{fb}}$ and ${\omega }_{\mathrm{D}}^{\text{fb}}$ are increasing functions of η.

However, a severe limitation does appear, as shown using Problem 1 for illustration. For any angular frequency, there exists a unique optimal choice η′(ω) that yields the lowest $\tilde{\mathfrak{c}}$ value,

$$\begin{array}{c}\hfill {\eta }^{\prime }(\omega )=\frac{\omega \tau +c_{\infty }\sqrt{1-c_{\infty }^2+{\omega }^2{\tau }^2}}{1-c_{\infty }^2}\end{array}$$(98)

The resulting intra-generation ‘reactivity’, $(\tilde{\mathfrak{c}}-1)/\tilde{\mathfrak{c}}$ is shown in Figure 7 for η = 1 and η = η′(ω). In particular, the optimal reactivity is negative, but tends rapidly to 0 at high frequencies. Generations can be simulated, yet their simulation time explodes dramatically. Therefore, adjusting the η factor is not efficient for arbitrarily high frequencies, for which the flight-based method is recommended.

7. Conclusions

In this work, we have investigated the convergence properties of Monte Carlo games for the linear neutron noise equation. We have proposed a general framework to characterize the convergence of a broad class of Monte Carlo games, including complex-weighted Monte Carlo games and generation-decomposition strategies, that allows for a thorough description of neutron noise algorithms using flight-based and collision-based methods. This framework was successfully confirmed by numerical experiments on analytically solvable toy models.

Relying on the theory developed in our work, we were able to complete Rouchon’s description of the convergence regimes of the neutron noise algorithms. Besides, we gained valuable insight that allowed formulating practical recommendations and improvements. On one hand, we showed that the novel norm-based Russian roulette is to be preferred systematically. In particular, it significantly extends the applicability of the flight-based method to all frequencies, without need of weight cancellation. On the other hand, we provided formulas for the cut-off frequencies of the collision-based algorithm in a point kinetic approximation in equation (68), which can be easily computed beforehand, during the power iteration part of the noise computation. These findings point at what value of η parameter to choose for the collision-based method, or even whether a flight-based method might be preferred altogether.

The contribution of weight cancellation has been addressed. Although the framework developed in this work does not allow for a formal treatment of this process, the hints gathered in our investigations suggest under which conditions weight cancellation becomes mandatory. Remark, however, that weight cancellation has been more broadly proven efficient for variance reduction, with a large impact on the figure of merit. Therefore, our general recommendation is to enforce weight cancellation whenever possible. In some circumstances, weight cancellation might turn out to be inefficient, e.g. in higher dimensions or if the generation-decomposition involves some non-isotropic reactions in addition to (isotropic) fission; in such cases, a proper knowledge of the extended eigenvalues would be of interest. Besides, weight cancellation has been observed to reduce variance, but this effect is yet to be quantified theoretically.

Moreover, although we exhibit eigenvalues, being able to estimate the figure of merit would be more valuable in view of selecting an algorithm, for instance choosing between the flight-based and collision-based method. A figure of merit requires the estimate of the expected computing time and of the variance. These quantities depend on the sampling strategy used in a given game, which has to be described in a broader formalism encompassing not only the location of the particle in the phase space and the weight, but also the multiplicity in branching processes and potential correlations, which is beyond the scope of this paper. Besides, even with a complete representation of a game, the second moment equations are generally hard to exploit [12].

Future work will focus on probing the present framework to real-world applications. Notably, Monte Carlo algorithms for neutron noise have been seldom applied to continuous-energy problems, and in particular never outside the plateau region, to our best knowledge. Besides, the question of computation time and variance reduction could be addressed. Finally, one should also question the physical relevance of the equation in the first place: the orthodox linearization is known to fail in the case of mechanical vibrations [32]. Monte Carlo algorithms capable of going beyond the linearized equations, for which preliminary work has recently been proposed [33], would surely be of utmost interest.

---

Funding

This research received no external funding.

Conflicts of interest

The authors have nothing to disclose.

Data availability statement

Data associated to this paper can be obtained from the corresponding author, upon reasonable request.

Author contribution statement

Conceptualization, A.R., D.M., and A.Z.; Methodology, A.F., A.R., D.M., and A.Z.; Formal Analysis, A.F.; Software, A.F.; Validation, A.F.; Investigation, A.F., A.R., D.M., and A.Z.; Writing – Original Draft Preparation, A.F.; Writing – Review & Editing, A.F., A.R., D.M., and A.Z.

References

Appendix

Alternatives choices of kernels

In Section 4.1, it is noted that the choice of kernels Z_μ and extended source $\tilde{q}$ is in fact not unique. Other sets of kernels and functions may be used, provided that they faithfully represent multiplication by a complex unitary factor e^iμ , and source in the extended phase space. This is formally written as

$$\begin{array}{c}\hfill \forall {\theta }^{\prime }\in {\mathcal{S}}_1,{\int }_{{\mathcal{S}}_1}{e}^{i\theta }{Z}_{\mu (\theta ,{\theta }^{\prime })}d\theta =e^{i(\mu +{\theta }^{\prime })}\end{array}$$(A.1)

$$\begin{array}{c}\hfill \forall x\in \mathcal{D},{\int }_{{\mathcal{S}}_1}{e}^{i\theta }\tilde{q}(x,\theta )d\theta =q(x)\end{array}$$(A.2)

For illustration, in a previous work we exhibited a four-species system, for the four cardinal directions of the complex plane {1, i, −1, −i}, which corresponds to restricting θ to multiples of π/2 [19]. This system requires, upon multiplication by a complex number, to split each particle in two, one for the real part and one for the imaginary part. Writing for readability e^iμ  = a + ib, the corresponding kernels are in matrix form ⁸

$$\begin{array}{c}\hfill Z_{\mu }=\left(\begin{array}{cccc}{a}_{+}& b_{-}& a_{-}& b_{+}\\ b_{+}& a_{+}& b_{-}& a_{-}\\ a_{-}& b_{+}& a_{+}& b_{-}\\ b_{-}& a_{-}& b_{+}& a_{+}\end{array}\right).\end{array}$$(A.3)

Here, subscripts + and − denote the positive and negative parts, with

$$\begin{array}{c}\hfill \forall x\in \mathbb{R},x_{+}=max(x,0),x_{-}=max(-x,0).\end{array}$$(A.4)

The quantity Z_μ is a Toeplitz, circulant matrix, whose eigenvalues are |a| + |b|, a + ib, a − ib and |a| − |b|.

However, this system consistently showed poorer convergence properties than the usual complex-weight system. This can be readily understood: the spectral radius of this matrix is |a| + |b| ≥ |a+ib| = 1, so that on average the expected total statistical weight of the children particles is larger than the original particle weight. In other words, because of the splitting into real and imaginary parts of the particle, the multiplication by a unitary complex number becomes a multiplicative process.

This finding is of general scope. Let us denote by Z_μ the kernels defined in equation (26), and let Z_μ′ be another set of kernels satisfying equation (99), associated to the extended operators $\tilde{\mathcal{K}}$ and $\widetilde{\mathcal{K}\prime }$ respectively. We can show that the spectral radius $\tilde{\mathfrak{K}}\prime$ of $\widetilde{\mathcal{K}\prime }$ is larger or equal to the one $\tilde{\mathfrak{K}}$ of $\tilde{\mathcal{K}}$, in the functional space $L^1(\tilde{\mathcal{D}})$. Observe that

$$\begin{array}{c}\hfill \forall \theta \in {\mathcal{S}}_1,{\int }_{{\mathcal{S}}_1}{Z}_{\mu }(\theta ,{\theta }^{\prime })d\theta =1,\end{array}$$(A.5)

while, according to the triangular inequality, we have

$$\begin{array}{c}\hfill \forall \theta \in {\mathcal{S}}_1,{\int }_{{\mathcal{S}}_1}{Z}_{\mu }^{\prime }(\theta ,{\theta }^{\prime })d\theta \ge 1.\end{array}$$(A.6)

Let us introduce, for any $\tilde{\psi }\in L^1(\tilde{\mathcal{D}})$, the function $[\tilde{\psi }]\in L^1(\mathcal{D})$, with

$$\begin{array}{c}\hfill [\tilde{\psi }]:x\mapsto {\int }_{{\mathcal{S}}_1}\left|\tilde{\psi }(x,\theta )\right|d\theta .\end{array}$$(A.7)

Then, it is straightforward to show that,

$$\begin{array}{c}\hfill \forall \tilde{\psi }\in L^1(\tilde{\mathcal{D}}),[\tilde{\mathcal{K}}\tilde{\psi }]\le \mathcal{K}[\tilde{\psi }]\le [{\tilde{\mathcal{K}}}^{\prime }|\tilde{\psi }|].\end{array}$$(A.8)

By immediate induction, for any such $\tilde{\psi }$ and n ∈ ℕ, we have

$$\begin{array}{c}\hfill [{\tilde{\mathcal{K}}}^n\tilde{\psi }]\le {\mathcal{K}}^n[\tilde{\psi }]\le [{({\tilde{\mathcal{K}}}^{\prime })}^n|\tilde{\psi }|].\end{array}$$

Therefore, integrating over 𝒟 and dividing by the norm of $\tilde{\psi }$, we obtain

$$\begin{array}{c}\hfill \frac{||{\mathcal{K}}^n\tilde{\psi }||}{||\tilde{\psi }||}\le \frac{||{({\tilde{\mathcal{K}}}^{\prime })}^n|\tilde{\psi }|||}{||\tilde{\psi }||}\le |||{({\tilde{\mathcal{K}}}^{\prime })}^n|||.\end{array}$$

Taking the supremum on all $\tilde{\psi }\in L^1(\tilde{\mathcal{D}})$, we have

$$\begin{array}{c}\hfill |||{\tilde{\mathcal{K}}}^n|||\le |||{({\tilde{\mathcal{K}}}^{\prime })}^n|||.\end{array}$$

Using Gelfand’s formula, we finally obtain

$$\begin{array}{c}\hfill \tilde{\mathfrak{K}}\le {\tilde{\mathfrak{K}}}^{\prime }.\end{array}$$(A.9)

The choice of the Z_μ kernel from equation (26) leads to the smallest possible spectral radius, in $L^1(\tilde{\mathcal{D}})$, and therefore to the fastest convergence. We can also infer that it minimizes the variance of the tallies, although no demonstration is provided.

Appendix

Statistical weights

The extension procedure and the discussion on convergence in Section 4.1 and 4.6 can be generalized to account for statistical weights, i.e. the usual, real positive weights, as well as the Russian roulette. Let 𝒲 ⊂ ℝ₊ and consider the kernel K to be a positive-weighted sum of non-negative kernels, namely

$$\begin{array}{c}\hfill K(x,x^{\prime })=\sum _j{a}_j(x,x^{\prime })K_j(x,x^{\prime }).\end{array}$$(B.1)

Introducing a set of kernels W_a(w, w′) on the dimension 𝒲, the same procedure as above can be repeated mutatis mutandis. In particular, any choice of kernels that satisfies

$$\begin{array}{c}\hfill \forall w^{\prime }\in \mathcal{W},{\int }_{\mathcal{W}w}{W}_a(w,w^{\prime })dw=a{w}^{\prime }\end{array}$$(B.2)

is admissible. The extended equation only involves non-negative quantities, and is thus suitable for Monte Carlo.

The subtlety is that one should then use the measure wdwdx to define norms on $L^1(\tilde{\mathcal{D}})$. This ensures that $\tilde{\mathcal{K}}$ having a spectral radius strictly smaller than 1 still warrants convergence of the tally, and that the equality $\tilde{R}=R$ of the tallies still holds. Moreover, one can then show that the physical and extended operators have the same spectral radii, i.e.

$$\begin{array}{c}\hfill \mathfrak{K}=\tilde{\mathfrak{K}},\mathfrak{c}=\tilde{\mathfrak{c}}.\end{array}$$(B.3)

In other words, substituting positive numbers by kernels on positive weights does not impact convergence.

B.1. Standard Russian roulette

A standard Russian Roulette converts statistical weights into probability of existence. It does not move particles in the other dimensions of the phase space. It can be encoded in the same way as a $W_x^{\mathit{RR}}(w,w\prime )$ kernel, acting on the w dimension only but possibly depending on the other variables, for instance via a position-dependent threshold. This kernel must verify

$$\begin{array}{c}\hfill \forall x\in \mathcal{D},\forall w^{\prime }\in \mathcal{W},{\int }_{\mathcal{W}}w{W}_x^{\mathit{RR}}(w,w^{\prime })dw=w^{\prime }.\end{array}$$(B.4)

Similarly as substituting positives numbers, operators with and without the Russian roulette have the same spectral radius, and consequently the same convergence properties. However, Russian roulette does impact the feasibility of the game, which Lux and Koblinger define as the property by which the Monte Carlo game ends with probability 1 [12].

Multiplicative processes and splitting procedures can also be taken into account. Again, a standard splitting procedure converts statistical weights into a number of identical and independent particles, and should not move the particles in the other dimensions of the phase space. One should then consider operators on the particle density, i.e. the sum of each particle probability density: then it will be apparent that such a procedure does not impact the spectral radius of the operators. A full description is beyond the scope of this discussion.

B.2. Non-analog complex-weighted games

Non-analog complex weighted games, with weights w ∈ ℂ, can be treated by considering the complex weight modulus as a statistical weight. The extended equation then has two additional dimensions: the phase θ = argw and the statistical weight $\tilde{w}=\left|w\right|$.

An immediate consequence is that the norm-based roulette becomes then a standard roulette. Therefore, using non-unitary complex weights, non-analog sampling strategies or the norm-based roulette will not affect the spectral radii: analysis of convergence properties can focus on the θ dimension only.

On the contrary, the two-part Russian roulette is not a standard roulette: it acts simultaneously on the two additional dimensions, moving particles around in the θ dimension. It can be modeled by an integral operator acting on both $\tilde{w}$ and θ. Remembering that the statistical weight is the norm of the complex weight, elementary calculations show that the expected value of the statistical weight after the roulette is larger or equal than the statistical weight before the roulette. Therefore, the two-part roulette acts as a multiplicative process: although on average it preserves the total complex weight, it increases the total statistical weight, and consequently increases the spectral radius of the operator. We can also conjecture that it might increase the variance of the tally.

Appendix

Expressions of noise operators

For completeness, we provide the expressions of the different physical operators 𝒯, 𝒮 and ℱ that form $\widehat{{\mathcal{B}}_0}$ in the linear noise theory. Their extended counterparts are derived as explained before. The flight-based or collision-based method correspond to assigning the iω/v term to either the operator ℒ or 𝒮 in equation (35), and to the corresponding choices of Σ. Using the flight-based method, the unknown is ${\psi }_{\text{fb}}=⟨{\mathrm{\Sigma }}_t⟩\widehat{\delta \phi }$, and the operators are

$$\begin{array}{c}\hfill {\mathcal{T}}_{\text{fb}}=\int ⟨{\mathrm{\Sigma }}_t⟩(\overrightarrow{r})e^{-i\frac{\omega }{v}s}{e}^{-{\int }_0^s⟨{\mathrm{\Sigma }}_t⟩(\overrightarrow{r}-s^{\prime }\overrightarrow{\mathrm{\Omega }})d{s}^{\prime }}ds,\end{array}$$(C.1)

$$\begin{array}{c}\hfill {\mathcal{S}}_{\text{fb}}=\int \int \frac{⟨{\mathrm{\Sigma }}_{s{f}_s}⟩}{⟨{\mathrm{\Sigma }}_t⟩}d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime },\end{array}$$(C.2)

and

$$\begin{array}{cc}\hfill {\mathcal{F}}_{\text{fb}}& =\frac{1}{4\pi k_{\text{eff}}}\int \int \frac{⟨{\chi }_p{\nu }_p{\mathrm{\Sigma }}_f⟩}{⟨{\mathrm{\Sigma }}_t⟩}d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }\hfill \\ \hfill & +\frac{1}{4\pi k_{\text{eff}}}\sum _{j=1}^J\frac{{\lambda }_j}{{\lambda }_j+i\omega }\int \int \frac{⟨{\chi }_{d,j}{\nu }_{d,j}{\mathrm{\Sigma }}_f⟩}{⟨{\mathrm{\Sigma }}_t⟩}d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }.\hfill \end{array}$$(C.3)

Conversely, using the collision-based method, the unknown is ${\psi }_{\text{cb}}={\tilde{\mathrm{\Sigma }}}_t^{\text{cb}}\widehat{\delta \phi }$, and the operators are

$$\begin{array}{c}\hfill {\mathcal{T}}_{\text{cb}}=\int {\tilde{\mathrm{\Sigma }}}_{t^{\text{cb}(\overrightarrow{r})}}{e}^{-{\int }_0^s{\tilde{\mathrm{\Sigma }}}_{t^{\text{cb}}}(\overrightarrow{r}-s^{\prime }\overrightarrow{\mathrm{\Omega }})d{s}^{\prime }}ds,\end{array}$$(C.4)

$$\begin{array}{c}\hfill {\mathcal{S}}_{\text{cb}}=\frac{\eta -i}{\eta }\frac{{\mathrm{\Sigma }}_{\eta }}{{\tilde{\mathrm{\Sigma }}}_{t^{\text{cb}}}}+\int \int \frac{⟨{\mathrm{\Sigma }}_{s{f}_s}⟩}{{\tilde{\mathrm{\Sigma }}}_{t^{\text{cb}}}}d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime },\end{array}$$(C.5)

and

$$\begin{array}{cc}\hfill {\mathcal{F}}_{\text{cb}}& =\frac{1}{4\pi k_{\text{eff}}}\int \int \frac{⟨{\chi }_p{\nu }_p{\mathrm{\Sigma }}_f⟩}{{\tilde{\mathrm{\Sigma }}}_{t^{\text{cb}}}}d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }\hfill \\ \hfill & +\frac{1}{4\pi k_{\text{eff}}}\sum _{j=1}^J\frac{{\lambda }_j}{{\lambda }_j+i\omega }\int \int \frac{⟨{\chi }_{d,j}{\nu }_{d,j}{\mathrm{\Sigma }}_f⟩}{{\tilde{\mathrm{\Sigma }}}_{t^{\text{cb}}}}d{E}^{\prime }d{\overrightarrow{\mathrm{\Omega }}}^{\prime }.\hfill \end{array}$$(C.6)

The extended counterparts of these operators are obtained, as explained in Section 4.1, by swapping the unitary complex numbers

$$\begin{array}{c}\hfill {\mu }_{\text{fb}}(s)=e^{-i\frac{\omega }{v}s},{\mu }_{\text{cb}}=\frac{\eta -i}{|\eta -i|},{\mu }_j=\frac{|{\lambda }_j+i\omega |}{{\lambda }_j+i\omega }\end{array}$$(C.7)

by the corresponding Z_μ kernels.

All Tables

All Figures