© G. Biot et al., Published by EDP Sciences, 2025

1. Introduction

Nuclear data encompass the fundamental properties of atomic nuclei and their interactions, which form the basis for simulating a wide array of nuclear applications. This data includes essential parameters such as nuclear cross-sections, decay modes, half-lives, and other basic nuclear properties. Therefore, nuclear data are an essential piece of information needed to understand all the physical processes that underlie most nuclear technologies. Applications of nuclear data include areas of nuclear science and technology such as fission energy, fusion energy, nuclear fuel cycles, waste management and decommissioning and accelerator driven systems. There are also non-energy applications such as medicine (diagnosis and therapy), production of radioisotopes for medical and industrial applications, dosimetry and radiation safety, nuclear security, materials analysis and astrophysics. Quantifying the uncertainties associated with nuclear data is critical for the reliability and safety of nuclear simulations. These uncertainties can significantly influence reactor physics parameters such as the neutron multiplication factor (k_eff), power distributions, and safety margins. Precise uncertainty quantification not only aids in optimizing reactor designs and enhancing safety protocols but also ensures the economic viability of nuclear technologies [1]. In addition, uncertainty quantification is crucial for the deployment of advanced nuclear systems, such as Generation IV and small modular reactors, fusion or even accelerator-driven systems. These technologies require high precision in nuclear data to ensure their safety, efficiency, and sustainability. Moreover, understanding nuclear data uncertainties can guide experimental efforts to improve the quality of nuclear data and refine theoretical models. Traditionally, two principal methods have been employed to propagate nuclear data uncertainties: sensitivity analysis using perturbation theory and the Total Monte Carlo (TMC) method. Perturbation theory utilizes sensitivity profiles and covariance data to linearly propagate input uncertainties to output quantities, providing an efficient but approximate deterministic approach. While widely adopted, this technique assumes linear relationships and depends on the availability of detailed covariance data. In contrast, the TMC method relies on random variations in the physical model parameters to generate perturbed nuclear data libraries. By performing independent simulations for each random realization and applying batch statistics, TMC bypasses the need for linearity assumptions and directly produces probability distributions for reactor parameters, offering a more comprehensive uncertainty characterization. Nowadays, hybrid techniques are used by generating random cross sections based on covariances and running Monte Carlo simulations for each random sample. Moreover, innovative methods using machine learning algorithms are emerging such as Lasso Monte Carlo presented in [2]. A new algorithm called Embedded Monte Carlo (EMC) was recently developed in [3] and presented in [4]. It is based on estimating the history moments to obtain the distribution of uncertainties of the simulation. The goal of this paper is to present the progress of this new algorithm on the two benchmarks used for validation and nuclear data uncertainty quantification by evaluating the impact of the fission source and spatial dependence on its performance.

2. Background

2.1. Monte Carlo statistics

In Monte Carlo codes, there are two different types of statistical estimators that are commonly implemented: history and batch. Assuming the neutrons are independent from each other, history and batch statistics are equivalent. The mean of a tally estimated using individual particle histories is equal to the mean obtained using batch statistics. However, while the history-based variance uses all available history data, the batch-based variance estimate uses the variance of the batch means composed of these same histories, leading to a reduction in degrees of freedom [5].

In a fixed source problem, neutrons are not correlated if the source is uncorrelated. This means that there is no difference between history and batch statistics beyond the number of samples used in estimating a sample variance. Batch statistics are typically preferred for computational performance and their ability to provide more rigorous confidence intervals. The Student-t distribution, commonly used for estimating confidence intervals, assumes that the samples are independent and identically distributed (i.i.d.) [6]. Since individual histories do not necessarily produce identically distributed tally contributions, grouping them in batches and relying on the Central Limit Theorem tends to produce a better distribution of tally contributions. As a result, the variance estimate aligns better with the i.i.d. assumption underlying the Student-t distribution, leading to a more reliable confidence interval estimate.

In eigenvalue problems, inter-cycle correlations and spatial dependencies arise from the stochastic fission process and the common fission source updating approach that links each generation of neutrons to the previous one. These correlations can hinder the statistical convergence of tallies, distort variance estimates, and compromise the reliability of uncertainties. It has been demonstrated in [7] that in systems such as a PWR, neutrons tend to cluster spatially due to the asymmetry between their uniform loss throughout the geometry and their localized birth at fission sites. These clusters drift over successive cycles, producing strong spatial and temporal fluctuations even after the multiplication factor has converged, making it difficult to reach a stationary flux distribution with low sample counts. While higher sample count will mitigate the appearance of clustering, it does not reduce the presence of correlation. The effect intensifies in larger systems, creating inter-cycle and spatially-correlated noise that violates the assumptions of independence required in commonly implemented variance estimators. It has been shown in [8] that the variance of tallies deviates significantly from the ideal 1/N scaling due to autocorrelation and that the deviation is more pronounced in larger tally regions.

It is important to recall that in the presence of correlation, the general and more accurate formula for the variance estimate of the mean is the following: [9]:

$$\begin{array}{c}\hfill \mathrm{Var}(\overline{x})=\frac{1}{N(N-1)}·[\sum _{i=1}^N\mathrm{Var}(x_i)+\sum _{i=1}^N\sum _{i\ne j}\mathrm{Cov}(x_i,x_j)]\end{array}$$(1)

where N corresponds to the number of independent realizations (particle histories or batches). The covariance of two different samples from the same population (thus μ being the mean for both samples) is defined as Cov(x_i, x_j)=𝔼[(x_i − μ)(x_j − μ)]. If the assumption that all batches are independent, then the second term in Equation (1) vanishes and this quantity becomes the unbiased variance estimator of the mean computed in Monte Carlo codes. This paper will not attempt to resolve the presence on inter-cycle correlations, but as the nuclear data uncertainty quantification process greatly relies on the use of variance estimators, it is important to understand when these methods may fail and lead to an overestimation of the nuclear data uncertainty caused by the underestimation of the statistical uncertainty.

3. Methodology

3.1. Total Monte-Carlo (TMC)

The uncertainty methodology with Monte Carlo is based on the use of random samples of nuclear data libraries and performing separate calculations for each random sample. This approach does not use any approximation, no first-order or linear approximation is required, but due to its stochastic nature, it carries this inherent statistical uncertainty which is not present when TMC is used with deterministic codes. The idea is to simulate a large number of calculations all alike but with different random nuclear data in each of them, the results consist of a probability distribution function from which different moments can be extracted. Two type of uncertainties arise from these moments. The first one is the aleatoric uncertainty, which results from the randomness of the calculation procedure performed by the Monte Carlo sampling process. This output uncertainty σ_statistical can be reduced or even eliminated by increasing the number of neutron histories or the number of batches in the simulation. The second one is the nuclear data uncertainty σ_nuclear-data that remains even with a high number of neutrons or batches. This method can be flexible in the use of covariance files but in [1], the TMC method generates random nuclear data from fundamental theoretical nuclear quantities with the help of a nuclear reaction code like TALYS [10]. For each random set of ENDF-6 files¹ and n simulations using these random samples, n different k_eff values or any other tallied quantity with their statistical uncertainty are obtained. If we assume that the samples are uncorrelated, the total variance of the simulation σ_total² can be separated into two parts [11]:

$$\begin{array}{c}\hfill {\sigma }_{\text{total}}^2\approx {\sigma }_{\text{statistical}}^2+{\sigma }_{\text{nuclear-data}}^2\end{array}$$(2)

where σ_statistical² corresponds to the statistical uncertainty derived from the number of histories used in the simulation and σ_nuclear-data² is due to the different random nuclear data files between calculations. It induces a spread in the distribution of the calculated reactor quantity, which can unequivocally be assigned to the spread of nuclear data (e.g. cross sections, angular distributions, nubar, double differential data, prompt fission neutron spectrum). The statistical uncertainty typically varies as $1/\sqrt{N}$ where N is the number of independent realizations (particle histories or batches) [12].

3.1.1. Fast TMC and GRS

One of the main drawback of Total Monte Carlo is the computation time required in comparison to a single calculation. Indeed, the minimization of the statistical uncertainty is making the uncertainty propagation inefficient compared to deterministic methods. The idea of fast TMC, defined in [13], is to obtain a Monte Carlo simulation with the same calculation time as a single run by reducing the number of histories for each run and using different seeds in each run. Uncertainties could then be provided for a Monte Carlo simulation within the same calculation time as a single simulation. The main differences between TMC and fast TMC are that the number of total histories is reduced for each sample and that each run has a different random seed in order to make sure there is no correlation coming from the pseudo-random number generator. It is important to note that this reduced number of histories does not impact the number of random samples since increasing this number does not lead to a more precise answer.

Another MC method for propagating faster nuclear data uncertainty is called fast GRS where two sets of simulations with different random seeds are compared [14]. Each set represents a collection of n simulations (corresponding to the number of random samples) performed with a small number of particle histories. The covariance between these two sets is equal to the variance of the distribution quantifying the epistemic part of sampling uncertainty alone after having eliminated the effect of aleatoric sampling uncertainty. However, because using very few neutron histories can cause bias or convergence issues, it is recommended to validate these results against a long, high-statistics simulation to ensure the reliability of uncertainty estimates.

3.2. Embedded Monte-Carlo (EMC)

The Embedded Monte Carlo (EMC) method rests on the observation that the output quantity of interest Z is an expectation value of the form Z(X)=𝔼_Y[Y|X], where X are input parameters and Y is the stochastic process being simulated (neutron random walk). EMC then samples different possible nuclear data parameter values X from their uncertainty distribution (e.g. normal distribution), and is able to converge towards accurate uncertainty estimates (Z(X) distribution) despite running only a few nuclear simulations with large variance estimates for each possible nuclear data sample (non-converged Z(X) for each X), thereby considerably reducing computational cost. EMC achieves this by estimating the moments of the output Z(X) distribution (forward propagation), and then reconstructing the distribution using the principle of maximum entropy², as explained in [3] and [15]. This means that instead of running a full MC simulation for each sample of the uncertain parameters distribution X to obtain directly the nuclear data uncertainties (histogram of Z(X)), the EMC method runs only one MC simulation assuming a fixed source approach where each batch has a different random sample of the nuclear data X sampling from the same fission source. The idea is to compute history statistics from each batch of neutrons. To fully grasp EMC, it is important to recall how all the estimators are computed. As in TMC, for each random sample M of input parameter X, N neutrons are simulated which correspond to N independent stochastic processes Y_n|X_m (X_m corresponds to the m^th random sample of parameter X, drawn from its uncertainty distribution). For each sample X_m, EMC estimates the mean $\widehat{Z}(X_m)$ and variance $\widehat{{\mathrm{\Sigma }}_Y}(X_m)$ between the N stochastic processes (also called “inner loop”) using Equations (3) and (4) respectively (equations and notation from [3]):

$$\begin{array}{cc}\hfill \widehat{Z}(X_m)& =\frac{1}{N}\sum _{n=1}^N{Y}_n|X_m\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill \widehat{{\mathrm{\Sigma }}_Y}(X_m)& =\frac{1}{N-1}\sum _{n=1}^N(Y_n|X_m-\widehat{Z}(X_m){)}^2.\hfill \end{array}$$(4)

The inner loop of the EMC algorithm thus runs the neutrons Y_n, while the outer loop samples the input parameters uncertainty distribution X_m. After running all M random samples, EMC computes the following estimators to find the nuclear data uncertainty (equations and notation from [3]):

$$\begin{array}{cc}\hfill ⟨\widehat{Z}⟩& =\frac{1}{M}\sum _{m=1}^M\widehat{Z}(X_m)\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill ⟨{\sigma }_Z^2⟩& =\frac{1}{M-1}\sum _{m=1}^M(\widehat{Z}(X_m)-⟨\widehat{Z}⟩{)}^2\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill ⟨\widehat{{\mathrm{\Sigma }}_Y}⟩& =\frac{1}{M}\sum _{m=1}^M\widehat{{\mathrm{\Sigma }}_Y}(X_m)\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill ⟨\widehat{{\mathrm{\Delta }}^{2}Z}⟩& =⟨{\sigma }_{\widehat{Z}}^2⟩-\frac{\widehat{⟨{\mathrm{\Sigma }}_Y⟩}}{N}\hfill \end{array}$$(8)

where $⟨\widehat{Z}⟩$ and $⟨\widehat{{\mathrm{\Delta }}^{2}Z}⟩$ estimate the mean and variance of the distribution of outputs Z(X)=𝔼_Y[Y|X]. These are unbiased estimators of the moments of the output quantity of interest, say k_eff = Z(X). $⟨\widehat{{\sigma }_Z^2}⟩$ represents the total variance, mixing both inner loop (aleatoric uncertainty from stochastic process Y) and outer loop (epistemic uncertainty from input parameter X). $⟨\widehat{{\mathrm{\Sigma }}_Y}⟩$ is the average of each history variance over the different input parameters. It is possible to see that Equation (6) is the same definition as σ_total² and that Equation (7) is identical to σ_statistics². Therefore, Equation (8) defines the way EMC computes the nuclear data uncertainty which is similar to the TMC method. Nevertheless, the main difference between TMC and EMC should be emphasized is that the mean of the statistical uncertainty in EMC is divided by the number of neutrons, not by the number of batches as defined in TMC. This can be explained by the fact that the statistical uncertainty is based on the mean variance of thepopulation. In MC codes such as OpenMC or MCNP, Equations (3) and (4) are computed based on the number of batches of the simulation and not the number of particle histories. The aim of EMC is able to get correct statistical estimates in order to run M samples with only a single batch of N neutrons for each sample.

4. EMC implementation in OpenMC

OpenMC is a community-developed Monte Carlo neutron and photon transport simulation code. It is capable of performing fixed source or k-eigenvalue multiplication calculations on models built using either a constructive solid geometry or CAD representation. OpenMC supports both continuous energy and multi-group transport. This open-source software simulates neutral particles (presentlyneutrons and photons) moving stochastically through an arbitrarily defined model that represents a real-world experimental setup. Knowing the behavior of neutrons allows one to determine how often and where reactions occur. By simulating many neutrons (millions or billions), it is possible to determine the average behavior of these neutrons (or the behavior of the energy produced, or any other quantity of interest) very accurately [16].

4.1. Practical challenges in EMC implementation

One limitation of implementing EMC as defined in OpenMC (i.e. using history statistics) is its increased memory usage for processing tallies. This issue comes from how EMC works: instead of accumulating a single running sum and sum of squares across all batches (which is memory-efficient), EMC needs the mean and variance within each individual batch. To achieve this, EMC requires at least two particles per batch and performs two independent tallies of the same observable given by Y₁ and Y₂, along with their squares Y₁² and Y₂². These are then used to compute the batch-level mean and variance, enabling estimation of both the observable’s central value and its statistical uncertainty within that batch. This method is especially powerful when each batch is run with different input parameters, such as perturbed nuclear data. In that context, EMC makes it possible to separate the total variance into the intrinsic Monte Carlo variance coming from the random walk of neutrons and the variance due to input parameter uncertainty such as uncertainties on the average neutron multiplicity. This ability to compute unbiased, per-batch variance estimates is the central contribution described in [3]. However, the added statistical resolution comes at a cost: the need to store and process multiple tallies per observable can significantly increase computational expense, even though the exact impact is not quantified in this study.

Another challenge in implementing EMC as originally defined arises from the parallelization architecture of OpenMC. This software supports both distributed-memory (MPI) and shared-memory (OpenMP) parallelism. When using OpenMP, multiple threads simulate particle histories independently, which is highly memory-efficient. However, this thread-level parallelism introduces a problem: when the function accumulate_tallies are called after each particle history, multiple threads may attempt to update the same tally bin simultaneously. This leads to race conditions where only one of the concurrent updates are applied, causing incomplete or incorrect tally accumulation. Specifically, the mean values obtained from history statistics are systematically lower, reflecting lost contributions during tally updates. A short-term solution is to run history statistics in single-threaded mode, which eliminates race conditions and ensures correct tally accumulation. A longer-term solution would involve modifying OpenMC’s code architecture to use atomic memory locks or tally data duplication. However, both options come with trade-offs in memory efficiency and parallel performance, highlighting the need for a different EMC implementation, which is presented inthis work.

4.2. Fission source convergence

Figure 1 presents the optimization scheme for EMC in OpenMC. The first step is to converge the fission source using the best estimate parameters. It is important to recall how the eigenvalue mode is structured in OpenMC. As many other MC codes, the fission bank stores the fission sites that will be sampled in the next generation to produce a source bank corresponding to the number of particles per batch. This new implementation stores the stationary fission bank at the end of the inactive cycles so that each active batch will have the same fission bank as a starting point. In this case, there is no positive correlations between the random samples but direct correlation from the initial stationary fission bank since all random samples are using the same initial distribution. This approach was deemed preferable since it insulates the estimation of nuclear data uncertainties from the spatial-correlation complexities.

Fig. 1. EMC scheme implemented in OpenMC.

In this implementation, an assumption is made that the fission source distribution as a small impact on nuclear data uncertainty. Indeed, the fission source is converged only once for EMC whereas it is converged for each random sample with TMC. Results will focus on comparing the output distributions of both fission sources and k_eff using these two test cases:

• Perform standard TMC and look at the distributions for both fission sources and k_eff.

• Converge the fission source using unperturbed nuclear data and load it as a source for each random sample running a full MC simulation and compare k_eff with the previous case.

In this latter case, two observations should be verified: the fixed fission source is close to the mean of the fission source distribution obtained from the first case and that the uncertainty distribution on k_eff for the second case is similar to the full TMC simulation.

4.3. Statistical uncertainty estimate

Once the fission source is converged, one active cycle is performed to compute the variance of the statistical process for the fission source obtained from the best estimate nuclear data. The estimate of the statistical uncertainty is made using the first active batch with the generations_per_batch capability of OpenMC. Indeed, the assumption made here is that the statistical uncertainty should be similar across different batches as long as the same number of particles is run.

It is essential to validate the assumption of using the same statistical uncertainty by examining the potential impact of nuclear data perturbations on the statistical uncertainty. In the TMC method, this uncertainty is computed as the average of the standard deviations from each simulation. Consequently, it is possible to analyze the distribution of the standard deviations from each random sample to assess whether the nuclear data perturbations influence the statistical uncertainty. In Equation (2), the two quantities are clearly separated; however, a sanity check is presented in Section 6.1.1 to ensure that this assumption holds true in these benchmarks.

4.4. Total uncertainty estimate

Instead of directly computing history-based means, the idea is to extract and store the mean values of k_eff at the end of each batch to compute the total uncertainty. One such modification involves dynamically loading new HDF5 files for each new batch during the simulation, which has now been successfully implemented in OpenMC. This is accomplished by leveraging the openmc.lib module, which binds the Python API to the underlying C++ functions defined in the OpenMC shared library, enabling the loading of new nuclides before the start of each active batch. Consequently, it becomes possible to utilize a unique set of nuclear data for each batch within a single Monte Carlo simulation. This feature can also be extended to materials, facilitating changes to the density of the fission isotope. Moreover, this capability holds significant promise for depletion calculations, where the material composition evolves over time while the geometry, tallies, and simulation settings remain unchanged.

It is important to note that the eventual goal of the EMC method remains to provide the capability of sampling nuclear data using covariance data on-the-fly before each batch. Indeed, HDF5 files do not contain the covariance data necessary for creating random samples. And even if the covariance data were present in the files, the processing code NJOY [17] would still be required to handle and assemble the complex elements of covariance data in the ENDF files (combining different sections of the covariance matrix from different energy grids, handling absolute and relative uncertainties or even combining uncertainties from cross sections MF33 and from resonance parameters MF32). Future work will look at integrating this approach with the windowed multipole method, and also developing data models that can represent data parameters that can be easily sampled for other nuclear parameters of interest. While not yet complete, this EMC implementation is a first step in enabling us to obtain nuclear data uncertainty using a single Monte Carlo simulation with a reduced number of particles.

5. The analytic benchmark for neutron transport with downscattering

This benchmark provides an analytically solvable reference problem for neutron transport in energy, addressing the critical challenge of accurately modeling neutron slowing down through resonance energy ranges [18]. By considering an infinite homogeneous medium with isotropic down-scattering, and using simplified cross sections for ²³⁹Pu, ²³⁸U and ¹H (scattering isotope), it is possible to derive closed-form solutions for both the time-dependent and steady-state Boltzmann transport problems. This benchmark provides exact expressions for neutron flux and reaction rates that capture the effects of resonance self-shielding typically resolved only through costly Monte Carlo simulations or approximate multi-group methods. This benchmark serves as a powerful tool to verify and validate the energy resolution, flux accuracy, and eigenvalue convergence of modern neutron transport solvers, particularly in criticality calculations.

In the original analytic benchmark, the average neutron multiplicity was defined as constant with $\overline{\nu }=2.88$ at all neutron energies whereas the PFNS was defined as uniform with equal probability of emitting neutrons at any energy. Indeed, the PFNS was defined as independent of incident and outgoing neutron energy using the unit step function (also called the Heaviside function) as defined in [19]:

$$\begin{array}{c}\hfill {\chi }_i={\chi }_{i,0}·H(E_0)+\end{array}$$(9)

$$\begin{array}{cc}& ({\chi }_{i,1}-{\chi }_{i,0})·H(E_1)+\dots +({\chi }_{i,N}-{\chi }_{i,N-1})·H(E_N)\hfill \end{array}$$(10)

where i and j are the indices over the incident and outgoing energy range, respectively and χ_i, j is the prompt fission spectrum for an incident neutron between E_i and E_i + 1 producing an outgoing neutron in the range between E_j and E_j + 1. Recently, a more realistic representation was produced for two parameters with covariance data based on the evaluations of ENDF/B-VIII.0 [20]. These ENDF files with covariance data for $\overline{\nu }$ and χ are used for testing EMC for nuclear data uncertainty quantification. Table 1 presents the values for $\overline{\nu }$ and the uncertainties $\delta \overline{\nu }$ with four energy groups going from 10⁻⁵ eV until 20 ⋅ 10⁶ eV. For the PFNS, the uncertainties were too small to be able to propagate them with any noticeable effect on the analytical benchmark. The impact of the prompt fission spectrum is very small on k_eff since the values for χ were starting at 100 eV which corresponds to a higher energy than the capture resonance of the benchmark. A non-physical PFNS distribution was built to obtain uncertainties that could be verified with a Monte-Carlo simulation for this test case. This distribution is made using two-group PFNS with equal numbers of neutrons produced below and above the capture resonance and will be used in this work.

6. Verification and results

6.1. The analytic benchmark

Based on these ENDF files, it is possible to propagate the uncertainties from the neutron multiplicity and the prompt neutron fission spectrum. These files represent the three isotopes given in the analytical benchmark with the fission (first resonance in ²³⁹Pu), capture (first resonance in ²³⁸U) and scattering isotopes. They were processed using NJOY to create ACE files that were converted to the OpenMC HDF5 format using the Python module openmc.data. For the density of ²³⁹Pu, since there is no covariance data available, a Gaussian perturbation is applied to the model and the modifications were made directly in the material.xml file for each batch. For the average neutron multiplicity and the prompts fission spectrum, a set of random HDF5 file were pre-generated using the covariance data and stored with a cross_sections.xml file. The OpenMC input is composed of a homogeneous mixture with reflective boundary conditions to represent the infinite medium. Table 2 presents the number density ratios of the analytic benchmark used in the simulations. The simulations were performed using 32 CPU cores (16 cores per socket with 64-bit CPU architecture) with OpenMC working in eigenvalue mode for each random sample. The output parameter k_eff is used to extract the moments of the distributions in order to compute the nuclear data uncertainty. For this benchmark, only the comparison between TMC and EMC is performed, since no spatial variation is present due to the benchmark’s infinite and homogeneous configuration.

6.1.1. Statistical uncertainty check

Before testing EMC, it is important to check that the assumption made in Section 4.3 is valid. Table 3 presents the mean and standard deviation of the statistical uncertainties computed from 500 random samples using the TMC method for the analytical benchmark. For the density and nubar cases, the average runtime is around 226 minutes, while the PFNS and combined cases were completed in roughly half that time. This difference can be explained by the number of energy groups used for the PFNS in the different cases as described in Section 5. These results demonstrate that performing simulations with the same number of particles leads to very similar standard deviations, regardless of which nuclear data is perturbed. Consequently, the same number of particles is employed in both the fast TMC and EMC methods to ensure a consistent basis for comparing the nuclear data uncertainty and computational efficiency between the two approaches. Additionally, it is important to recall that the statistical uncertainty σ_stat² should not exceed 5% of the total uncertainty σ_tot², as outlined in [13].

6.1.2. Average neutron multiplicity

For the average neutron multiplicity, 250 random samples were produced with 250 active batches and 10⁵ particles per batch for the Total Monte Carlo method. Figure 2a presents the neutron multiplication factor distribution when perturbing the average neutron multiplicity and three normal distributions with moments computed by different methods: TMC – blue dashed line, EMC – black dotted line, and in red the sensitivity of the k-eigenvalue to $\overline{\nu }$ computed analytically in [19] using the adjoint flux. The histogram shows the output distribution generated from 250 random samples using the TMC method, which serves to validate the calculated uncertainty on the neutron multiplication factor for the extended benchmark. Good agreement can be observed for the mean values across all three methods with differences well within the statistical uncertainty as seen in Table 4.

Fig. 2. keff distribution for the analytic benchmark for $\overline{\nu }$, PFNS, density of 239Pu and all parameters combined using 500 random samples and respectively 108 and 5 ⋅ 106 neutron histories for TMC and EMC respectively. (a) Average neutron multiplicity simulation. (b) PFNS simulation. (c) Density of 239Pu case. (d) All input parameters simulation.

6.1.3. Density of ²³⁹Pu

For the density of ²³⁹Pu, a random distribution is used to create random samples with mean μ = 1.0 and σ = 0.005 for all methods. For TMC, 10³ batches with 10⁵ neutrons per batch are used whereas only 5 ⋅ 10⁶ neutrons are used for the EMC method. Only the mean and total uncertainty are computed for each value of $\overline{\nu }$ and the reference statistical uncertainty from the first active batch is used since the same number of particles is used in all batches. It is important to note that there is no sensitivity computed for the density of the fission isotope since there is no uncertainties on this parameter coming from the extended benchmark. One could note that the mean values and nuclear data uncertainties are in good agreement between the two methods with differences within the statistical uncertainty.

6.1.4. Prompt fission neutron spectrum

In this case, 400 random samples with 10³ active batches and 10⁵ particles are used with the TMC method. Using the non-physical PFNS, Figure 2b presents the k_eff distribution with the nuclear data uncertainty computed with the same three different methods described in Section 6.1.2. It can be observed that both methods agree once again well with the deterministic value.

6.1.5. Combined simulation

Finally, the last case considers perturbing all the input parameters, as illustrated in Figure 2d. The agreement between TMC and EMC is once again within the statistical uncertainty. Furthermore, the density appears to have the largest impact on the combined value, as the mean values, given in Table 4, are very similar to those obtained by perturbing only the density of the fission isotope, but the largest contribution from uncertainty comes from neutron multiplicity. This observation is consistent with the nature of this benchmark, which relies entirely on number density ratios to ensure an exact k_eff of 1.00 for the physical PFNS data. One interesting point is to study the possibility of adding up the uncertainties coming from different parameters such that:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{ND}}^2={\sigma }_{\overline{\nu }}^2+{\sigma }_{\chi }^2+{\sigma }_{\rho }^2.\end{array}$$(11)

With the results obtained in this analytic benchmark, it is possible to show that for TMC, there is a 3.5 pcm difference between the results obtained in Table 4 and with Equation (11) which is within the statistical uncertainty. One should note that adding the uncertainties individually up is valid as long as the uncertainty sources are established to be independent [21]. Concerning the simulation time, it is possible to observe that for EMC, the average simulation time is about 34 minutes for the density and nubar cases and about 16 minutes for the PFNS and combined cases. The analytic benchmark does not need a lot of cycles for the fission source to converge since it is an infinite homogeneous mixture of three isotopes thus only providing a reliable energy distribution of the source sites. This benchmark serve thus as a verification of the implementation, but additional benchmarking is necessary to validate other assumptions.

Fig. 3. keff distribution for the spherical benchmark for $\overline{\nu }$, χ, density of 239Pu and all parameters combined using 400 random samples and respectively 107, 106 and 106 neutron histories (except for PFNS) for TMC, fast TMC and EMC respectively. (a) Average neutron multiplicity simulation. (b) PFNS simulation. (c) Density of 239Pu case. (d) All input parameters simulation.

6.2. The spherical benchmark

Since the analytic benchmark does not have spatial variations, it is desirable to evaluate EMC on a more realistic benchmark such as a finite spherical benchmark similar to the Godiva benchmark [22]. The OpenMC input is composed of a homogeneous mixture of ²³⁹Pu, ²³⁸U and ¹H with vacuum boundary conditions and modified number density ratios to obtain criticality in this benchmark. The output parameter k_eff is used to extract the moments of the distributions in order to compute the nuclear data uncertainty.

6.2.1. Density of ²³⁹Pu

For the density of ²³⁹Pu, a random distribution is employed to generate random samples, with a mean of μ = 1.0 and two different standard deviations are used: σ₁ = 0.001 and σ₂ = 0.01, across all methods. The second case is performed to test the limits of EMC even if 1% error on the density of the fission isotope is not realistic since it will induce a large change in the curvature of the solution. For the first case, a total of 400 random samples, using 500 active batches and 10⁵ neutrons per batch, are used for the TMC method. Figure 3c shows the output distribution with the normal distribution computed using three methods: TMC with the red line, EMC with the black dotted line, and fast TMC with the dashed blue line. Table 5 presents the results obtained for the three methods using σ₁ = 0.001. Due to the low level of uncertainty for this case, the number of particles required for EMC and fast TMC simulations can be reduced by only one order of magnitude. Despite this reduction, all methods agree well with results well within statistical uncertainty.

Fig. 4. keff distribution for the second case of 239Pu and all parameters combined using 400 random samples and respectively 108, 106 and 106 neutron histories for TMC, fast TMC and EMC respectively. (a) Density of 239Pu case. (b) All input parameters simulation.

6.2.2. Average neutron multiplicity

As for the average prompt neutron multiplicity, 400 random samples are used with 1000 active batches and 10⁵ neutrons per batch for the reference case. Figure 3a illustrates the nuclear data uncertainty computed using the two different methods. Due to the very small ratio of $\frac{{\sigma }_{\mathrm{stat}}^2}{{\sigma }_{\mathrm{tot}}^2}=0.036$, it is possible to reduce the number of neutrons by two orders of magnitude for EMC and fast TMC, reaching a ratio of $\frac{{\sigma }_{\mathrm{stat}}^2}{{\sigma }_{\mathrm{tot}}^2}=0.33$ with only 10³ neutrons per batch. The results are also in agreement with the TMC results and within the bounds of the statistical uncertainty terms.

6.2.3. Prompt fission neutron spectrum

For the PFNS, 400 random samples with 1000 active batches and 10⁵ neutrons per batch for the reference case. Similar to the previous case, the ratio between the statistical and the total uncertainties is very small $\frac{{\sigma }_{\mathrm{stat}}^2}{{\sigma }_{\mathrm{tot}}^2}=0.038$, it is possible to reduce the number of neutrons by two orders of magnitude for EMC and fast TMC while still showing good agreement with the TMC results.

6.2.4. Combined simulation

For the combined case, 400 random samples with 1000 active batches and 10⁵ neutrons per batch for the reference case. Similar to the previous case, since the ratio between the statistical and the total uncertainties is very small $\frac{{\sigma }_{\mathrm{stat}}^2}{{\sigma }_{\mathrm{tot}}^2}=0.021$, it is possible to reduce the number of neutrons by two orders of magnitude for EMC and fast TMC to get a ratio of $\frac{{\sigma }_{\mathrm{stat}}^2}{{\sigma }_{\mathrm{tot}}^2}=0.17$ yielding once again good agreement across the methods as shown in Figure 3d. To fully compare the different methods, Table 6 presents the simulation time of each method for the perturbed input parameters.

It is possible to observe a significant difference in total simulation time between TMC and the two other methods. This can be explained by the reduction in particle count by one or two orders of magnitude for fast TMC and EMC depending on the parameter being perturbed. It is important to note all the k_eff values are higher than 1.0 since all the simulations were made using the non-physical data representing the PFNS for consistency. Indeed, the original ENDF file contained a more physical representation of the PFNS with 23 outgoing energy groups and an average energy of 2.1228 MeV but couldn’t be used for all cases.

One current drawback of EMC is the convergence of the fission source since OpenMC is running N particles with B generations per batch which is not necessary in this benchmark taking more than half of the total simulation time. From Table 6, it is possible to see that fast TMC can perform better than EMC by obtaining uncertainties on nuclear data two times faster when accounting for the inactive cycles of EMC. To address this, OpenMC will include an option to modify the number of particles per generation per batch during inactive cycles, enabling faster simulations for EMC. Nevertheless, for problems with fast fission source convergence and minimal spatial heterogeneity, EMC could become more computationally efficient than fast TMC, which requires fission source convergence for each random sample.

6.2.5. Large density perturbation

A second case is ran with an uncertainty of 1% on the fission isotope density to test the limits of the EMC method. Table 7 presents the results obtained for the density and the combined cases for the three different methods. It is possible to see significant discrepancies between the TMC and the EMC methods even for the total uncertainty. This difference can be explained by the fact that changing the density of the fission isotope is modifying the fission source substantially. Since the new fission source is not converged at each perturbation, the estimation of the nuclear data uncertainty is not correct as observed in Figure 4. This test demonstrates the impact of fission source changes on the current EMC implementation. One possibility would be to make EMC look even more like fast-TMC by adding inactive batches at each perturbation or implementing a super-history approach that would extend the history of each source site. Both of these options would impact runtime. Future work will look at the impact of source convergence on a suite of realistic benchmarks and at ways to evaluate the uncertainty related to spatial variations.

7. Conclusion

To enable the advancement of future nuclear technologies such as space exploration, medical diagnostics, and carbon-neutral power generation, it is essential to perform highly accurate nuclear simulations supported by reliable confidence intervals. At the core of these simulations lie nuclear data, which originate from experimental measurements and theory and inherently carry uncertainty. Properly accounting for these uncertainties are critical to ensure both the accuracy and safety of nuclear systems. In fact, nuclear data uncertainty has become a limiting factor in the development of innovative reactor technologies, including Small Modular Reactors (SMRs), Generation IV systems, and fusion reactors. The Total Monte Carlo (TMC) method and its accelerated variant, fast TMC, are state-of-the-art techniques for nuclear data uncertainty quantification. By generating a large ensemble of randomly sampled nuclear data files, these methods allow for a comprehensive assessment of how uncertainties propagate through nuclear simulations, thereby enhancing system reliability and safety. This study contributes to the development of the Embedded Monte Carlo (EMC) method by evaluating its performance against two benchmarks in eigenvalue simulations. The results indicate that EMC performs well in propagating uncertainties for the average neutron multiplicity, the prompt fission neutron spectrum, the density of the fissile isotope used in these studies (²³⁹Pu). Among these, perturbing the average neutron multiplicity has the biggest impact on the multiplication factor. Many improvements are still needed to optimally integrate EMC in OpenMC, but this work provides an initial path at evaluating key assumptions needed as well as the development of reading new nuclear data libraries at each batch. Moreover, higher moments will be computed in order to correctly estimate the distribution’s tail behavior and for assessing the likelihood of extreme values. Indeed, EMC is able to compute the skewness (third moment) and the kurtosis (fourth moment) of the output distribution in order to obtain correct and unbiased estimators of the output distribution as defined in [3]. This study shows the possibility of integrating nuclear data uncertainty quantification in a single run. It should also be noted that fast-TMC and EMC are very similar in nature and that relaxing some of the stated assumptions of EMC will bring them even closer.

Funding

This research is being performed using funding received from the DOE Office of Nuclear Energy’s Nuclear Energy University Programs.

Conflicts of interest

The authors have nothing to disclose.

Data availability statement

This article has no associated data generated and/or analyzed/Data associated with this article cannot be disclosed due to legal/ethical/other reason.

Author contribution statement

Methodology, G. Biot and P. Ducru and B. Forget; Validation, A. Lewis and V. Sobes; Formal Analysis, G. Biot; Resources, A. Lewis and V. Sobes; Data Curation, A. Lewis and V. Sobes; Writing – Original Draft Preparation, G. Biot; Writing – Review & Editing, B. Forget and P. Ducru and A. Lewis; Supervision, B. Forget and P. Ducru; Funding Acquisition, B. Forget.

References

All Tables

All Figures

	Fig. 1. EMC scheme implemented in OpenMC.
In the text

	Fig. 2. keff distribution for the analytic benchmark for $\overline{\nu }$, PFNS, density of 239Pu and all parameters combined using 500 random samples and respectively 108 and 5 ⋅ 106 neutron histories for TMC and EMC respectively. (a) Average neutron multiplicity simulation. (b) PFNS simulation. (c) Density of 239Pu case. (d) All input parameters simulation.
In the text

	Fig. 3. keff distribution for the spherical benchmark for $\overline{\nu }$, χ, density of 239Pu and all parameters combined using 400 random samples and respectively 107, 106 and 106 neutron histories (except for PFNS) for TMC, fast TMC and EMC respectively. (a) Average neutron multiplicity simulation. (b) PFNS simulation. (c) Density of 239Pu case. (d) All input parameters simulation.
In the text

	Fig. 4. keff distribution for the second case of 239Pu and all parameters combined using 400 random samples and respectively 108, 106 and 106 neutron histories for TMC, fast TMC and EMC respectively. (a) Density of 239Pu case. (b) All input parameters simulation.
In the text