# Adaptive adjoint-based population-control methods for kinetic simulations in TRIPOLI-4

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

► To cite this version:

Cecilia Montecchio, Davide Mancusi, Vincent Lamirand, Wilfried Monange, Andrea Zoia. Adaptive adjoint-based population-control methods for kinetic simulations in TRIPOLI-4. EPJ N - Nuclear Sciences & Technologies, 2025, 11, pp.74. 10.1051/epjn/2025069. hal-05442429

HAL Id: hal-05442429

https://asnr.hal.science/hal-05442429v1

Submitted on 5 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
REGULAR ARTICLE

OPEN ACCESS

# Adaptive adjoint-based population-control methods for kinetic simulations in TRIPOLI-4

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

Received: 20 May 2025 / Received in final form: 15 September 2025 / Accepted: 29 September 2025

**Abstract.** Time-dependent Monte Carlo simulations for reactor kinetics applications require special variance-reduction and population-control techniques in order to efficiently cope with the typically huge imbalance between the respective time scales and population sizes of neutrons and precursors. Building upon the legacy implementation of the algorithms devoted to kinetics in the Monte Carlo code TRIPOLI-4, in this work we propose an adaptive adjoint-based population-control method that considerably improves the behaviour of time-dependent simulations. Thanks to a time-dependent importance-sampling scheme, based on the solution of the adjoint point-kinetics equations, neutron and precursor weights are continuously adjusted, which paves the way towards the simulation of previously unattainable reactor transients involving long times and large reactivity excursions. The computational effectiveness of the newly developed method is evaluated in terms of Figure of Merit (FoM) over a set of time-dependent scenarios encompassing the Flattop-Pu, SPERT III E-core and CROCUS benchmarks.

## 1 Introduction

Monte Carlo methods are traditionally considered the gold standard for stationary neutron-transport problems in reactor physics applications. In recent years, they have been also successfully extended to reactor kinetics (including delayed neutron precursors contributions), a domain previously handled exclusively by deterministic solvers [1–9]. This improvement has been made possible thanks to the increasing availability of computational resources, combined with advancements in variance-reduction techniques [10,11], prompted by the pioneering work by Sjenitzer and Hoogenboom [12–17]. The development of efficient Monte Carlo methods for time-dependent problems is essential for the numerical validation of deterministic kinetic solvers, which rely on the discretization of the phase space and thus involve model-dependent approximations; in this respect, kinetic simulation capabilities have been implemented in the TRIPOLI-4® Monte Carlo code [18,19], developed at CEA, including several variance-reduction and population-control methods [3,10], and extensive benchmark comparisons have been carried out [2,4,20].

Despite these advances, kinetic Monte Carlo methods are still in their infancy and face some serious limitations in terms of Figure of Merit (FoM), since they typically

require extremely long computer times in order to achieve acceptable statistical accuracy on the sought response functions. These considerations severely hinder the simulation of transients involving large reactivity insertions or withdrawals. One of the toughest challenges in time-dependent Monte Carlo simulations arises from the imbalance between neutron and precursor population sizes, which stems from their largely different respective time-scales [13]. Indeed, the characteristic time-scale of the precursor population (i.e. the average decay time $1/\bar{\lambda}$) is much longer than the typical neutron time-scale (i.e. the mean generation time $\bar{\Lambda}$). This is mirrored in the ratio $\beta/(\bar{\lambda}\Lambda)$ between the two population sizes at critical conditions: for a typical PWR, a snapshot of the core at equilibrium would show $10^4$ more precursors than neutrons; the ratio would exceed $10^5$ for a typical fast reactor [13]. This substantial difference, possibly coupled to the unbounded growth or disappearance of the particles during a transient, imposes significant constraints on kinetic Monte Carlo simulations in terms of memory burden and computer time.

Extreme care is therefore required to manage both the population size and its evolution while addressing the neutron-precursor disparity. This is typically accomplished using population-control methods, which keep the overall particle population within reasonable bounds while optimizing the allocation of computational

<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), 1015 Lausanne, Switzerland</footnote>

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

<footnote>* e-mail: cecilia.montecchio@cea.fr</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Sketch of the current implementation of the kinetic algorithm in TRIPOLI-4, taken from Ref. [3]. In the kinetic simulation (green) neutrons and precursors are sampled during the last iteration of the power iteration (white) and their subsequent time evolution is followed on a time grid, during which forced precursor decay and population control algorithms are applied.</figcaption></figure>

resources [10]. Building on previous investigations concerning the use of the combing strategy [21], we have recently implemented a preliminary version of an adaptive population control algorithm in TRIPOLI-4 [22]. In this work, we considerably extend our findings and show how combing can be made highly efficient by using an adjoint-based importance-sampling scheme: a thorough numerical assessment of the importance-dependent strategy on a collection of benchmark configurations will substantiate our claims.

This manuscript is organized as follows: in Section 2 we briefly recall the key features of the legacy implementation of the kinetic algorithm in TRIPOLI-4, focusing primarily on the population control methods. In Section 3 we introduce an importance-dependent adaptive population control strategy, which considerably enhances the effectiveness of the method proposed in Ref. [22]. Then, in Section 4 we examine the FoM of the legacy and importance-dependent algorithms relying on the Flattop-Pu [23] and on the SPERT III E-core benchmark [24]. Further analyses are carried out in order to probe the performance of the population-control methods for scenarios corresponding to large reactivity and long transients, relying on the Flattop-Pu and on the CROCUS benchmark [25]. Conclusions are finally drawn in Section 5.

## 2 Overview of the kinetic capabilities of TRIPOLI-4

We begin by briefly recalling the key features of the legacy algorithm currently implemented in TRIPOLI-4. For an exhaustive description of the kinetic capabilities of the code, the reader is referred to the original work by Faucher et al. [3].

In TRIPOLI-4, the kinetic simulation scheme is divided into two main stages, as illustrated in Figure 1. The first stage is used to sample the critical distribution of neutrons and precursors, if needed. It consists of a power iteration over a number of cycles, called inactive cycles, because they do not contribute to the final tallies. Their

number must be large enough to ensure convergence of the fission source. The power iteration is performed on the equilibrium configuration of the simulated system; once convergence is attained, the neutron and precursor populations are sampled within a further power iteration cycle. These populations then become the initial condition for the second stage, where the particles are simulated in time-dependent conditions, over a user-defined time mesh used for tallying, and the geometry and the material compositions of the system are updated according to the chosen kinetic scenario. This two-stage scheme is repeated for a sufficiently large number of batches, and the final result reported by TRIPOLI-4 is obtained by averaging over them. In principle, each batch should be an independent calculation including both a source-convergence and a kinetic phase, but this is computationally demanding. To reduce this cost, the fission neutrons produced in the last inactive cycle of one batch are reused to initialize the next. This approach reduces overhead but introduces correlations between batches, which are mitigated by adding a few additional de-correlation inactive cycles. In parallel runs, each task executes the scheme shown in Figure 1.

Due to the combined effect of nuclear data libraries and technological data, the effective multiplication factor $k_{eff}$ of a simulated system may deviate from unity for a configuration that is experimentally observed to be critical. Correspondingly, the associated fundamental eigenmode $\phi_k$ will also deviate with respect to the actual stationary flux. When performing a kinetic simulation starting from equilibrium conditions, it is therefore very important to correct this bias, otherwise we would observe a departure from the stationary condition for the neutron and precursor population of the supposedly critical system. For this purpose, in TRIPOLI-4 we normalize the average number of neutrons generated by fission by $k_{eff}$, forcing the system to be critical.

### 2.1 Non-analog precursor sampling

Precursors are treated as a distinct particle. Currently, precursors retain the spatial coordinates of the neutrons from which they are created, i.e. their migration is neglected. The large difference in the time-scales of neutrons and precursors, which would lead to intractable issues in a naive time-dependent simulation, is taken into account by suitably altering the sampling laws of precursors, and correspondingly adjusting their statistical weights in order to preserve the unbiasedness of the tallies, based on a strategy proposed by Legrady and Hoogenboom [12]. Precursors belonging to different families are grouped together in a single representative precursor particle, which reduces the variance due to the dispersion of the decay constants $\lambda_j$. Furthermore, the decay time distribution of such representative precursor is forced to be uniform within every simulated time step $[t_i, t_{i+1}]$, namely,

$$p_{\text{decay}}(t) = \frac{1}{t_{i+1} - t_i} \chi(t_i, t_{i+1}), \quad (1)$$

where $\chi(a, b)$ is the marker function of the interval $[a, b]$, and $i$ is the index of the time step. This enhances the emissions
sion of a delayed neutron in each time interval, compared to the analog exponential distribution of the decay times. To preserve an unbiased Monte Carlo game, the statistical weight of the delayed neutrons emitted upon precursor decay must be modified accordingly, i.e.

$$W_{\text{delayed}}(t|t_0) = (t_{i+1} - t_i) \sum_j \frac{\lambda_j \beta_j}{\beta} e^{-\lambda_j(t-t_0)} w_C, \quad (2)$$

where $w_C$ is the initial weight of the emitting precursor, i.e. the weight of the neutron from which the precursor was generated, $\beta_j$ is the fraction associated with the delayed neutron precursor family j (with $\beta = \sum_j \beta_j$), and $t_0$ is the precursor creation time. Upon decay, the precursor is not killed but instead added to a buffer of particles to be followed in the next time step $i+1$, with a weight equal to:

$$W(t_{i+1}|t_0) = w_C \sum_j \frac{\beta_j}{\beta} e^{-\lambda_j(t_{i+1}-t_0)}. \quad (3)$$

## 2.2 Population control algorithms

If the particles are left free to evolve during a kinetic transient, the population size over time would either grow unbounded if the system is supercritical, or shrink to zero if the system is subcritical, both conditions being unacceptable for the Monte Carlo simulation. Even a critical system, where the population stays constant on average, would lead in practice to intractable conditions for the simulation, since the variance would increase linearly in time because of the fluctuations induced by the fission chains [26]. Furthermore, the significant disparity between neutron and precursor time-scales creates an imbalance among the two populations, favoring precursors over neutrons; since tallies in the simulation are carried by neutrons, this situation is again detrimental to the overall variance of the sought responses.

The set of approaches designed to tackle these issues is collectively referred to as population control strategies. The goal of such techniques is to maintain a suitable number of particles to be simulated while also addressing the pronounced disparity between neutrons and precursors, thus ensuring acceptable statistics. In TRIPOLI-4, population control strategies are applied on the same time grid as forced decay, minimizing the need to define new grids. For convenience, we can distinguish two main procedures in a population control algorithm:

1. The scaling of the neutron and precursor populations according to an importance-sampling scheme;

2. The application of either combing to both neutrons and precursors, or Russian roulette specifically to precursors to control the actual population size at the end of each time step.

The role of the first procedure is to promote the allocation of computer resources for the neutron population in the transport simulation, in view of the previous considerations about tallies being sensitive to neutrons only. In this work we have extensively modified this procedure

with respect to the legacy implementation in TRIPOLI-4 described in Ref. [3]: a thorough description of the novel implementation is provided in Section 3.

The role of the second procedure is to enforce the control of the population sizes, in order to restore a convenient amount of particles for the subsequent time step, based on the particles that survived until the end of the current time step.

If combing is chosen, the number of particles is set to a target value $M$, while their total weight and the proportion of respective weights (in the case of multiple particle species) are preserved. The new $M$ particles are sampled among the $K$ current particles by using two equally spaced weight intervals, one for neutrons and one for precursors, over the total available weight. It has been shown that this algorithm is unbiased, provided that a random offset is imposed when selecting the first particle [21].

If Russian roulette for precursors is chosen, then the population size of precursors is controlled through the comparison with a reference weight: if the weight of the particle is lower than a predefined threshold, it is subjected to a random survival test. The roulette is applied on the basis of the expected delayed neutron weight $w_{\text{delayed}}$ that can be emitted during the following time step $[t_i, t_{i+1}]$, which follows from [27]

$$w_{\text{delayed}}(t_i, t_{i+1} | t_0) = w_C \sum_j \frac{\beta_j}{\beta} \left( e^{-\lambda_j(t_i - t_0)} - e^{-\lambda_j(t_{i+1} - t_0)} \right). \quad (4)$$

This quantity is compared to a threshold value, which in TRIPOLI-4 is set by default to 0.8. A detailed comparison of the performance of combing versus Russian roulette for population control will be provided in Section 4.

In addition to either combing or Russian roulette at the end of each time step, in TRIPOLI-4 Russian roulette and splitting are also applied to neutrons within each time step after every collision, as a standard variance-reduction technique.

## 2.3 Legacy importance-sampling scheme for neutrons and precursors

The kinetic algorithm of TRIPOLI-4 features an importance-sampling scheme, whose legacy implementation was described in Refs. [3] and [10]. In a recent work, we have introduced a preliminary version of an enhanced importance-sampling scheme for kinetics calculations [22], whose comprehensive description is provided in Section 3. The key motivation behind such importance-sampling scheme is to adjust the weights of particles and tame the imbalance between neutron and precursors, prioritizing neutrons in memory allocation. For the sake of clarity, we stress that the notation used in this paper differs from the one adopted in Faucher et al. [3,10]: we use $1/I_n$ instead of $I_n$, and $1/I_c$ instead of $I_c$, respectively. The importance ratio $R = I_n/I_c$ is still defined as the ratio between the neutron importance and the precursor importance, but the values of the absolute importance are
altered. The reason for this change, which was suggested in Ref. [28], is that the new definition yields importance ratios significantly larger than one for short time steps, which aligns better with the intuitive understanding that neutrons play a more significant role than precursors on short time-scales. Moreover, this definition is more consistent with the interpretation of importance as the expected contribution to the detector response and as the solution to the adjoint Boltzmann equation.

We have then to distinguish physical and simulation weights. The statistical weights assigned to neutrons and precursors in the simulation, denoted respectively as $w_n$ and $w_c$, are referred to as physical weights. However, to optimize the computational handling of neutron and precursor populations, we replace them with the simulation weights, $\tilde{w}_n$ and $\tilde{w}_c$. In the legacy algorithm, the simulation weights $\tilde{w}_n$ and $\tilde{w}_c$ are derived from the physical weights $w_n$ and $w_c$ through the following relationships:

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

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

More specifically, the legacy algorithm is structured so that $I_c = 1$ and $I_n = R = I_n/I_c$, meaning that only the neutron weights are effectively corrected, while the precursor weights remain unmodified. The sum of the physical weights is required to be equal to the sum of the simulation weights, ensuring the conservation law:

$$ w_n + w_c = \tilde{w}_n + \tilde{w}_c. \qquad (6) $$

Using the constraint in equation (6) and referring to the definitions of $\tilde{w}_n$ and $\tilde{w}_c$ in equation (5), we can derive the explicit expressions of the neutron importance $I_n$ and of the precursor importance $I_c$ as a function of the importance ratio $R$, namely:

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

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

$$ R = \frac{I_n}{I_c}. \qquad (7c) $$

The subscript tot denotes the sum of all particle weights. The importance ratio $R$ represents a degree of freedom in this importance-sampling scheme. In the legacy version of TRIPOLI-4, $R$ is a user-defined parameter. When a list of importance ratios is specified, along with the corresponding time values at which changes must be applied, the algorithm adjusts the weights at the beginning of each time step by satisfying two conditions: (1) the equivalence of the sum of the physical weights from time step $i$ to the subsequent time step $i+1$ and (2) the preservation of the physical weights for both neutrons and precursors, namely

$$ w_{n,i} + w_{c,i} = w_{n,i+1} + w_{c,i+1} \qquad (8a) $$

$$ \frac{\tilde{w}_{n,i+1}}{I_{n,i+1}} = \frac{\tilde{w}_{n,i}}{I_{n,i}} \qquad (8b) $$

$$ \frac{\tilde{w}_{c,i+1}}{I_{c,i+1}} = \frac{\tilde{w}_{c,i}}{I_{c,i}}. \qquad (8c) $$

In the legacy algorithm, the transformation of physical weights into simulation weights are carried out at the end of the inactive cycles of the power iteration, for the particles that have been sampled during the last iteration of the critical phase and represent the initial condition for the kinetic simulation [3].

During the simulation of the time steps, each precursor created by a neutron is assigned a weight equal to the neutron simulation weight divided by the importance ratio, namely $\tilde{w}_{c,i} = \tilde{w}_{n,i} \cdot I_c/I_n$. Conversely, the weight of a precursor decaying into a neutron is multiplied by the importance ratio, namely $\tilde{w}_{n,i} = \tilde{w}_{c,i} \cdot I_n/I_c$. If the importance ratio changes from one step to the next, the weights of neutrons and precursors are adjusted according to equations (8b) and (8c), respectively. Finally, during the scoring phase, any tally involving neutron weights is divided by the neutron importance to ensure unbiased results. A graphical representation of this scheme, considering two time steps for which two different importance ratios are defined, is shown in Figure 2. If only one scalar value is provided for $R$ across the whole time grid, the conversion from the physical weights to the simulation weights (as per Eq. (5)) is performed only once, at the beginning of the kinetic simulation, based on the particle weights inherited from the critical source (Eqs. (7a) and (7b)). Selecting a global value for $R$ for the entire transient might be a poor choice, especially for cases evolving over long time-scales and/or with large reactivity excursions, where the weights of the particles can change significantly. As a result, this approach is not particularly effective to deal with such scenarios.

Despite its apparent simplicity and its easy application, the legacy implementation of the importance-sampling scheme suffers from two major drawbacks.

First, a comprehensive investigation has shown that there exists an optimal value of the importance ratio $R$, which maximises the simulation performance. However, this value is system-dependent and varies with the reactivity introduced into the system and with the simulation time steps used for applying the population control algorithm [10,20]. A non-homogeneous time mesh with varying time step sizes might be extremely convenient in many practical problems. Indeed, in the simulation of a critical system a very fine grid for the population control is not needed, since the population is expected to remain stable. On the contrary, when the reactivity excursion is very fast (e.g. during a prompt jump or drop) one might be interested in controlling the population more often; in such case, finer grids might be necessary. Thus, determining the optimal value of $R$ for each time step, while also accounting for reactivity changes, is a challenging task that often relies on a trial-and-error approach [3]. Optimizing algorithm performance would require extensive verification to establish the most suitable set of importance ratios across the entire simulation domain.

Second, more subtly, the population importance values $I_n$ and $I_c$ used in the algorithm are stochastic quantities. In equation (5), the weights are multiplied by a stochastic value, itself depending on the particle weights; thus, equation (5) is a non-linear operation, and has been pointed out to possibly lead to a bias [22], similarly to what happens
<figure><img src="image_3.png" /><figcaption>Fig. 2. Schematic representation of the importance sampling strategy implemented as a variance reduction technique in the kinetic algorithm in TRIPOLI-4. The scheme is applied to a kinetic simulation with two time steps, considering two different values $R_i$ and $R_{i+1}$ for the importance ratio, one for each time step.</figcaption></figure>

during the normalization procedure used in the regular power iteration [29,30]. In particular, the expected value $E[f(\tilde{w}_n)]$ of any estimator involving a nonlinear function $f$ will be different from the true value $f(E[\tilde{w}_n])$.

### 3 A refined adjoint-based importance-sampling scheme

In order to address the aforementioned issues, we have implemented an improved version of the importance-sampling scheme. The conversion from physical weights to simulation weights is performed according to equation (5), as in the legacy version. The major change in the algorithm concerns the elements appearing in the definition of the importance functions, shown in equation (7), which are required for the weight conversion. Instead of estimating these quantities based on the neutron and precursor weights obtained along the Monte Carlo simulation, we rely on the system of deterministic one-family point-kinetics equations

$$ \frac{dn(t)}{dt} = \frac{\rho_k(t) - \beta_{\text{eff}}}{\Lambda_{\text{eff}}} \cdot n(t) + \bar{\lambda} \cdot c(t) \quad (9a) $$

$$ \frac{dc(t)}{dt} = \frac{\beta_{\text{eff}}}{\Lambda_{\text{eff}}} \cdot n(t) - \bar{\lambda} \cdot c(t), \quad (9b) $$

where $\Lambda_{\text{eff}}$ is the effective mean generation time, $\beta_{\text{eff}}$ the effective delayed neutron fraction, $\bar{\lambda}$ the average precursor decay time and $\rho_k$ the static reactivity. The parameters $\Lambda_{\text{eff}}$, $\beta_{\text{eff}}$, and $\bar{\lambda}$, which are required in equations (9), are estimated by TRIPOLI-4 through an eigenvalue calculation using the IFP method [31]. This calculation must be performed prior to the kinetic simulation. It is implicitly assumed that these values are only weakly affected by the system evolution during the transient, so that it is safe to consider them to be time-independent. On the contrary, the system reactivity is estimated using the definition of static reactivity $\rho_k = 1/k_{\text{initial}} - 1/k_{\text{perturbed}}$, running a series of eigenvalue calculations for each time-step configuration where the reactivity changes. It is worth noting

that reactivity values can be obtained from power iteration calculations, which are much less computationally expensive than the target kinetic simulation<fnref n="1" />. The analytical solution of the point-kinetics equations (9) is provided in the Appendix. The decision to collapse all precursor families into a single effective family is driven by the need of a tractable analytical solution. The investigation of the impact of considering the full system of equations with multiple precursor families is left for future work.

Based on the solution of equations (9), we can determine the importance values for each time step, namely

$$ I_n(t_i) = R(t_i) \cdot \frac{n(t_i) + c(t_i)}{R(t_i) \cdot n(t_i) + c(t_i)} \quad (10a) $$

$$ I_c(t_i) = \frac{n(t_i) + c(t_i)}{R(t_i) \cdot n(t_i) + c(t_i)}. \quad (10b) $$

These expressions are equivalent to equations (7a) and (7b), the only difference being that the physical weights of the particles are replaced by the deterministic estimates $n(t_i)$ and $c(t_i)$, which are evaluated at the end of each time step. While equations (9) provide a highly simplified and approximate representation of the kinetic response of the system, we will later show that they are sufficient in this context to capture the overall time evolution and contribute to an effective correction of particles weights.

In the improved version of the importance-sampling scheme, the importance values obtained from equations (10) are non-stochastic, which suppresses the emergence of the aforementioned bias due to the non-linear operation involving statistical weight multiplications. This is a remarkable advantage over the legacy implementation. Furthermore, the proposed approach naturally takes into account the evolution of the importance values over time, which allows for a built-in dynamic adjustment of physical weights to simulation weights, a feature which is not achievable with the legacy algorithm, since the definition

<footnote n="1">¹ However, additional challenges may arise when coupling TRIPOLI-4 with an external thermal-hydraulics solver.</footnote>
of the importance values given in equations (7a) and (7b) remains constant over time and changes only if the user provides a different value of R.

Equations (10) depend on the importance ratio $R(t_i)$, which should be assigned for each time step. In the improved version of the importance-sampling scheme, we have implemented a feature that enables to estimate optimal values of $R(t_i)$ automatically, in order to alleviate the difficulties of the legacy algorithm mentioned in Section 2.3. The key idea is to leverage the equivalence between the importance function and the solution of the adjoint time-dependent transport equations [32]: the adjoint neutron density $n^\dagger$ and the adjoint precursor density $c^\dagger$ physically represent the importance of a neutron (respectively, precursor) with respect to a detector. We define the detector response as the integral over the time interval $[t_i, t_{i+1}]$ of the product between the adjoint neutron source $Q^\dagger$, interpreted as the detector response function, and the neutron density $n$. Here the adjoint source $Q^\dagger$ is equal to 1 between $t_i$ and $t_{i+1}$, and zero elsewhere. It is assumed that precursors do not directly contribute to the detector response. The importance ratio can be thus derived by taking the ratio between the solutions of the adjoint transport equations, i.e., $R = n^\dagger/c^\dagger$. By analogy with the strategy used to infer the forward neutron and precursor densities $n$ and $c$, we assume that the quantities $n^\dagger$ and $c^\dagger$ are determined by solving the set of adjoint point-kinetics equations, namely

$$-\frac{dn^\dagger(t)}{dt} = \frac{\rho_k(t) - \beta_{\text{eff}}}{\Lambda_{\text{eff}}} n^\dagger(t) + \frac{\beta_{\text{eff}}}{\Lambda_{\text{eff}}} c^\dagger(t) + Q^\dagger(t) \quad (11a)$$

$$-\frac{dc^\dagger(t)}{dt} = \bar{\lambda}n^\dagger(t) - \bar{\lambda}c^\dagger(t). \quad (11b)$$

The exact solutions of equations (11) are recalled in the Appendix; for a thorough discussion on the role of the adjoint flux in zero-variance schemes, see e.g. Ref. [11]. The evolution of $n^\dagger$ and $c^\dagger$ with respect to time depends on the physical features of the simulated system (i.e., $\Lambda_{\text{eff}}$, $\beta_{\text{eff}}$ and $\bar{\lambda}$), on the reactivity $\rho_k$, and on the simulation parameters (i.e., the time step $\Delta t$ for applying the population control strategy). The solution of the system of equations (11) provides the values $n^\dagger(t_i)$ and $c^\dagger(t_i)$, whose ratio at each time step $i$ serves as an estimate of the importance ratio $R(t_i) = n^\dagger(t_i)/c^\dagger(t_i)$.

At the initial time, assuming a critical reactor condition, the precursor-to-neutron population ratio is given by $\beta_{\text{eff}}/(\Lambda_{\text{eff}}\bar{\lambda})$, which is a very large number for typical thermal and fast systems. Correspondingly, the ratio between the neutron and precursor importance is expected to be very large, reflecting the imbalance in the respective particle contributions to the response, i.e., the integral of the neutron density $n$ over $[t_i, t_{i+1}]$. When a reactivity change is induced in the system, the importance of neutrons and precursors adjusts accordingly, depending on the sign of the reactivity change. For a positive insertion, the net population of both neutrons and precursors increases, reducing the relative importance of a single particle in contributing to the response, since the overall larger population makes the individual contribution less meaningful. Conversely, for a negative reactivity insertion,

the population decreases, and the importance of individual particles increases as their contributions will have a larger impact. Furthermore, if the response involves the simulation of very small time steps, the importance of a neutron will be significantly higher than in simulations with larger time steps.

The time dependence of the importance ratio $R(t_i)$ computed according to this scheme ultimately reduces to a dependence on the size of the time step $[t_i, t_{i+1}]$. Indeed, as easily inferred from equation (A.2) in the Appendix, the solutions of equations (11) are evaluated at $t = t_i$, meaning that they depend only on the time-step size $t_{i+1} - t_i$. In the improved version of the importance-sampling scheme, it is still possible to override the importance ratio determined from solving the adjoint point-kinetics equations (11) and impose instead a constant user-provided value. However, even when a constant value of R is used, the new method differs from the legacy approach, as the importance values in the new algorithm depend on the neutron and precursor concentrations computed using the forward point-kinetics equations (10), whereas in the legacy algorithm they are based on the weights given in equation (7). Therefore, the improved importance-sampling is better tailored to follow the kinetic evolution of the simulated scenario and is expected to perform better than the legacy algorithm.

In conclusion, the main features of the new importance-sampling scheme calculation are:

(i) The definition of the importance values depends on the importance ratio and on estimates coming from the solution of the point-kinetics equations. This dependence is particularly advantageous as it enables a step-by-step tracking of the system evolution. Consequently, it provides a finely tailored update of the importance values, taking into account the changing sizes of the neutron and precursor populations;

(ii) The user is not required to carry out a trial-and-error procedure to determine the best importance ratio for the system: a set of values is automatically determined based on the solution of the adjoint point-kinetics equations;

(iii) The importance ratio depends on specific parameters of the system and varies according to changes in the time-step size or in the reactivity value.

## 4 Testing the new kinetic capabilities of TRIPOLI-4

When comparing the effectiveness of various Monte Carlo simulations, a commonly used metric is the figure of merit (FoM), namely

$$\text{FoM} = \frac{1}{\sigma^2 T}, \quad (12)$$

where $\sigma^2$ is the estimated variance for a given observable and $T$ is the duration of the simulation. In this section we investigate the population control strategy proposed in the previous section using the FoM as a comparative metric. Specifically, we evaluate the performance of the
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Graphical representation of the Flattop-Pu benchmark in its critical configuration.</figcaption></figure>

two strategies described in Section 2.2, which are combing and Russian roulette on precursors, when combined with the importance-sampling scheme, exploiting different values of the importance ratio. All FoM calculations presented hereafter are conducted using alternatively the legacy and the refined algorithm and they are executed as single-processor simulations, in order to avoid potential distortions associated with parallelism. We emphasize that our analysis is strictly limited to configurations that are within reach of the legacy algorithm; otherwise, a meaningful comparison would not be possible. All the simulation results presented in this section are based on the ENDF/B-VIII.0 nuclear data library [33].

## 4.1 A critical system: Flattop-Pu

We start by considering the PU-MET-FAST-006 (Flattop-Pu) benchmark [23], which has been previously used for the assessment of the preliminary implementation of the revised importance sampling algorithm in TRIPOLI-4 [22]. This benchmark consists of a plutonium sphere surrounded by a natural uranium reflector, and is characterized by a fast spectrum. For illustration, see Figure 3; the kinetics parameters of Flattop-Pu are provided in Table 1. Despite its apparent simplicity, the Flattop-Pu benchmark provides an extremely challenging scenario for kinetics [3], because of the very low value of the mean generation time, $\Lambda_{\text{eff}} \sim 10$ ns, which induces a huge imbalance between the neutron and precursor population sizes at equilibrium.

Our first FoM comparison is based on the following reference simulation: the kinetic evolution is performed over a single time step, during which the system is maintained in its critical state<fnref n="2" />. The observable chosen for the FoM analysis is the ensemble-averaged time integral of the energy and volume-integrated neutron flux. We run this simulation using the legacy version of TRIPOLI-4 and

<table><caption>**Table 1.** $k_{\text{eff}}$ and adjoint-weighted kinetics parameters for the nominally critical configuration. All the parameters have been obtained by running a TRIPOLI-4 simulation with the Iterated Fission Probability (IFP) method [31].</caption><thead><tr><th>Configuration</th><th>Critical</th></tr></thead><tbody><tr><td>$k_{\text{eff}}$</td><td>0.99856 ± 4 × 10<sup>-5</sup></td></tr><tr><td>$\Lambda_{\text{eff}}$ [ns]</td><td>13.123 ± 0.008</td></tr><tr><td>$\beta_{\text{eff}}$ [pcm]</td><td>279.12 ± 0.98</td></tr><tr><td>$\beta_{\text{eff},1}$ [pcm]</td><td>7.65 ± 0.16</td></tr><tr><td>$\beta_{\text{eff},2}$ [pcm]</td><td>52.23 ± 0.42</td></tr><tr><td>$\beta_{\text{eff},3}$ [pcm]</td><td>45.05 ± 0.40</td></tr><tr><td>$\beta_{\text{eff},4}$ [pcm]</td><td>98.15 ± 0.58</td></tr><tr><td>$\beta_{\text{eff},5}$ [pcm]</td><td>56.02 ± 0.44</td></tr><tr><td>$\beta_{\text{eff},6}$ [pcm]</td><td>20.02 ± 0.26</td></tr></tbody></table>

the refined algorithm. To evaluate the optimal population control strategy, we investigate the following schemes:

(i) Russian roulette combined with the importance-sampling scheme, exploiting one by one six different values of importance ratio, namely $R = \{1, 10^3, 10^6, 10^9, 10^{12}, 10^{15}\}$. For the new algorithm, we test also the value of $R$ computed solving the adjoint point-kinetics equations, as described in Section 3.

(ii) Combining combined with the importance-sampling scheme, exploiting one by one the same six values of importance ratio as listed in (i) and, for the new algorithm, the values computed solving the adjoint point-kinetics equations.

All these simulations are repeated for three different sizes of the simulated time step, namely $\Delta t = \{10^{-8}, 10^{-9}, 10^{-10}\}$ s. The values of $\Delta t$ have been selected based on the value of $\Lambda_{\text{eff}}$. Additionally, since this analysis encompasses simulations over a wide range of importance ratio values, using larger time steps would be impractical due to the increased computation time.

For each kinetic simulation, the number of particles and cycles is adjusted to ensure that the estimated standard error on the total neutron flux remains below 20%. When using lower importance ratio values, simulations require a larger number of histories to achieve the target accuracy, due to the very large imbalance between neutrons and precursors. In this kinetic scenario, involving only a single time step over which the system is critical, the differences between the legacy and the new algorithm are not expected to be significant for the same value of importance ratio. This occurs because here the only discrepancy between the importance-sampling schemes lies in how the importance values are computed at the initial time to obtain the *simulation* weights for the single time step: the legacy algorithm will use weight estimates from the power iteration, whereas the improved algorithm will use the point-kinetics equations. Consequently, this investigation of a critical system cannot fully showcase the improvements of the revised importance-sampling scheme in dynamically adjusting the importance values according

<footnote n="2">² The system is nominally critical, but the criticality simulation with ENDF/B-VIII.0 yields a reactivity of the order of $-145$ pcm ± 4 pcm, as reported in Table 1. This bias may be attributed to uncertainties in nuclear data or geometry, as discussed in Section 2.1. To address this issue, the simulation was performed by rescaling the value of $\nu$ by the value of $k_{\text{eff}}$, which makes the system critical by construction, up to the statistical uncertainty on the simulated value of the effective multiplication factor.</footnote>
<figure><img src="image_5.png" /></figure>

**Fig. 4.** Analysis of the FoM gain associated to different population control strategies and choices of the population importance ratio for the kinetic simulation of a single time step of the Flattop-Pu benchmark in its critical configuration. The comparison is repeated for three different choices of the time step size: $\Delta t = 10^{-8}$ s, $\Delta t = 10^{-9}$ s and $\Delta t = 10^{-10}$ s. For each time step size all the computed FoM values are normalized by the FoM value of a reference simulation where Russian roulette is employed and no importance-sampling scheme is applied ($R = 1$). The results of the legacy algorithm are displayed in red, while those of the new algorithm are shown in black. The circle in the curve of the new algorithm represents the value of FoM gain associated with the importance ratio based on the solution of the adjoint point-kinetics equations.

to the system evolution; nonetheless, the results illustrated in the following clearly demonstrate the impact of the chosen population control strategy on the FoM and confirm the existence of an optimal range of values for the importance ratio.

In Figure 4 we show the FoM gain associated with different population control strategies and importance ratios used in the importance-sampling scheme, for three different sizes of the time step. For each time step, the FoM gain is always measured relative to a reference simulation in which Russian roulette is employed, and the importance-sampling scheme is not applied ($R = 1$). Notably, the FoM gain shows a strong dependence on the selected value of importance ratio $R$, whereas no significant differences are observed between using combing over Russian roulette. The observed behavior of the FoM gain as a function of the importance ratio, characterized by a steep initial increase followed by a plateau, is consistent across all tested time step sizes and reflects the fact that the simulation is sub-optimal for low values of $R$. This is due to the disparity in the weight values between precursors and neutrons, leading to an under-sampling of neutrons. Precursors decay within the time step based on the forced decay scheme described in Section 2.1: without an appropriate importance adjustment, they will produce delayed neutrons with extremely low weights. Most of these neutrons are subsequently eliminated through the Russian roulette that is enforced after each collision in TRIPOLI-4. Therefore, choosing an importance ratio that increases the weight of neutrons are crucial and clearly improves the performance of the simulation by various orders of magnitude. The FoM gain associated with the importance ratio relying on the adjoint point-kinetics equations are located within the plateau region, indicating that the computed importance ratio provides an optimal FoM gain.

Previous studies have shown that the FoM gain also depends on the time step size [10,20]. Specifically, it has been observed that shorter time steps allow for a higher maximum gain achievable through the importance-sampling scheme, for both combing and Russian roulette. However, in the simulations conducted in this work, the observed differences in FoM gain are relatively small, suggesting that the previously identified trend is not strongly pronounced for our choices of time step sizes. Overall, as briefly outlined above, the computational efficiency trends of the legacy and improved implementations are rather consistent; the importance-sampling scheme has a significant impact on improving the FoM, and the value of the importance ratio derived from the solution of the adjoint point-kinetics equations is highly effective.

## 4.2 Analysis of a reactivity insertion in the Flattop-Pu benchmark

The next step of our analysis addresses more challenging scenarios for kinetic simulations, where the system transitions from a critical state to a supercritical or sub-critical state. For this purpose, we consider a positive reactivity insertion into the initially critical Flattop-Pu benchmark, inspired by previous investigations in Ref. [22]. After 0.5 µs we apply a step-wise modification of the volume of the inner plutonium sphere, from an unperturbed radius $R_1 = 4.5332$ cm to a perturbed radius $R_1 = 4.5650$ cm, by keeping its density fixed, which results in a reactivity insertion of approximately $2.2$. The radius of the external
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Left: evolution of the neutron flux in the Flattop-Pu benchmark after a positive reactivity insertion of $2.2$, according to the solution of the point-kinetics equations with one family of precursors. Middle: Evolution of the importance ratio for every time-step used for the kinetic simulation. The population control strategy is performed on a grid that becomes finer as the neutron flux increases. The values of $I_n$ and $I_c$ are calculated according to equations (A.2); consequently, they vary with the simulated time-step length and with the considered reactivity. Right: Evolution of the importance values ($I_n$ and $I_c$) for the same simulation.</figcaption></figure>

sphere remains constant. The inner sphere is then brought back to its original size after 2.5 µs, which restores criticality. Based on estimates derived from point kinetics, using the parameters shown in Table 1, we find that this transient causes the flux to increase to approximately 2.7 times its initial value.

We assess the FoM across six distinct time intervals of different length. In the first step ([0–0.5] µs), the system remains critical. During the second ([0.5–2] µs), third ([2–2.75] µs), and fourth ([2.75–3] µs) steps, the system is kept in the supercritical state described above. Finally, in the fifth ([3–3.25] µs) and last step ([3.25–4] µs), the critical configuration is restored. Time steps of different size have been chosen in order to increase the frequency of population control during the supercritical excursion. For the performance assessment, we estimate the FoM($t_i$) value for each time step $t_i$, according to equation (12). We rely on the same observable as in Section 4.1 (i.e. the time, energy and volume-integrated neutron flux) tal- lied for the specific time step. Regarding the simulation time, we consider the cumulative time required to per- form the dynamic simulation up to and including the time step for which the FoM is assessed. The simulations are performed alternatively with both importance-sampling schemes, combined with combing and Russian roulette. The importance ratio is either computed by solving the adjoint point-kinetics equations, or imposed by choosing values in $R = \{1, 10^3, 10^6, 10^9, 10^{12}\}$.

The point-kinetics evolution of the kinetic scenario is illustrated in Figure 5, using the parameters provided in Table 1. We also display the evolution of the importance ratio computed using the improved importance-sampling scheme described in Section 3, along with the evolution of the neutron and precursor importance. The importance ratio depending on the adjoint point-kinetics equations changes whenever there is a variation either in the time step size or the static reactivity value, which leads to an adaptive time-dependent control of the neutron popula- tion through the factor $I_n$. In Figure 6 we present the

evolution of the FoM gain as a function of the importance ratio across the six time steps analyzed, combined with either combing (Fig. 6i) or Russian roulette (Fig. 6ii). The improved algorithm is compared to the legacy implemen- tation. For each of the six time steps, the FoM gain is assessed with respect to a reference simulation using the importance ratio $R = 10^3$, combined with either comb- ing (Fig. 6i) or Russian roulette (Fig. 6ii). The case of $R = 1$ (i.e. no importance-sampling) is not used for our comparisons, since numerical tests have shown that for this configuration it is extremely difficult to achieve sta- tistical uncertainties below 20% for all time steps, due to the strong under-sampling of the neutron population.

This analysis suggests that variance-reduction tech- niques are not only beneficial for improving statistics but are in practice mandatory for obtaining meaning- ful average results in the first place. When combing is employed, the trend across the time steps shows a decrease in the maximum achievable FoM gain. In contrast, with Russian roulette, the FoM gain remains relatively con- stant across time steps, except for the very first one, where it reaches a slightly higher value. In both cases, the use of the importance-sampling scheme is crucial for optimiz- ing the simulation performance. In general, after a sharp increase, a plateau of the FoM is observed, extending over a wide range of values of importance ratio; beyond the plateau, for high values of the importance ratio, a steep decrease occurs. This behaviour is likely due to an extreme increase of the neutron weight, causing the sampled neu- trons to undergo excessive splitting, which in turn leads to strong correlations that reduce the FoM [20,28]. Overall, for this benchmark configuration the performance of the improved algorithm is comparable to that of the legacy implementation, except for the first time step, where the improved algorithm performs better in the plateau region. The FoM gain associated with the importance ratio value obtained from the adjoint point-kinetics equations falls within the plateau region, confirming also in this case its optimality.
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Evolution of the FoM gain for every time step simulated in the kinetic scenario presented in Figure 5. The displayed FoM gain values for every time step are computed choosing combing (top) or Russian roulette on precursors (bottom) as the population control strategy. In both cases we investigate a broad range of importance ratio values, spanning from $10^3$ to $10^{12}$. The circle represents the FoM gain obtained using the value of *R* determined from the solution of the adjoint point-kinetics equations.</figcaption></figure>

## 4.3 Analysis of the SPERT III E-core benchmark

In order to further broaden our analysis, we consider a scenario corresponding to a negative reactivity insertion, using the Special Power Excursion Reactor Test III (SPERT III) E-core benchmark [34–36], based on a small pressurized-water research reactor built and operated in the United States in the 1960s to probe the effects of power excursions on fuel and structural materials. The core contains 60 assemblies, including 48 fuel assemblies with 25 (5 × 5) fuel rods, 4 assemblies with 16 (4 × 4) fuel rods, and 4 pairs of borated steel control rods. The SPERT III E-core benchmark was previously used in a

series of investigations carried out with TRIPOLI-4, to which we refer for a thorough description of the reactor components [4,24,37–39]. The corresponding model is illustrated in Figure 7.

Here we will consider a shut-down scenario: starting from a critical (cold zero power) configuration, all the control rods are inserted into the core, leading to an instantaneous reactivity withdrawal of ~$19. Following the control rod insertion, the neutron population first undergoes a prompt drop and then a slower decrease driven by the precursor decay. In view of the position of the control rods in the core, the reactivity insertion is radially symmetric. This kind of kinetic simulation has been
<figure><img src="image_8.png" /><figcaption>Fig. 7. Radial and axial cut of the SPERT III E-core benchmark, as modelled with TRIPOLI-4.</figcaption></figure>

previously investigated using the legacy kinetic algorithm in TRIPOLI-4 and a constant importance ratio [3], and it is thus interesting to revisit it using the improved algorithm, with a customized time grid. Specifically, the time grid is divided as follows:

$$ \Delta t = \begin{cases} 100 \mu s, & \text{for } 0 \le t < 1 \text{ ms,} \\ 1 \text{ ms,} & \text{for } 1 \text{ ms} \le t < 0.011 \text{ s,} \\ 10 \text{ ms,} & \text{for } 0.011 \text{ s} \le t < 0.111 \text{ s,} \\ 0.1 \text{ s,} & \text{for } 0.111 \text{ s} \le t < 1.111 \text{ s,} \\ 1 \text{ s,} & \text{for } 1.111 \text{ s} \le t < 10.111 \text{ s.} \end{cases} \quad (13) $$

This choice is motivated by the fact we want to ensure a finer time grid during the early phases of the simulation, allowing frequent application of population control to maintain a roughly constant and convenient number of particles during the rapid transient. At later times, the grid becomes coarser as the decrease of the flux is regulated by slower time-scale of the delayed neutrons precursors.

To test the simulation capabilities of TRIPOLI-4 on the transient described above, we rely on the importance-sampling scheme with the adaptive calculation of R implemented in the improved algorithm. In this scenario, the improved algorithm automatically adapts the importance ratio according to the solution of the adjoint point-kinetics equations, reflecting the fact that in the first part of the grid (with shorter time steps) the neutrons are more important than the precursors, while in the second part (with longer time steps) the contribution from the precursors is the most relevant.

For comparison, the legacy algorithm offers two approaches for handling the two time-scales of this scenario. We can either use a single constant value for the importance ratio R, or assign a list of importance values for each time step. The use of a constant value is immediately faced with a serious issue: choosing a large value of R to boost the number of neutrons in the first phase of the transient (where time steps are very short and the contribution from precursor decay is not significant) will allocate excessive computational resources to neutrons in the sec-

<figure><img src="image_9.png" /><figcaption>Fig. 8. Importance ratio as a function of the time-step size used in the dynamic simulation. In the new algorithm these values are computed by the adjoint point-kinetics scheme, whereas in the legacy algorithm they are provided by the user.</figcaption></figure>

ond phase, potentially leading to a catastrophic memory burden. Conversely, if a low value of R is chosen in view of optimizing the second phase of the transient, then in the first phase neutrons will be significantly under-sampled, leading to unacceptable statistical uncertainty. In order to overcome the problems posed by using a constant R for the simulation of a long transient in SPERT III, in Ref. [3] two different values of importance ratio were determined through trial and error, one to deal with the prompt drop and one for the rest of the transient. For reference, we will reproduce this configuration in our investigation, setting $R = 10^6$ from $t = 0 \text{ s}$ to $t = 0.001 \text{ s}$, and $R = 10^2$ thereafter.

The values of the importance ratio computed with the improved algorithm for this scenario are shown in Figure 8, as a function of all the time steps sizes listed in equation (13). The guess values used in the legacy algorithm are also displayed. It is apparent that the value of R changes by several orders of magnitude over the transient, significantly more than what was observed for the positive transient described in Section 4.2, where the ratio varies only by a factor of two. This suggests that the adaptive calculation of the importance ratio performed by the improved algorithm might be particularly efficient for this scenario.

The comparison between the improved and the legacy algorithm is shown in Figure 9, using the neutron flux as the observable and the FoM($t_i$) values, computed at each time step $t_i$, as metrics. We consider two cases: combing (Fig. 9i) and Russian roulette (Fig. 9ii). In both cases simulations are performed using 100 cycles and 500 neutrons per cycle.

Our results confirm that the adaptive procedure for the importance ratio calculation provides the most efficient strategy for the FoM throughout the whole simulation. This trend is observed regardless of the choice of the
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** A comparison of the neutron flux (left) and the standard error of the neutron flux (middle) is presented between the legacy algorithm, using two different values of $R$, and the new algorithm, where $R$ is adapted at each change in time step size according to the adjoint-based scheme. The figure on the right displays the FoM ratio between the new and legacy algorithms. This analysis is provided for both Russian roulette (top) and combing (bottom).</figcaption></figure>

population control algorithm. However, while the improved algorithm demonstrates comparable results and uncertainties below 5% in the flux score when using either combing or Russian roulette, the legacy algorithm appears to yield more accurate results when combing is applied, with a maximum uncertainty of approximately 7% compared to 22%. In conclusion, the improved algorithm proves to be a highly effective variance-reduction technique, particularly when combined with the adjoint-based calculation of the importance ratio.

## 4.4 Handling strong reactivity excursions

After evaluating the FoM achieved by the improved algorithm in TRIPOLI-4, our goal is to further demonstrate its practical capabilities. In particular, we would

like to assess its impact on the simulation of strong flux excursions and long transients, which are two of the major challenges for kinetic Monte Carlo simulations. For this purpose, we reconsider the Flattop-Pu benchmark addressed in Section 4.1. Reactivity is again inserted by varying the inner Pu sphere radius. After maintaining a critical state for 1 µs, we inject $2.2 \times 10^5$ in four distinct scenarios, the system is kept in supercritical conditions for (i) 24 µs, (ii) 34 µs, (iii) 44 µs, and (iv) 51 µs, respectively. Adjusting the time duration over which the system remains supercritical, we can generate neutron flux excursions of different amplitudes, spanning several orders of magnitude. According to the point-kinetics equations, the neutron population is expected to increase by a factor of the order of $10^3$ in the first case, $10^4$ in the second case, $10^5$ in the third case, and $10^6$ in the fourth case. At the end of each reactivity insertion, the system is then brought
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Neutron flux obtained on the Flattop-Pu benchmark with the new kinetic algorithm, juxtaposed to the evolution predicted by the point-kinetics equations. The transient is initiated by enlarging the internal Pu sphere, and it consists of a positive reactivity insertion of $2.2$ on a system initially kept in its critical state for $1$ µs. Different durations of the supercritical state are investigated, namely (i) $24$ µs, (ii) $34$ µs, (iii) $44$ µs, (iv) $51$ µs. The evolution of the importance values and the importance ratio used in the new importance-sampling scheme is also presented, demonstrating the highly adaptive and dynamic nature of this new variance reduction strategy.</figcaption></figure>
<table><caption>**Table 2.** Total CPU time, in $kh_{CPU}$, for different kinetic simulations of the transient simulated with the Flattop-Pu benchmark. The system, initially in a critical state, evolves to a supercritical state which is kept for different durations. These simulations were carried out on AMD Milan CPUs (64 physical cores), with a total of 180 independent processes.</caption><thead><tr><th>Duration of supercriticality</th><th>24 µs</th><th>34 µs</th><th>44 µs</th><th>51 µs</th></tr></thead><tbody><tr><th>kh<sub>CPU</sub></th><td>5705</td><td>5917</td><td>7526</td><td>10944</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Neutron flux obtained on the Flattop-Pu benchmark with the legacy kinetic algorithm, juxtaposed to the evolution predicted by the point-kinetics equations. The transient is initiated by enlarging the internal Pu sphere, and it consists of the maintenance of a positive reactivity insertion of $2.2$ for $34~\mu\text{s}$ on a system initially kept in its critical state for $1~\mu\text{s}$. The importance values and importance ratio used in the legacy importance-sampling scheme in this case remain constant, as the legacy algorithm does not support their dynamic adaptation. As a consequence of this less tailored variance reduction technique, the kinetic results from TRIPOLI-4 significantly underestimate the system's evolution.</figcaption></figure>

back to criticality, and the neutron population is tracked over a total simulation time of $100~\mu\text{s}$.

In the first, second and third cases, the time step size used for the application of the population control algorithm is kept equal to $\Delta t = 1/3~\mu\text{s}$, while in the last case a finer grid with $\Delta t = 0.25~\mu\text{s}$ is employed, given the extremely steep increase of the flux during the supercritical regime. All cases were simulated with 1000 cycles and 10000 neutrons per cycle. Our simulation results are shown in Figure 10, where the flux obtained with the TRIPOLI-4 kinetic simulation is compared to the solution of the point-kinetics equations, for reference. The upper and lower bounds of the error band on the point-kinetics curve correspond to solutions obtained with $\rho_k = 1/k_{\text{critical}} - 1/k_{\text{perturbed}}$ varied within its uncertainty, which in turn derives from the statistical uncertainty of the effective multiplication factor estimated through the power iteration algorithm in TRIPOLI-4.

In the same figure, we also present the evolution of the importance ratio and the neutron and precursor importance values, normalized to their initial values. In all the

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Radial cut of the TRIPOLI-4 model of the CROCUS core.</figcaption></figure>

tested cases, the importance ratio remains relatively stable, as there are no significant changes in the time step size or reactivity. In contrast, the evolution of the neutron and precursor importance values spans several orders of magnitude, following the duration of the supercritical state. The significant variation in the importance values, despite the importance ratio remaining nearly constant, arises from the fact that these values depend on the neutron and weight proportions, which in the improved algorithm are continuously updated based on the direct point-kinetics equations (10a) and (10b). This is in contrast to the legacy algorithm, where these proportions are exclusively based on the initial weights (Eqs. 7a and 7b) and are not updated during the simulation. The computational cost of these simulations, expressed in $kh_{CPU}^3$, is presented in Table 2.

The newly developed variance-reduction technique is ideally suited for this challenging scenario, as witnessed by the evolution of the importance values to closely follow the transient. To illustrate the impact of the features added in the new algorithm, Figure 11 presents the results obtained using the legacy kinetic algorithm for one of the cases just investigated, specifically the scenario where the system remains in a supercritical state for $34~\mu\text{s}$. In this configuration, the importance ratio and importance values are not dynamically adjusted during the simulation but remain constant, resulting in a less effective variance reduction and a significant underestimation of the neutron flux.

## 4.5 Handling long transients

To further substantiate our analysis and probe the behaviour of the improved algorithm with respect to scenarios involving long transients, we consider finally the CROCUS benchmark [25], based on a thermal two-region uranium-fueled and light-water-moderated zero-power reactor located at the École Polytechnique Fédérale de Lausanne (EPFL) in Switzerland. The CROCUS core is loaded with two types of fuel rods, using respectively uranium oxide ($UO_2$) at $1.806\%_{\text{wt}}^{235}\text{U}$ enrichment, and

<footnote n="3">³ The $kh_{CPU}$ is the unit used to quantify the total CPU time, defined as the product of the elapsed wallclock time and the total number of processes employed during the simulation.</footnote>
<table><caption>**Table 3.** Details of the three supercritical benchmark configurations analyzed in this work.</caption><thead><tr><th>CROCUS Config.</th><th>Water level [cm]</th><th>Experimental inverse period α [s<sup>-1</sup>]</th></tr></thead><tbody><tr><td>H<sub>1</sub></td><td>96.51</td><td>0.0</td></tr><tr><td>H<sub>2</sub></td><td>98.51</td><td>1.3422 × 10<sup>-2</sup></td></tr><tr><td>H<sub>3</sub></td><td>99.00</td><td>1.8180 × 10<sup>-2</sup></td></tr><tr><td>H<sub>4</sub></td><td>99.51</td><td>2.3392 × 10<sup>-2</sup></td></tr></tbody></table>

<table><caption>**Table 4.** Total CPU time, in $kh_{CPU}$, for different kinetic simulations of the transient tailored on the CROCUS benchmark, where the system, initially in a critical state, evolves to different supercritical states ($H_2$, $H_3$, $H_4$). These simulations were run on AMD Milan CPUs (64 physical cores), with a total of 62 independent processes.</caption><thead><tr><th></th><th>H<sub>2</sub></th><th>H<sub>3</sub></th><th>H<sub>4</sub></th></tr></thead><tbody><tr><td>JEFF-3.3</td><td>2.90</td><td>3.60</td><td>4.10</td></tr><tr><td>ENDF/B-VIII.0</td><td>3.80</td><td>6.82</td><td>9.54</td></tr></tbody></table>

uranium metal at 0.947%$_{wt}$ $^{235}$U enrichment. The experimental configurations that we have used for our Monte Carlo simulations are documented in the international OECD/NEA benchmark collection [25], and have been previously modelled with the TRIPOLI-4 code [40,41] (see Fig. 12).

In this work we consider three transient scenarios, where the system departs from a critical configuration ($H_1$) and attains three different supercritical configurations, called $H_2$, $H_3$, and $H_4$, respectively, by varying the water level inside the core tank. Therefore, the reactivity insertion for these scenarios are again radially symmetric. The water levels, along with the corresponding experimentally measured inverse reactor period $\alpha$, are shown in Table 3. In our simulations, the system is initially critical for 10 s, then suddenly becomes supercritical by changing the water level and stays supercritical for 1 minute; finally, criticality is restored and kept for extra 10 s. We used 3000 neutrons per cycle for a total of 150 cycles. Population control is performed on a time grid with $\Delta t = 0.5$ s.

## 4.6 Comparison between the results obtained with the kinetic algorithm and point kinetics

The neutron fluxes estimated by TRIPOLI-4 are shown in Figure 14 and compared to the solution of the point-kinetics equations, for reference. As in Section 4.4, the uncertainty on the point kinetics curve is derived from the uncertainty of the effective multiplication factors used to compute the inserted static reactivity. Exclusively for this benchmark study, the JEFF-3.3 nuclear data library [42] has been tested alongside ENDF/B-VIII.0 to evaluate the impact of nuclear data on the results. Overall, a good

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Evolution of the CROCUS system from a critical state ($H_1$) to a super critical configuration ($H_2$), which is kept for one minute before going back to criticality. The results of the kinetic simulation with TRIPOLI-4 are compared with the behavior predicted by the point-kinetics equation, using the static reactivity $\rho_k = 1/k_{critical} - 1/k_{perturbed}$ and the k-adjoint weighted kinetics parameters (green) but also the dynamic reactivity $\rho_\alpha$ and the $\alpha$-adjoint weighted kinetics parameters (red).</figcaption></figure>

agreement is observed between the results of the kinetics simulations in TRIPOLI-4 and point-kinetics for the configurations $H_3$ and $H_4$, for both libraries. This is not surprising, given that the CROCUS reactor is a highly compact system where spatial effects are minimal, particularly when reactivity is introduced homogeneously via water level control in the tank. As a general trend, the flux estimated using the ENDF/B-VIII.0 nuclear data library is higher than that estimated using JEFF-3.3, consistent with the fact that the injected reactivity is estimated to be larger when using the former library. This trend is also reflected in the time required to complete the simulation, measured in $kh_{CPU}$, for the three $H_2$, $H_3$, $H_4$ supercritical states, which are consistently higher for simulations using ENDF/B-VIII.0, as shown in Table 4.

The $H_2$ configuration shows the largest discrepancy between the point-kinetics and TRIPOLI-4 kinetics results, as illustrated in Figure 14i. Indeed, in this case the point-kinetics equations seem to underestimate the flux in the system compared to kinetic simulations. Stimulated by this puzzling behaviour, we
<figure><img src="image_15.png" /><figcaption>Fig. 14. Transient simulation, going from the critical configuration called $H_1$ to different supercritical states ($H_2$, $H_3$, $H_4$) by adjusting the water level in the tank. The results of the kinetic simulation in TRIPOLI-4 are compared with the solution of the point-kinetics equation, using the adjoint-weighted kinetics parameters. Results have been obtained both with the ENDF/B-VIII.0 (left) and the JEFF-3.3 (right) nuclear data library.</figcaption></figure>

resorted to a complementary analysis: instead of relying on the static reactivity $\rho_k$ and adjoint k-weighted kinetics parameters, we solved the point-kinetics equations using the dynamic reactivity $\rho_\alpha$, defined by the

equation

$$ \rho_{\alpha} = \alpha \left( \Lambda_{\alpha} + \sum_{j=1}^{6} \frac{\beta_{\alpha,j}}{\alpha + \lambda_{j}} \right), \quad (14) $$
where $\alpha$ is the dominant time eigenvalue, and the adjoint $\alpha$-weighted kinetics parameters $\Lambda_{\alpha}$ and $\beta_{\alpha,j}$, which can be computed using TRIPOLI-4 [43–46]. The corresponding point-kinetics results are illustrated in Figure 13 and align more closely with the behavior predicted by the kinetic simulation for both nuclear data libraries. For the $H_3$ and $H_4$ configurations, the good agreement already observed with the standard point-kinetics approach is also maintained with the dynamic reactivity $\rho_{\alpha}$ and the $\alpha$-weighted kinetics parameters (not shown here). These considerations clearly demand further investigations, since a systematic comparison between k-point kinetics and $\alpha$-point kinetics could provide valuable insights, especially in view of the possibility of using kinetic simulations as a reference, together with experimental results (when available).

## 5 Conclusions

In this manuscript we have presented the development of an improved importance-sampling scheme for time-dependent simulations in TRIPOLI-4, in order to effectively tame the imbalance between the neutron and precursor populations. In the first part of this work we have provided a general overview of this scheme, highlighting the differences with respect to the legacy implementation. We have then discussed a thorough analysis of the algorithm over a set of benchmark configurations and various reactivity insertion and withdrawal scenarios, using the FoM as a metric for the performance evaluation. In the last portion of the paper we have showcased the enhanced kinetic capabilities enabled by the new importance-sampling scheme, demonstrating the simulation of longer and more steep transients previously unattainable with the legacy version of the algorithm. Our findings contribute thus to establishing Monte Carlo methods as a reference simulation tool for the analysis of non-stationary configurations in reactor physics applications. A possible shortcoming of our findings is that we only considered scenarios where the cores are tightly coupled and the reactivity perturbation is radially symmetric, which improves the agreement between point kinetics and full three-dimensional kinetics; this might not be the case for more general scenarios.

Future work will concern the validation of the improved time-dependent capabilities of TRIPOLI-4 against experimental data stemming from experimental campaigns that are currently under way in the CROCUS reactor, encompassing spatially symmetric and asymmetric reactivity insertions and withdrawals.

## Acknowledgments

TRIPOLI-4® is a registered trademark of CEA.

## Funding

This work was partially funded by EDF.

## Conflicts of interest

The authors declare no conflict of interest.

## Data availability statement

The datasets generated and/or analysed during the current study are not publicly available due to them being generated using export-controlled codes, but are available from the corresponding author on reasonable request.

## Author contribution statement

Cecilia Montecchio implemented the new version of the importance-sampling scheme, performed the simulations and wrote the manuscript. Davide Mancusi and Andrea Zoia supervised the methodologies, verified the algorithm implementation, and wrote the manuscript. Vincent Lamirand and Wilfried Monange reviewed the article and co-supervised the work.

## References

1. V. Sanchez-Espinoza et al., The McSAFE project – High-performance Monte Carlo based methods for safety demonstration: From proof of concept to industry applications, EPJ Web Conf. **247**, 943 (2020), https://doi.org/10.1051/epjconf/202124706004

2. D. Ferraro et al., Serpent and Tripoli-4 transient calculations comparisons for several reactivity insertion scenarios in a 3D PWR minicore benchmark, in Proc. M&C2019 (Portland, USA, 2019)

3. M. Faucher, D. Mancusi, A. Zoia, New kinetic simulation capabilities for Tripoli-4: Methods and applications, Ann. Nucl. Energy **120**, 74 (2018)

4. D. Mancusi, M. Faucher, A. Zoia, Monte Carlo simulations of the SPERT III E-core transient experiments, Eur. Phys. J. Plus **137**, 127 (2022)

5. A. Levinsky et al., Modeling of the SPERT transients using Serpent 2 with time-dependent capabilities, Ann. Nucl. Energy **125**, 80 (2019), https://doi.org/10.1016/j.anucene.2018.09.038

6. D. Ferraro et al., Serpent/SUBCHANFLOW pin-by-pin coupled transient calculations for the SPERT-IIIE hot full power tests, Ann. Nucl. Energy **142**, 107387 (2020), https://doi.org/10.1016/j.anucene.2020.107387

7. J. Kang et al., Direct whole core modeling and simulation of the SPERT III E-Core experiments by nTRACER, Prog. Nucl. Energy **139**, 103824 (2021), https://doi.org/10.1016/j.pnucene.2021.103824

8. X. Guo et al., Kinetic methods in Monte Carlo code RMC and its implementation to C5G7-TD benchmark, Ann. Nucl. Energy **151**, 107864 (2021), https://doi.org/10.1016/j.anucene.2020.107864

9. B. Molnar, G. Tolnai, D. Legrady, A GPU-based direct Monte Carlo simulation of time dependence in nuclear reactors, Ann. Nucl. Energy **132**, 46 (2019), https://doi.org/10.1016/j.anucene.2019.03.024

10. M. Faucher et al., Variance-reduction methods for Monte Carlo kinetic simulations, in *International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering, M&C 2019* (2019), pp. 130–139

11. D. Mancusi, A. Zoia, Zero-variance schemes for kinetic Monte Carlo simulations, Eur. Phys. J. Plus **135**, 401 (2020), https://doi.org/10.1140/epjp/s13360-020-00387-8

12. D. Legrady, J.E. Hoogenboom, Scouting the feasibility of Monte Carlo reactor dynamics simulations, in *PHYSOR'08: International Conference on the Physics of Reactors* (Interlaken, Switzerland, 2008)
13. B.L. Sjenitzer, J.E. Hoogenboom, Dynamic Monte Carlo method for nuclear reactor kinetics calculations, Nucl. Sci. Eng. **175**, 94 (2013)

14. B.L. Sjenitzer, Ph.D. thesis, Delft University of Technology, Netherlands, 2013

15. D. Legrady et al., Population-based variance reduction for dynamic Monte Carlo, Ann. Nucl. Energy **149**, 107752 (2020), https://doi.org/10.1016/j.anucene.2020.107752

16. I. Variansyah, R.G. McClarren, Analysis of population control techniques for time-dependent and eigenvalue Monte Carlo neutron transport calculations, Nucl. Sci. Eng. **196**, 1280 (2022), https://doi.org/10.1080/00295639.2022.2091906

17. K. Fröhlicher et al., Improving variance estimation for time-dependent detectors in Monte Carlo dynamic calculations using adaptive sampling of neutrons, in *M&C 2023 – The International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering, American Nuclear Society* (Niagara Falls, Ontario, Canada, 2023)

18. F.X. Hugot et al., Overview of the TRIPOLI-4 Monte Carlo code, version 12, EPJ Nucl. Sci. Technol. **10**, 17 (2024), https://doi.org/10.1051/epjn/2024018

19. E. Brun et al., Tripoli-4, CEA, EDF and AREVA reference Monte Carlo code, Ann. Nucl. Energy **82**, 151 (2015)

20. M. Faucher, Ph.D. thesis, Université Paris-Saclay, France, 2019

21. T. Booth, A weight (charge) conserving importance-weighted comb for Monte Carlo, in *Proc. RPSD Topical Meeting* (Falmouth, USA, 1996)

22. C. Montecchio et al., Towards a highly efficient and unbiased population-control algorithm for kinetic Monte Carlo simulations, EPJ Web Conf. **302**, 09006 (2024), https://doi.org/10.1051/epjconf/202430209006

23. OECD/NEA, International Handbook of Evaluated Credibility Safety Benchmark Experiments, 1995

24. A. Zoia, E. Brun, Reactor physics analysis of the SPERT III E-core with Tripoli-4®, Ann. Nucl. Energy **90**, 71 (2016), https://doi.org/10.1016/j.anucene.2015.11.032

25. OECD/NEA, Benchmark on the Kinetics Parameters of the CROCUS Reactor, Tech. Rep. NEA 4440, 2007

26. C. de Mulatier et al., The critical catastrophe revisited, J. Stat. Mech. Theory Exp. **2015**, 08021 (2015), https://doi.org/10.1088/1742-5468/2015/08/p08021

27. B. Sjenitzer, J. Hoogenboom, A Monte Carlo Method for calculation of the dynamic behaviour of nuclear reactors, Prog. Nucl. Sci. Technol. **2**, 716 (2011), https://doi.org/10.15669/pnst.2.716

28. D. Mancusi, HDR thesis, Université Paris-Saclay, France, 2023

29. E.M. Gelbard, R.E. Prael, Monte Carlo Work at argonne national laboratory, in *Proc. NEACRP Meeting of a Monte Carlo Study Group, ANL-75-2* (1974)

30. R. Brissenden, A. Garlick, Biases in the estimation of Keff and its error by Monte Carlo methods, Ann. Nucl. Energy **13**, 63 (1986)

31. G. Truchet et al., Computing adjoint-weighted kinetics parameters in Tripoli-4 by the iterated fission probability method, Ann. Nucl. Energy **85**, 17 (2015)

32. G.I. Bell, S. Glasstone, *Nuclear Reactor Theory* (Van Nostrand Reinhold Inc., US, 1970)

33. D.A. Brown, M.B. Chadwick, R. Capote et al., ENDF/B-VIII.0: The 8th major release of the nuclear reaction data library with CIELO-project cross sections, new standards and thermal scattering data, Nucl. Data Sheets **148**, 1 (2018)

34. R.E. Heffner, T.R. Wilson, SPERT-III Reactor Facility, Tech. Rep. IDO-16721, Idaho Operations Office, U.S. Atomic Energy Commission, 1961

35. J.E. Houghtaling, J.A. Norberg, J.C. Haire, Addendum to the SPERT-III Hazards Summary Report: Low-Enrichment Oxide Core, Tech. Rep. IDO-17003, Idaho Operations Office, U.S. Atomic Energy Commission, 1965

36. R.K. McCardell, D.I. Herborn, J.E. Houghtaling, Reactivity Accident Test Results and Analyses for the SPERT-III E-CORE – A Small Oxide-Fueled, Pressurized-Water Reactor, Tech. Rep. IDO-17281, Idaho Operations Office, U.S. Atomic Energy Commission, 1969

37. A. Dervaux et al., Benchmarking of the SPERTIII E-core experiment with the Monte Carlo codes TRIPOLI-4®, TRIPOLI-5® and OpenMC, EPJ Web Conf. **302**, 13011 (2024), https://doi.org/10.1051/epjconf/202430213011

38. A. Zoia et al., Analysis of dynamic reactivity by Monte Carlo methods: The impact of nuclear data, Ann. Nucl. Energy **110**, 11 (2017), https://doi.org/10.1016/j.anucene.2017.06.012

39. R. Hourcade, A. Jinaphanh, A. Zoia, Reactor Period Calculations by Monte Carlo Methods: Analysis of the SPERT III E-Core Reactivity Excursions, in *Proceedings of the International Conference on Physics of Reactors (PHYSOR 2022)* (May 15-20, Pittsburgh, PA, 2022)

40. A. Zoia et al., Monte Carlo analysis of the CROCUS benchmark on kinetics parameters calculation, Ann. Nucl. Energy **96**, 377 (2016), https://doi.org/10.1016/j.anucene.2016.06.024

41. Y. Nauchi, A. Jinaphanh, A. Zoia, Analysis of time-eigenvalue and eigenfunctions in the CROCUS benchmark, EPJ Web Conf. **247**, 04011 (2021), https://doi.org/10.1051/epjconf/202124704011

42. A.J.M. Plompen et al., The joint evaluated fission and fusion nuclear data library, JEFF-3.3, Eur. Phys. J. A **56**, 181 (2020), https://doi.org/10.1140/epja/s10050-020-00141-9

43. A. Zoia, E. Brun, F. Malvagi, Alpha eigenvalue calculations with TRIPOLI-4, Ann. Nucl. Energy **63**, 276 (2014)

44. A. Zoia et al., Monte Carlo methods for reactor period calculations, Ann. Nucl. Energy **75**, 627 (2015), https://doi.org/10.1016/j.anucene.2014.09.014

45. V. Vitali et al., Comparison of direct and adjoint k and α-eigenfunctions, J. Nucl. Eng. **2**, 132 (2021)

46. N. Terranova, A. Zoia, Generalized iterated fission probability for Monte Carlo eigenvalue calculations, Ann. Nucl. Energy **108**, 57 (2017), https://doi.org/10.1016/j.anucene.2017.04.014
## Appendix A: Exact solutions of the forward and adjoint point-kinetics equations

We report the exact solution of the forward point kinetics equations, considering only one family of precursors:

$$n(t) = n_0 \left( \frac{\exp(\alpha_1 t) - \exp(\alpha_2 t)}{\Lambda_{\text{eff}}(\alpha_1 - \alpha_2)} + \frac{\alpha_2 \exp(\alpha_1 t) - \alpha_1 \exp(\alpha_2 t)}{\alpha_2 - \alpha_1} \right) \quad (\text{A.1a})$$

$$c(t) = \frac{\beta_{\text{eff}}}{\lambda \Lambda_{\text{eff}}} \left( \frac{\alpha_2 \exp(\alpha_1 t) - \alpha_1 \exp(\alpha_2 t)}{\Lambda_{\text{eff}}(\alpha_2 - \alpha_1)} \right). \quad (\text{A.1b})$$

We also provide the exact solution of the adjoint point kinetics equations evaluated at $t_i$:

$$n^{\dagger}(t_i) = v \left( \frac{\bar{\lambda}\alpha_2 (e^{\alpha_1(t_{i+1}-t_i)} - 1)}{\alpha_1 - \alpha_2} + \frac{\bar{\lambda}\alpha_1 (1 - e^{\alpha_2(t_{i+1}-t_i)})}{\alpha_1 - \alpha_2} + \frac{\alpha_1\alpha_2 (e^{\alpha_1(t_{i+1}-t_i)} - e^{\alpha_2(t_{i+1}-t_i)})}{\alpha_1\alpha_2(\alpha_1 - \alpha_2)} \right) H(t_{i+1} - t_i) \quad (\text{A.2a})$$

$$c^{\dagger}(t_i) = v \bar{\lambda} \left( \frac{\alpha_1 (1 - e^{\alpha_2(t_{i+1}-t_i)})}{\alpha_1 \alpha_2 (\alpha_1 - \alpha_2)} + \frac{\alpha_2 (e^{\alpha_1(t_{i+1}-t_i)} - 1)}{\alpha_1 \alpha_2 (\alpha_1 - \alpha_2)} \right), \quad (\text{A.2b})$$

where $v$ is the neutron velocity, which cancels out when taking the ratio $n^{\dagger}(t_i)/c^{\dagger}(t_i)$, and $H$ is the Heaviside step function. Here $\alpha_1$ and $\alpha_2$ are the roots of the characteristic polynomial associated with the Nordheim equation, namely:

$$\alpha_{1,2} = \frac{(\rho_k - \beta_{\text{eff}} - \Lambda_{\text{eff}}\bar{\lambda}) + \sqrt{(\Lambda_{\text{eff}}\bar{\lambda} - \rho_k \pm \beta_{\text{eff}})^2 + 4\rho_k \bar{\lambda}\Lambda_{\text{eff}}}}{2\Lambda_{\text{eff}}} \quad (\text{A.3})$$

Among all the parameters appearing in equation (A.3), only the reactivity $\rho_k$ changes over the time grid, mirroring the changes in the simulated configuration. All the other parameters are assumed to stay constant during the transient.