# Towards a highly efficient and unbiased population-control algorithm for kinetic Monte Carlo simulations

Cecilia Montecchio, Vincent Lamirand, Davide Mancusi, Wilfried Monange, Andrea Zoia

<figure><img src="image_1.png" /></figure>

## To cite this version:

Cecilia Montecchio, Vincent Lamirand, Davide Mancusi, Wilfried Monange, Andrea Zoia. Towards a highly efficient and unbiased population-control algorithm for kinetic Monte Carlo simulations. International Conference on Supercomputing in Nuclear Applications + Monte Carlo, SNA + MC 2024, Oct 2024, Paris (France), France. pp.09006, 10.1051/epjconf/202430209006. hal-04889126

HAL Id: hal-04889126

https://hal.science/hal-04889126v1

Submitted on 15 Jan 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Towards a highly efficient and unbiased population-control algorithm for kinetic Monte Carlo simulations

Cecilia Montecchio<fnref n="1" />,*, Vincent Lamirand<fnref n="2" />,** Davide Mancusi<fnref n="1" />,***, Wilfried Monange<fnref n="3" />,****, and Andrea Zoia<fnref n="1" />,†

**Abstract.** Population-control methods are key to non-stationary Monte Carlo simulations of multiplying systems: they prevent either the unbounded growth or the disappearance of neutrons, occurring respectively in supercritical and subcritical conditions; furthermore, they contribute to an efficient allocation of computational resources by addressing the unbalance between the neutron and the precursor populations. In this paper, we present two alternative population-control algorithms: the legacy implementation in TRIPOLI-4®, the Monte Carlo code developed at CEA, and an improved version that is currently under investigation, based on the use of a simplified point-kinetics solver. We assess the performance of these methods through the simulation of a $2.2$ step reactivity insertion in a fast system (Flattop-Pu), leading to an increase of the neutron population by a factor 200, which is benchmarked against point kinetics. We show that the new implementation not only suppresses the slight bias that was present in the legacy method due to a stochastic normalization factor, but also outperforms the previous algorithm in terms of variance reduction and improvement of the figure of merit.

## 1 Introduction

Thanks to the increasing availability of computational resources and the development of ingenuous variance-reduction techniques, Monte Carlo codes have recently broadened their range of application and can now address non-stationary neutron-transport problems, which historically have been handled solely by deterministic methods. Pioneering work in this domain was conducted by Sjenitzer and Hoogenboom, who were the first to demonstrate the feasibility of simulating nuclear reactor kinetics by using Monte Carlo methods [1]. Kinetic simulation capabilities have recently been integrated also in the code TRIPOLI-4® [2]. For a

<footnote n="1">¹Université Paris-Saclay, CEA, Service d'Etudes des Réacteurs et de Mathématiques Appliquées, 91191, Gif-sur-Yvette, France</footnote>

<footnote n="2">² École polytechnique fédérale de Lausanne (EPFL), Switzerland</footnote>

<footnote n="3">³ Institut de Radioprotection et de Sûreté Nucléaire, 31 avenue de la Division Leclerc, Fontenay-aux-Roses, France 92260</footnote>

<footnote>*e-mail: cecilia.montecchio@cea.fr</footnote>

<footnote>**e-mail: vincent.lamirand@epfl.ch</footnote>

<footnote>***e-mail: davide.mancusi@cea.fr</footnote>

<footnote>****e-mail: wilfried.monange@irsn.fr</footnote>

<footnote>†e-mail: andrea.zoia@cea.fr</footnote>
comprehensive overview of the computational strategy adopted in this code to solve the time-
dependent Boltzmann equation, coupled with the precursor balance equations, the reader is
referred to the description provided by Faucher et al. [3]. Extensive tests of TRIPOLI-4®'s ki-
netic capabilities have been carried out in various configurations, often considering thermal-
hydraulics feedback, with satisfactory performance recorded, including inter-code compar-
isons [4–6]. However, kinetic Monte Carlo codes still face some notable limitations, mainly
concerning the simulation time and the statistical accuracy of the results, which severely hin-
der the simulation of severe transients (e.g. high reactivity insertions) in large reactor cores.
One of the toughest challenges encountered, which is the focus of our paper, is induced by
the imbalance among the sizes of the neutron and precursor populations due to their largely
different respective timescales [1]. Indeed, the characteristic timescale of the precursor pop-
ulation (i.e. the average decay time $1/\lambda$, ranging from 0.1 s to 10 s) is much longer than the
typical neutron timescale (i.e. the mean generation time $\Lambda$, approximately 25 µs in a typical
PWR). This is mirrored in the ratio $\beta/(\lambda\Lambda)$ between the two populations at critical conditions:
since the delayed neutron fraction $\beta$ is approximately 0.7% for a thermal $^{235}$U-based system,
a snapshot of the core at equilibrium would show $10^4$ more precursors than neutrons [1].

This substantial difference (which would be even larger for fast, Pu-based systems, with smaller $\beta$ and $\Lambda$), combined with the potentially unbounded growth or disappearance of the particles during a transient, has cogent consequences for kinetic Monte Carlo simulations. Since the number of particles that can be handled in a simulation is necessarily smaller than those circulating in a reactor core, extreme care must be taken to manage the size and the evolution of the sampled populations while also addressing the disparity between precursors and neutrons. This task is typically achieved using population-control methods, which aim at preventing either the unbounded growth or the disappearance of the particles, and at more efficiently allocating the computational resources in order to improve the Figure of Merit (FOM) of the simulations. In this work we illustrate the legacy method implemented in TRIPOLI-4® and a new method that is currently under investigation: the merits and drawbacks of these strategies are assessed on the Flattop-Pu benchmark configuration [7].

This paper is organized as follows: in Sec. 2, we present the two population-control strate-
gies. Section 3 outlines the benchmark used to assess the performance of these algorithms.
Section 4 presents the outcomes of applying the algorithms and compares them with the time
evolution of the point-kinetics model. Special emphasis is given to the fact that the new al-
gorithm is unbiased, contrary to the legacy one. Conclusions and perspectives are addressed
in Sec. 5.

## 2 Population-control algorithms

In this section, we describe two algorithms for controlling the neutron and the precursor
populations in a kinetic Monte Carlo simulation. Both algorithms consist of two phases,
which are consistent in essence and differ only in their application. They are applied during
the kinetic phase of the simulation, at each time $t_i$ of the time grid provided by the user, after
the system has been prepared in its initial (typically critical) state. These phases are:

i) *The scaling of the neutron and precursor populations according to the definition of the population importance ratio.* This importance-sampling scheme allows managing the unbalance between neutrons and precursors stemming from the different timescales. The idea behind this phase is straightforward: although precursors are more abundant, simulations should prioritize neutrons as they usually are the ‘carriers’ of the responses of interest (i.e. reaction rates, energy deposition, etc.), at least if the observation time scale
is short. This involves rescaling the neutron and precursor weights by appropriate factors, and compensating for these factors when tallying the sought scores. By defining a global importance value for the neutron ($I_n$) and the precursor ($I_c$) populations, denoting the ratio between the two as $R = I_n/I_c$, we can replace the physical particle weights ($w$) by the simulation weights ($\tilde{w}$), as follows [8]:

$$ \tilde{w}_n = w_n \cdot I_n \qquad (1a) $$

$$ \tilde{w}_c = w_c \cdot I_c. \qquad (1b) $$

The evaluation of $I_n$ and $I_c$ differs between the two algorithms. In both cases, the user is required to provide an importance ratio $R$. Although the choice of $R$ is arbitrary, it has been shown that there is an optimal value that maximises the FOM and depends on the size of the time-step of the scoring time grid [6].

ii) The application of the combing technique to control the actual population size (i.e. the number of particles transported within a time-step $\delta t_i$ of the simulation). This method, initially proposed by Booth [9] for stationary problems, involves sampling $N$ particles out of a collection of $K$ particles while preserving the total statistical weight. This population resizing is very effective when applied frequently during kinetic Monte Carlo simulations, because it preserves a convenient amount of histories to be simulated, thus preventing catastrophic divergence or shrinking of the population [1, 3]. Following the population resizing process, we assign each of the $N$ newly generated particles an equal weight. Further adjustment of this weight is achievable by performing an additional normalization of the total simulation weight by a global factor, however this passage depends on the algorithm and will be more thoroughly explained below. For completeness, it is worth noting that an alternative technique to control the population size, based on Russian roulette, can be applied specifically to the precursors [3]. However, the two algorithms described in this paper focus on the use of combing, which has consistently provided better results in previous investigations [6].

## 2.1 Legacy algorithm

We address first the algorithm that is currently implemented in TRIPOLI-4®. For the importance-sampling scheme described in point i), we rescale the neutron and the precursor physical weights by the importance values at the beginning of the kinetic simulation, after the fundamental distribution of neutrons and precursors has been determined by power iteration. The values $I_n$ and $I_c$ are estimated according to the following expressions [6]:

$$ I_n = R \cdot \frac{w_{n,tot} + w_{c,tot}}{R \cdot w_{n,tot} + w_{c,tot}} \qquad (2a) $$

$$ I_c = \frac{w_{n,tot} + w_{c,tot}}{R \cdot w_{n,tot} + w_{c,tot}}, \qquad (2b) $$

which can be easily derived by imposing the equality among the sum of the physical weights and the sum of the simulation weights:

$$ w_{n,tot} + w_{c,tot} = \tilde{w}_{n,tot} + \tilde{w}_{c,tot}. \qquad (3) $$

The definition of the importance functions as in Eq. (2) is rather straightforward, but may conceal a bias due to the fact that the computed importance values are non-linear functions of the total simulation weights, which are random variables. The normalization by a random
variable, characterized by statistical fluctuations, is an operation that, when used in a Monte Carlo implementation, is known to possibly introduce a bias [10]. The occurrence of this bias has been extensively demonstrated e.g. in the context of the estimation of the effective multiplication factor $k_{\text{eff}}$ with the power iteration algorithm [10, 11]. Therefore, when computing the simulation weights according to the corrections outlined in Eqs. (1), the bias introduced by the importance factors is embedded in the resulting transported weights.

Concerning the second phase, described by ii), combing is applied in such a way that the population size is brought back to the initial size and each particle is assigned a weight equal to one, which is the most convenient transport condition. This results in the creation of $N$ particles with a total weight equal to $N$. In order to achieve this, as already anticipated in ii) within Sec. 2, the simulation weights of the newly generated particles undergo an additional normalization, later compensated during the tallying procedure. This normalization, applied at the beginning of every time-step $\delta t_i$ of the time grid, involves the following factor:

$$f_{t=t_i} = \frac{\tilde{w}_{tot,t=t_i}}{N} \qquad (4)$$

It is apparent that Eq. (4) is significantly affected by statistical fluctuations and varies from replica to replica for every time-step $\delta t_i$. Therefore, this phase of the population-control algorithm implies a normalization by a random variable, and thus also introduces a bias in the tallies, similar to the one already introduced with the importance-sampling scheme.

## 2.2 New algorithm

In view of the previous considerations, a new algorithm was designed with the primary goal of achieving an unbiased kinetic simulation. However, it should also be regarded as a refinement of the legacy method, encompassing additional features, as will be extensively demonstrated in Sec. 4. Concerning the first phase, which is the application of the importance-sampling scheme, the importance functions are defined in a similar way as the ones presented in Eq. (2). However, instead of using the simulation weights as estimates of the population sizes, we use the solutions of a point-kinetics model. The modified definition is summarized by the following set of equations:

$$\frac{dn}{dt} = \frac{\rho - \beta_{\text{eff}}}{\Lambda_{\text{eff}}} \cdot n + \bar{\lambda} \cdot c \qquad (5a)$$

$$\frac{dc}{dt} = \frac{\beta_{\text{eff}}}{\Lambda_{\text{eff}}} \cdot n - \bar{\lambda} \cdot c \qquad (5b)$$

$$I_n = R \cdot \frac{n + c}{R \cdot n + c} \qquad (5c)$$

$$I_c = \frac{n + c}{R \cdot n + c}, \qquad (5d)$$

where Eqs. (5a) and (5b) are the point-kinetics equations, with a single precursor family; $n$ and $c$ are respectively the neutron population size and the precursor population size, and the reference values of the kinetics parameters are estimated by TRIPOLI-4® using the IFP method [12]. When solving the equations of the point-kinetics model, the effective kinetics parameters $\Lambda_{\text{eff}}$ and $\beta_{\text{eff}}$, as well as the average decay time $\bar{\lambda}$, are treated as time-independent, while the (static) reactivity is computed as $\rho = 1/k_{\text{critical}} - 1/k_{\text{perturbed}}$. Eqs. (5c) and (5d) are used to compute the importance values within the importance-sampling scheme, and play the same role as Eqs. (2a) and (2b). The quantity $R$ is again the user-defined importance ratio, similarly as in Sec. 2. The main difference with the legacy formulation is that in this revised
version, instead of the (stochastic) total weights provided by the simulation, we use the (non-stochastic) solutions of the point kinetics, which can be seen as ‘reasonable’ deterministic estimates. This strategy avoids the normalization by a random variable in the definition of the correction factors, which prevents the emergence of a bias. In this new algorithm, the definition of the simulation weights according to Eq. (1), hence the correction based on the importance values in Eq. (5c) and (5d), is performed not only at the beginning of the kinetic simulation, as done in the legacy algorithm, but also at the start of every time-step $\delta t_i$. Indeed, the importance values are computed at the beginning of every time-step $\delta t_i$, using the values of n and c stemming from the analytical solutions of the point-kinetics equations at $t = t_i$. Therefore, with this new algorithm, the weights are progressively adjusted during the transient, following the scheme

$$w_{n,t_{i+1}} = w_{n,t_i} \cdot \frac{I_{n,t_{i+1}}}{I_{n,t_i}} \qquad (6a)$$

$$w_{c,t_{i+1}} = w_{c,t_i} \cdot \frac{I_{c,t_{i+1}}}{I_{c,t_i}}, \qquad (6b)$$

where $t_i$ and $t_{i+1}$ are the starting points of two subsequent time-steps. This correction does not introduce any bias, since the adjustment of the particles weights relies only on non-stochastic quantities.

Concerning phase ii) in Sec. 2, the application of combing provides only an additional refined control, since by using the new algorithm the weights are already adjusted at every time-step according to the scheme described above. In this case, the population size is restored to the original value N and we assign to the newly generated particles a weight equal to $w_{tot}/N$. The additional normalization to achieve a unit weight is not performed and the total simulation weight $w_{tot}$ is preserved, because it is expected that the point-kinetics based normalization should already bring $w_{tot}/N$ close to one. This phase, like the previous one, avoids the normalization by a random variable (i.e. the simulation weight $w_{tot}$) and does not introduce any bias.

For the sake of completeness, it should be borne in mind that kinetic simulations starting from a critical state suffer from another source of possible bias, independently of how population control and importance sampling are implemented. Indeed, the initial neutron and precursor distributions at equilibrium are sampled through a power iteration, a method known to introduce a bias in the estimation of the average tallies, due to normalization [10, 11], and of their variance, due to the correlation of neutrons in adjacent cycles [13]. Additionally, it is often assumed that the starting point of each Monte Carlo replica of the time-dependent transient corresponds to the final fission neutron distribution of the previous replica: this choice considerably reduces the required computational cost (power iteration is performed only once), at the expense of breaking the independence of replicas [3]. This effect is mitigated by adding a few extra power iteration cycles between successive samplings of the kinetic source. For a more in-depth discussion of this statement, the reader is referred to the explanation of the kinetic simulations provided by Faucher et al. [3].

## 3 Benchmark specifications

The benchmark chosen for the analysis of the population-control algorithms described in Sec. 2 is Flattop-Pu (PU-MET-FAST-006) [7], which consists of a plutonium sphere surrounded by a natural uranium reflector, characterized by a fast spectrum. For the nuclear data, we use the ENDF/B-VII.r1 library.
<figure><img src="image_3.png" /></figure>

<table><thead><tr><th>Configuration</th><th>Critical</th><th>Supercritical</th></tr></thead><tbody><tr><td><i>k</i><sub>eff</sub></td><td>1.00037 ± 4×10<sup>-5</sup></td><td>1.00656 ± 4×10<sup>-5</sup></td></tr><tr><td>Λ<sub>eff</sub> [ns]</td><td>13.285 ± 0.0008</td><td>13.255 ± 0.0008</td></tr><tr><td>β<sub>eff</sub> [pcm]</td><td>278.95 ± 0.98</td><td>278.37 ± 0.99</td></tr><tr><td>β<sub>eff,1</sub> [pcm]</td><td>8.05 ± 0.17</td><td>7.82 ± 0.16</td></tr><tr><td>β<sub>eff,2</sub> [pcm]</td><td>52.39 ± 0.43</td><td>52.52 ± 0.42</td></tr><tr><td>β<sub>eff,3</sub> [pcm]</td><td>44.11 ± 0.39</td><td>44.55 ± 0.39</td></tr><tr><td>β<sub>eff,4</sub> [pcm]</td><td>97.34 ± 0.58</td><td>97.82 ± 0.57</td></tr><tr><td>β<sub>eff,5</sub> [pcm]</td><td>57.39 ± 0.44</td><td>55.53 ± 0.44</td></tr><tr><td>β<sub>eff,6</sub> [pcm]</td><td>19.67 ± 0.26</td><td>2.012 ± 0.26</td></tr><tr><td>ρ [$]</td><td>~ 0</td><td>~ 2.2</td></tr></tbody></table>

**Figure 1 & Table 1:** Left: Graphical representation of the Flattop-Pu benchmark in its critical configuration. The transient is initiated by enlarging the radius $R_1$ of the inner Pu sphere. Right: Table reporting the adjoint-weighted kinetics parameters, both for the critical and supercritical configuration. All the parameters have been obtained by running a TRIPOLI-4® simulation with the Iterated Fission Probability (IFP) method [12] using the ENDF/B-VII.r1 nuclear data library.

The choice of Flattop-Pu is known to provide an extremely challenging scenario for kinetics [3]: as the neutron time-scale decreases, the ratio $β/(λΛ)$ increases, resulting in a lower ratio between the neutron and precursor population sizes. Therefore, in this configuration, achieving satisfactory variance reduction is a highly demanding task for Monte Carlo. The system is initially kept in its initial, critical state for 1 µs. From Tab. 1 it is evident that the configuration identified as critical does not exhibit a multiplication factor precisely equal to one, due to the effects of the nuclear data library. In kinetic simulations, even a slight deviation from $k_{eff} = 1$ causes the system to drift away from the stationary state. To minimize this effect, we rescale the average fission multiplicity $v$ by the estimated multiplication factor $k_{eff}$ of the critical configuration. This adjustment allows mitigating the departure from criticality and practically enforces a steady-state configuration [3]. After 1 µs, the transient is initiated by suddenly enlarging the sphere of Pu. The inner radius is increased from $R_1 = 4.5332$ cm to $R_1 = 4.5650$ cm; the outer radius is kept unchanged and it is equal to 24.1420 cm. This injects a positive reactivity step of roughly $2.2$. After 25 µs, the inner radius is reduced to its original value, bringing the system back to criticality. The evolution of the transient is followed up to 40 µs.

## 4 Simulation results

In this section, we present the simulation results obtained for the transient described in Sec. 3. The fiducial quantity chosen for the benchmark is the energy-integrated neutron flux, normalized to the stationary value (i.e. the neutron flux in the critical system), integrated over the whole geometry. The results of the two algorithms are compared against each other and also against the analytical solution of the point-kinetics equations, considering only one family of precursors. The scores are collected on the time-steps of a time grid with a regular spacing of 0.5 µs, matching the reference grid used for population control. For both algorithms, we fixed the value of importance ratio at $R = 10^7$, basing this choice on previous analyses.

We explored different combinations of simulation parameters, tuning the number of particles per replica (i.e. *N*) and the number of replicas (i.e. *M*). The parameters were adjusted so that the total number of histories remained constant. The probed configurations were: 1500 particles with 1000 replicas, 3000 particles with 500 replicas, and 6000 particles with 250 replicas. This was done for the purpose of emphasizing the presence of a bias in the tallies
<figure><img src="image_4.png" /><figcaption>**Figure 2:** Comparison of the normalized neutron flux collected in Flattop-Pu during the transient described in section 3. The two algorithms are applied for three different combinations of simulation parameters. The trend is compared to the analytical solution of the point-kinetics equations, which has been solved for $\rho$, $\rho + \delta\rho$, and $\rho - \delta\rho$ to provide a window of uncertainty based on the precision of the reactivity parameter. Error bars represent one standard error.</figcaption></figure>

estimated using the legacy algorithm, and its absence when using the new algorithm. Based on previous investigations, the bias induced by normalizing by a stochastic quantity is known to be (at least to first order) proportional to the variance associated to the factor itself [11]. Therefore, we expect that the bias will exhibit a behavior proportional to $1/N$, decreasing as the number of neutrons per replica increases.

These considerations are supported by the findings illustrated in Fig. 2. The legacy algorithm gives reliable and acceptable results, in good agreement with the point-kinetics solution,
only when the number of simulated particles is high, which is clearly a consequence of the bias introduced being a decreasing function of the number of particles. Indeed, for the configuration with the lowest number of simulated particles (i.e. 1500), the departure from the point-kinetics model becomes the most pronounced, particularly after 20 µs. Furthermore, the legacy algorithm exhibits notably higher uncertainties in the scores, which also vary significantly across the three configurations probed, suggesting substantial fluctuations in the variance estimation. On the contrary, the new algorithm exhibits a more consistent pattern, showing remarkably similar results and less variation in the magnitude of the uncertainties across the three configurations tested. The stability of the variance trend can be directly attributed to the enhanced reliability of the new results with respect to those returned by the legacy algorithm. A more detailed analysis of the behavior of the uncertainty (in % units) is shown in Fig. 3.

<figure><img src="image_5.png" /><figcaption>**Figure 3:** Standard error on the neutron flux score obtained with the application of the legacy algorithm (left) and the new algorithm (right) for the three combinations of simulation parameters. In the first case, the standard error is noticeably higher and shows significant variation among the three configurations tested, while in the second case a more stable trend, with substantially lower standard error magnitudes, is observed across the three cases studied.</figcaption></figure>

Figure 3 prompts two crucial observations: first, it confirms that, when the new algorithm is used, the uncertainty on the scores is very stable and it does not fluctuate once the total number of simulated histories is fixed, thus affirming the reliability of the results. Actually, a small improvement of around one percentile point is observed for the configuration with the largest number of particles per replica. This improvement is consistent and may potentially stem from the reduction of the bias on the variance of the scores caused by the power iteration method, which is also known to be proportional to $1/N$ to first order [10]. Second, the results indicate significantly improved statistical performance in the refined algorithm, characterized by uncertainties on the scores consistently below 8%. In contrast, the legacy algorithm exhibits notably higher uncertainties, with values reaching up to 50% for some time-steps. For the legacy algorithm, the worst uncertainties are generally observed for the simulation with the smallest number of particles per replica. However, establishing a clear trend becomes challenging, especially after the system is brought back to criticality at $t = 25$ µs, since the statistics on the scores is strongly affected by the reduction of the neutron weights (recall that in the legacy algorithm the unbalance between the neutron weights and the precursor weights is addressed only at the beginning of the simulation, through the importance values $I_n$ and $I_c$). It is worth specifying at this point that the significantly improved statistics of the new algorithm is achieved without witnessing a radical increase in simulation time. Preliminary tests
have demonstrated a substantial FOM gain for the new algorithm. However, further work is necessary to provide a thorough assessment of this improvement.

<figure><img src="image_6.png" /><figcaption>**Figure 4:** Deviation [%] between the flux returned by the analytical solution of the point-kinetics equations and the flux obtained with TRIPOLI-4® using the legacy algorithm (left) and the new algorithm (right) for the three combinations of simulation parameters tested. Error bands represent one standard error.</figcaption></figure>

Figure 4 provides the percentage deviation from the point-kinetics result for the three combinations of parameters that we have tested. The deviation is estimated by normalizing the relative difference between the results of TRIPOLI-4® and the results of the point kinetics. The new algorithm shows again a better behavior, demonstrating a very good and similar agreement with point kinetics for all the cases, with a peak deviation of approximately 9%. This discrepancy can be partially attributed to the uncertainties associated to the Monte Carlo results, and partially to the fact that the point-kinetics model is by definition only an approximation of the exact kinetics. Further investigations are needed to discriminate the specific contributions of these factors. Worse results are observed for the legacy algorithm, where the deviations between TRIPOLI-4® and the point-kinetics model, in addition to displaying a very irregular behavior, can reach absolute values up to 80%. This deviation tends to decrease as the number of particles simulated increases, mainly due to the reduction of the bias.

## 5 Conclusions

In this work we have focused on population-control techniques for kinetic Monte Carlo simulations. We have described two algorithms: the legacy method currently implemented in TRIPOLI-4®, and an improved, point-kinetics-based version that is under testing. We have assessed both methods on the Flattop-Pu benchmark by simulating a step reactivity insertion of $2.2$. The goal of our investigation was to compare the performances of these algorithms in handling both the population growth and the unbalance between the neutron and the precursor populations. The fiducial quantity for the benchmark was the normalized neutron flux, taking the solution of the point kinetics model as a reference. While in the legacy algorithm the presence of a bias makes the simulation results unreliable if the number of simulated particles is low, the new algorithm shows good statistical results for all the sets of tested parameters. Furthermore, the new algorithm overall provides a smaller uncertainty on the scores. Additionally, it exhibits a better agreement with the point-kinetics model, in contrast to the legacy algorithm, which yields much larger deviations. In conclusion, the new
algorithm appears significantly superior in controlling the particle populations compared to the legacy one, paving the way towards more challenging scenarios such as severe reactivity injections or shutdowns. The overall statistical improvement, combined with the possibility of a FOM gain, also opens possibilities for simulating longer transients or transients in larger systems. Future work will aim at enhancing the newly developed algorithm by investigating the impact of updating kinetics parameters at each time step rather than computing them only once at critical conditions. Additionally, other implementations of point kinetics, particularly the alpha-eigenvalue formulation (where static reactivity is replaced by dynamic reactivity [14]) will be also explored.

## Acknowledgements

TRIPOLI-4® is a registered trademark of CEA. The authors thank EDF and IRSN for partial financial support.

## References

[1] B. L. Sjenitzer, J. E. Hoogenboom, Nucl. Sci. Eng. **175**, 94–107 (2013).

[2] E. Brun et al., Annals of Nuclear Energy **82**, 151-160 (2015).

[3] M. Faucher, D. Mancusi, A. Zoia, Ann. Nucl. Energy **120**, 74-88 (2018).

[4] D. Ferraro et al., in Proc. M&C2019, Portland, USA (2019).

[5] D. Mancusi, M. Faucher, A. Zoia, Eur. Phys. J. Plus **137** (1), 127 (2022).

[6] M. Faucher, M. "Coupling between Monte Carlo Neutron Transport and Thermal-Hydraulics for the Simulation of Transients Due to Reactivity Insertions", PhD thesis, Université Paris-Saclay, France (2019).

[7] OECD/NEA, "International Handbook of Evaluated Criticality Safety Benchmark Experiments" (1995).

[8] D. Mancusi, "Reactors, reactions, and random numbers", Habilitation à Diriger des Recherches, Université Paris-Saclay, France (2023).

[9] T. Booth, in Proc. RPSD Topical Meeting, Falmouth, USA (1996).

[10] R. Brissenden, A. Garlick, Ann. Nucl. Energy **13** (2), 63–83 (1986).

[11] E. M. Gelbard, R. E. Prael, Proc. NEACRP Meeting of a Monte Carlo Study Group (ANL-75-2), Argonne National Laboratory (1974).

[12] G. Truchet et al., Ann. Nucl. Energy **85**, 17–26 (2015).

[13] E. Dumontel et al., Ann. Nucl. Energy **63**, 612–618 (2014).

[14] N. Terranova, A. Zoia, Ann. Nucl. Energy **108**, 57-66 (2017).