Better physics, less noise: inline multi-physics reduces correlation in Monte Carlo reactor simulations

1 Introduction

Nuclear reactor analysis is a multi-physics problem. For example, the production of heat by neutrons inducing fission affects local temperatures which, in turn, affects neutron moderation and thus their propensity to fission nuclei. The accurate resolution of this multi-physics problem is necessary in the design of high-performance nuclear systems while minimising conservatism.

Within the past two decades, coupling of different reactor-relevant physics domains in high-fidelity has become commonplace (Leppänen et al., 2022; Novak et al., 2022; Romano et al., 2021). In this paper we focus on neutronics and thermal-hydraulics in particular. The simplest and most common approach to this problem involves a Picard iteration: given an initial guess for local fission rates and/or fuel and coolant temperatures, neutronics is solved; the power distribution is passed to thermal-hydraulics which is then solved; the coolant density and fuel temperatures are passed to neutronics; and the cycle repeats. This is done until some degree of convergence is reached – many different metrics have been used. The two significant drawbacks of this Picard iteration are its relatively long runtime (depending on the fidelity of the neutronics and thermal-hydraulics solvers) and its conditional stability. The latter is typically handled by an under-relaxation, which may increase runtime if not applied wisely. These approaches and their challenges are very well surveyed by Gill et al. (2017).

One promising variant of Picard iteration is where neutronics is not fully solved before the power distribution is passed to thermal-hydraulics. Both theoretically (Kochunas et al., 2017) and in practice (Gill et al., 2017; Gill, 2021; Kreher et al., 2019; Kreher, 2023) this has been shown to provide unconditional stability to multi-physics calculations and is also much faster than a full Picard iteration. We will use the term proposed by Gill (2021) to denote this approach: inline coupling.

Inline coupling has been investigated alongside Monte Carlo neutronics on several occasions (Gill et al., 2017; Kreher et al., 2019; Gill, 2021; Kreher, 2023). The current work focuses on the sensitivity of the algorithm to choices an analyst might make: the number of neutrons per generation, the frequency of applying multi-physics updates, and whether to continue coupling during active cycles. Some of these questions were addressed in previous works, albeit not in concert, and some of the algorithmic choices or conclusions in previous works even appear to be partially conflicting.

The choice to investigate neutron population size is wholly novel. This is likely to be an important parameter because it strongly impacts the degree of correlation in Monte Carlo statistical estimates (in k-eigenvalue calculations). Although it is perhaps impossible for this choice to meaningfully introduce an instability in the coupling, it is conceivable that it might notably delay or limit convergence. This is partly done in analogy with behaviour observed during coupled Monte Carlo burn-up calculations (Cosgrove et al., 2020a) – the main difference being that burn-up calculations are only conditionally stable (Densmore et al., 2013; Cosgrove et al., 2020b; Cosgrove and Shwageraus, 2021), while inline coupling should be unconditionally stable (Kochunas et al., 2017).

Therefore, the essence of this work is a series of numerical experiments performed to understand the behaviour of coupling algorithms on reactor problems of engineering interest. This should inform the choice of multi-physics strategy when coupling Monte Carlo neutronics to other physics. These experiments are also compared against the traditional strategy of Picard iterations to contextualise the costs and benefits a reactor analyst might expect.

All investigations are performed using the SCONE Monte Carlo code (Raffuzzi et al., 2026) which was recently updated to support inline multi-physics coupling. Details of this implementation will also be provided in the paper.

2 Theory

2.1 Picard iteration

Coupling between thermal-hydraulics and neutronics is generally numerically unstable. This coupling most often takes the form of a Picard iteration. If x ( n ) is the neutronics and thermal-hydraulic solution vector on iteration n, the Picard iteration can be written as:where H is a non-linear operator consisting of both a thermal-hydraulic and neutronic update on the initial solution vector to obtain that for the next iteration. Various criteria have been used to determine when the iterations should terminate, e.g., various norms applied to power, coolant temperature, and fuel temperature vectors (Gill et al., 2017). In this work, where Picard iterations are used, they are performed for a fixed number of iterations rather than using more sophisticated stopping criteria. A schematic of Picard iteration is shown in Figure 1. If performed naïvely for a certain class of problems commonly observed in nuclear reactor analysis, this iteration does not converge. Notably, for a 3D fuel pin problem with realistic parameters, Gill et al. (2017) showed that the power distribution will tend to oscillate between iterations due to Doppler feedback. Fortunately, Picard iteration can be stabilised by applying an under-relaxation:where ω ∈ ( 0,1 ) is the under-relaxation factor. If chosen too small, many Picard iterations may be required to converge the coupled problem; if chosen too large, instabilities may persist. A value of 0.7 was found to be adequate for the pin problem of Gill et al. (2017).

x n + 1 = H x n ,

x n + 1 = 1 − ω x n + ω H x n , (1)

As discussed in the introduction, a further problem of Picard iteration is that it implies fully solving a given physics domain prior to passing its solution to the next solver. It is possible that some of these solvers might be exceedingly computationally expensive – such as Monte Carlo neutronics. This could be seen as a waste of computation given that (at least during the first few iterations) the initial conditions provided by other physics domains are generally not converged. This motivates searching for more economical coupling strategies.

2.2 Inline coupling

Inline coupling, in the context of neutronics calculations, is where the neutronics solver is not fully converged before power density information is passed to another solver. This has been empirically observed to stabilise coupled physics solves (Gill et al. 2017; Kreher et al. 2019; Gill, 2021) and has been theoretically shown to do so for deterministic neutronics solvers (Kochunas et al., 2017). The practical implementation of inline coupling consists of performing one or several neutronics iterations (also known as cycles or generations in the context of Monte Carlo) to partially converge the fission source before pausing the calculation, passing power density information to thermal-hydraulics, updating the temperature and density fields output from thermal-hydraulics, and then resuming neutronics from where it paused. A schematic of inline coupling is shown in Figure 2. As well as the mentioned theoretical benefit of numerical stability, this provides computational efficiency by requiring only a small fraction of a neutronics solve between physics updates. It also avoids the initialisation cost of neutronics which can sometimes be non-negligible, e.g., loading in nuclear data during a burn-up calculation.

Most attention has been paid to inline coupling in the context of using a Monte Carlo neutronics solver. In part, this may be because deterministic solvers often require acceleration techniques, resulting in a more sophisticated multi-physics strategy where the low-order neutronics solver communicates with other physics (Shen et al., 2019; Shen and Kochunas, 2021; Shen et al., 2021) (a similar strategy has been attempted with Monte Carlo (Harper et al., 2022)). This might also be because there are additional ambiguities in the application of inline coupling when there is noise in the neutronics solution.

There is some choice in how inline coupling may be performed. The two related degrees of freedom are “How many neutronics iterations should be performed between physics updates?” and “How should one accumulate power density results to pass to thermal-hydraulics?” One could update once per power iteration or less often, and one could average results between updates or using a moving window (Kreher et al., 2019; Kreher, 2023). In the context of deterministic neutronics, varying either option appears to have, at best, a weak effect on the stability of the inline coupling. Cosgrove and Shwageraus (2022) analysed the stability of the typical xenon equilibrium iteration (Griesheimer, 2010; Isotalo et al., 2013) (which is performed analogously to inline thermal-hydraulics) considering results accumulation both over a batch of J cycles between updates and updating each cycle with results accumulated from a moving window of width J. In both cases the convergence of the coupled system was found (for more problematic low frequencies) to be unaffected by the value of J.

For Monte Carlo, however, one might generally expect some benefit from not updating frequently. This is because the estimate of the power density distribution obtained on a given cycle will generally be noisy. By accumulating these estimates over several cycles, the deleterious effects of noise may be mitigated. Gill (2021) examined the impact of varying the number of Monte Carlo cycles over which results are accumulated and after which the physics is updated, while Kreher et al. (2019) considered accumulating results over a moving window while performing a physics update every cycle. Gill (2021) concluded that the performance of inline coupling was essentially independent of the update frequency. Kreher et al. (2019) found that the moving window approach resulted in a smoother convergence to a reference solution than using a single-cycle estimate of power densities. Subsequently, Kreher (2023) found that a ‘ramp-up’ strategy (where the number of cycles per update was increased as the coupling proceeded) was stable and efficient.

The more expensive the other solvers (or the more expensive the communication), the less economical frequent inline coupling will become. An inline coupling approach is only practical when the neutronics solution is significantly more expensive than the solution of other physics domains. This would not be the case if using an especially expensive thermal-hydraulics solver; for example, one based on computational fluid dynamics (CFD). If the other physics is, in fact, comparably expensive with Monte Carlo then an inline approach will reduce the efficiency of calculations as the other solver will be called many times per full neutronics solve. In such a situation either a standard, underrelaxed Picard iteration may be more appropriate or one might use a low-order surrogate for the additional physics, e.g., a subchannel solver which preserves most of the neutronics-relevant physics of high-order CFD.

One aspect of inline coupling not investigated in previous works is the effect of the neutron population. A Monte Carlo simulation should ideally use a small number of neutrons each cycle in order to converge the fission source in the least computational time, assuming the cycle-wise convergence rate of the fission source is independent of the population size. In practice, however, using a small neutron population can result in several deleterious effects with respect to tallies: correlation (Ueki, 2002; Ueki et al., 2003; Herman, 2014) and neutron clustering (Dumonteil et al., 2014; Sutton and Mittal, 2017) tend to become more severe as the population becomes smaller, requiring disproportionately more cycles to obtain low uncertainty results. Hence, one might expect the power density tallies and the coupling convergence behaviour to be sensitive to this choice. The effect of population will be investigated by examining the convergence behaviour of a coupled 3D fuel pin problem with simulations performed using SCONE.

There is another aspect of inline coupling which has been examined but where the recommendation is unclear: whether coupling should continue into the active cycles. Coupling entails the costs of communicating results between solvers and running the other solver. Depending on the computational setup, this may be non-negligible. As such one would prefer to terminate coupling once some acceptable degree of convergence has been achieved. However, for inline methods, there may be significant noise in the solution passed to thermal-hydraulics. If the resulting temperature and density fields were then fixed, a bias may ensue. Gill (2021) produced all results by continuing coupling during the active cycles, while Kreher et al. (2019) suggests that the correlation during tallying is significantly worsened by inline coupling and thus deactivates it. However, Kreher (2023) highlights that an error was identified in the coupling methodology used by Kreher et al. (2019) that might affect the conclusions drawn with respect to correlation.

2.3 Implementation of multi-physics in SCONE

SCONE is a Monte Carlo code which supports neutron physics commonly used in fission reactor analysis (Raffuzzi et al., 2026). Recently, it has been modified to accept temperature and density fields super-imposed upon the geometry. This is similar to the multi-physics interface capability of Serpent (Leppänen et al., 2025). One difference, at present, is that surface tracking with an imposed physics field does not rely upon rejection sampling. Instead, the distance until a change in the field value occurs is explicitly computed and the neutron is moved to the next field region if that distance is shorter than the distance to collision and distance to the nearest surface. This limits the imposed fields to being piecewise constant and requires the geometry of the imposed field to be endowed with a distance-to-surface calculation. This may be relaxed in future for the Monte Carlo solver in SCONE to allow arbitrary fields (Leppänen, 2013), but this constraint is useful for upcoming multi-physics coupling using the random ray method as the neutronics solver (Tramm et al., 2017; Cosgrove and Tramm, 2024).

As well as imposed fields, a ‘coupling admin’ has recently been developed in SCONE. This is an object which manages the communication between SCONE’s physics (transport, tallying, imposition of temperatures and densities) and external coupling scripts. For loose coupling such a feature is unnecessary: communication between SCONE and other physics could be managed entirely by an external script, with SCONE terminating between Picard iterations. However, for inline coupling, the coupling admin allows significant control over SCONE without having to terminate the Monte Carlo calculation between physics updates. The frequency of updates can be controlled and one can specify arbitrary tallies used exclusively for coupling. The admin is also told where updated physics fields can be found and how to name output tally information to interact with external scripts. It can be specified whether coupling occurs only during inactive cycles or continues during the active cycles. The coupling can currently be terminated either by reaching some set number of physics updates or by receiving a signal from the external coupling script. All communication between SCONE and external codes is managed by exchanging communication files.

SCONE already supports on-the-fly temperature changes in different materials. In the fuel Doppler broadening is handled by target motion sampling (Viitanen and Leppänen, 2017) and in moderators S ( α , β ) data is handled by stochastic interpolation (Richards, 2020). In order that the majorant cross section need only be constructed on initialisation, the maximum expected temperature in the system must be an input to the calculation. Likewise, for S ( α , β ), an appropriate range of libraries at different temperatures must be supplied in the input such that the moderator temperature at any point does not fall outside of this range.

3 Numerical investigations and results

The behaviour of inline coupling is investigated using a fuel pin problem, comparable to that described in Gill et al. (2017) and Gill (2021), the specifications of which we repeat in Table 1. The original paper did not specify the fuel enrichment; this has been taken to be 2.2%. The table also includes several thermal-hydraulic properties of water which were assumed constant for the simple calculations. The neutronic boundary conditions are radially reflective and axially vacuum. This problem was chosen for two reasons. First, it is comparable with problems used in previous studies (Gill et al., 2017; Kreher et al., 2019; Gill, 2021; Kreher, 2023) – each of these used a PWR pin problem. Second, the studies performed in this work are computationally intensive: larger problems would require more inactive cycles to achieve convergence and larger neutron populations to achieve tightly converged tallies globally. Given the available computing resources, this is tractable for a fuel pin, but would be challenging for larger systems.

TABLE 1. Parameters of the fuel pin model.

Thermal-hydraulics is solved using a simple low-order script. This amounts to performing an enthalpy balance along the 36 axial nodes over which the pin is discretised. The coolant temperature and density are updated along these nodes given an imposed axial power profile. This is then used to calculate the fuel surface and centreline temperatures with the fuel conductivity given by the Babcock and Wilcox correlation (Todreas and Kazimi, 1990). The details of this calculation are given in Appendix. The fuel temperature in neutronics is approximated as radially uniform. The ‘average’ temperature is computed from the fuel surface and centreline temperatures using the Rowlands correlation (de Kruijf and Janssen, 1996) Table 1:

T = 5 9 T surface + 4 9 T centre .

The Monte Carlo simulation uses an imposed temperature and density field across 36 axial zones of equal height across the active fuel, identical to those used for thermal-hydraulics calculations. Power is tallied across the same regions, assuming local energy deposition. The calculations used ENDF/B-VIII.0 data (Brown et al., 2018).

The investigations begin by producing a reference solution. This is done by standard Picard iteration with an under-relaxation factor of ω = 0.7 (following Equation 1), iterated 10 times. For this problem, 10 iterations should be sufficient given the same problem achieved convergence with the same under-relaxation factor in Gill et al. (2017). Each Monte Carlo simulation uses 2 million neutrons per cycle, 200 inactive iterations, and 200 active iterations. The number of inactive cycles was chosen conservatively to achieve source convergence. The number of active cycles was arbitrarily chosen to produce a large number (20 billion) active neutrons each Picard iteration. This entire coupled calculation was repeated with 10 independent runs. The resulting mean power, temperature, and density profiles with 1 σ error bars are shown in Figure 3.

Given the reference solution, we investigate how the population settings and continuation of coupling into active cycles affect the convergence to this solution using inline coupling. We vary three parameters: the frequency of physics updates measured in Monte Carlo cycles, α, the fraction by which we scale the cycle-wise population relative to the reference population of 2 million, γ, and the binary choice of continuing coupling during active cycles. The values these may take are α ∈ { 1,2,5,10,20,50 } and γ ∈ { 0.01 , 0.02 , 0.05 , 0.1 , 0.2 , 0.5 , 1 }. The total number of active neutrons simulated is kept constant between all simulations by varying the number of active cycles depending on the neutron population; given a fixed number of inactive cycles, the cost of these will vary for different values of γ. Each of these calculations was repeated five times with different random number seeds. In addition to simulations using inline coupling, a set of Picard calculations were performed. These calculations use the same settings as the reference calculation, except that their population per cycle is scaled by γ. As for the inline calculations, each Picard run was also performed five times with different random number seeds.

We first show, in Figures 4, 5, the results from the set of Picard iterations. These demonstrate the extent of ensemble error one could expect from common practice multi-physics simulations. Despite each calculation using the same number of active neutrons, there is a significant increase in error (and the variance of the error) as the number of neutrons per cycle is reduced. Going from 20,000 neutrons per cycle to 1,000,000 neutrons per cycle, the maximum power error plus 1 σ decreases from 4% to just over 1%. In the same range of neutrons per cycle, the maximum fuel temperature error plus 1 σ decreases from about 1% to about 0.4%. This is a result of a higher true variance when fewer neutrons per generation are simulated; this increased variance becomes a bias when propagated through the Picard iteration.

We next examine the errors resulting from inline coupling. Across the range of calculations performed, we compute an integral error metric for ease of comparison. This is defined as:where x is the solution, x ref is the reference solution, μ denotes a mean vector (normalised to either total power or the axial average temperature), σ 2 denotes a relative variance vector, the m subscript denotes an element of the vector, and M is the total number of elements. The means and variances are computed over the ensemble of calculations using the same settings. The power errors over the range of α and γ values for the two inline approaches are shown in Figure 6. The corresponding fuel temperature errors are shown in Figure 7. The bottom row of each plot also shows the same metric evaluated for Picard iterations. We emphasise that γ = 1 for Picard iterations is statistically identical to the reference calculation.

ϵ = 1 M ∑ m = 1 M μ x , m − μ x ref , m 2 + σ x , m 2 + σ x ref , m 2 ,

Both power and temperature results suggest that there are significant benefits in continuing coupling during active cycles. For power and temperature, errors of the order of 1 0 − 4 and 1 0 − 3, respectively, are achieved with coupling during active cycles for all values of α and γ. Conversely, without continuing coupling, the achievable error can deteriorate by an order of magnitude, depending on the settings. Without coupling during active cycles, the error decreases as the number of neutrons per cycle is increased, and there is an optimal update frequency for a given population. The former effect is a result of having a thermal-hydraulics field less contaminated by noise. The latter effect is because lower uncertainty results are obtained with fewer updates, but, in the limit of very few updates, this can result in an unconverged coupled solution. The same error phenomenon does not visibly appear when coupling during active cycles because additional coupling updates occur; although the solution may not be fully converged when active cycles begin, it later converges and subsequently accumulated statistics appear to disguise the bias that would have resulted.

Without continuing coupling during active cycles, a Picard calculation is typically much more accurate than inline coupling for the same number of neutrons per generation. However, if active coupling is performed then inline calculations tend to be significantly more accurate than Picard when using relatively few neutrons per generation; otherwise their accuracy is comparable. For selected settings with both approaches, this is shown in terms of the axial power error in Figure 8. We emphasise that for each Picard calculation, the γ-equivalent inline calculation was approximately 10× faster. This is because only one transport solve is required by inline coupling whereas Picard calculations each performed ten full transport solves and the thermal-hydraulics solve and data exchange consumed negligible computing time.

This increase in accuracy for active inline coupling relative to Picard is unexpected: naïvely, for the same Monte Carlo population settings, one might expect comparable accuracy given a sufficient number of Picard iterations or inline updates to converge the physics. We hypothesise that the discrepancy arises from a reduction in inter-generational correlation under inline coupling.

In Picard iteration, the fission source evolves over many cycles against a fixed thermal-hydraulic field, allowing the well-documented correlation between successive generations in Monte Carlo criticality calculations (Ueki, 2002; Ueki et al., 2003; Herman, 2014; Dumonteil et al., 2014; Sutton and Mittal, 2017) to build up largely unimpeded. With fewer neutrons per generation, this produces a higher true variance result than if one used many neutrons per generation (for the same total number of active neutrons). As discussed earlier, a higher variance result fed to thermal-hydraulics induces a larger bias.

Under inline coupling, by contrast, the thermal-hydraulic field itself is repeatedly modified by the noisy power tallies during the course of the solve. We believe that this continual adjustment disrupts the build-up of inter-generational correlation, lowering the variance of the resulting estimates relative to a Picard calculation with an equivalent neutron population; perhaps this might be viewed as the physics ‘helping’ the power iteration to keep neutrons distributed according to the (coupled) fundamental mode. Similar effects have been observed when applying source convergence acceleration during active cycles (Li, 2017). This effect would be expected to diminish as the neutron population increases, since inter-generational correlation is itself weaker at high population (Herman, 2014) – consistent with the convergence of inline and Picard accuracy observed at high γ in Figures 6, 7.

We attempt to substantiate this hypothesis by presenting the autocorrelation coefficients (ACCs) of the spatially-averaged power production in the fuel pin. This is done both with a fixed thermal-hydraulic field and with inline coupling. Both cases have γ = 0.01, i.e., use 20,000 neutrons per cycle and 20,000 active cycles. The inline calculation has α = 1, i.e., a multi-physics update is performed each cycle. The fixed or uncoupled case corresponds to the standard Monte Carlo power iteration used in a single Picard iteration and uses the converged temperature and density fields from the reference calculation. The lag k ACC of a random variable, X, which has n samples, X n, and mean μ x is given by:a value of ρ k closer in magnitude to zero implies less correlation between realisations of X separated by k samples. We follow the approach described by Herman (2014) where the ACC is computed for each power tally in a given axial region, m, and then averaged, weighted by the tally’s variance on a given run, r, σ r m 2. For R independent runs and M regions, the average variance is then computed as:this gives a problem-averaged ACC:this can be related to the expected root-mean-square (RMS) error of the problem as (Herman, 2014):here n corresponds to the number of active generations simulated. This analysis assumes that the ACCs depend only on the lag, not on the generation, n ¹ . Ten independent calculations are run for the coupled and uncoupled cases to compute the ACCs. The results are shown in Figure 9 where the coupled case is clearly subject to less correlation across all lags.

ρ k = Cov k X Var X = E X n − μ x X n + k − μ x Var X n Var X n + k . (2)

σ ̄ 2 = 1 R M ∑ r = 1 R ∑ m = 1 M σ r m 2 . (3)

ρ ̄ k = 1 R M σ ̄ 2 ∑ r = 1 R ∑ m = 1 M σ r m 2 ρ rmk . (4)

E ϵ RMS = σ ̄ 2 n 1 + 2 ∑ k = 1 n − 1 1 − k n ρ ̄ k (5)

From the simulations, the RMS error (of power tallies normalised by the total power) is also computed as a function of the number of active neutrons. This is done by storing the cycle-wise power tallies of each calculation and post-processing to obtain the RMS error up to a given cycle. These results are shown in Figure 10, alongside the fitted expected RMS error from Equation 5. Although the results are subject to noise, the expected values agree relatively well with theory, suggesting that inline coupling does not suffer from correlation as much as power iteration. This implies that even if the precise thermal-hydraulic conditions are known, for the same Monte Carlo population per generation, one can obtain a more accurate result from an inline coupled simulation than from an uncoupled calculation. The expected error reduction for the same number of active particles is shown in Figure 11. Although there are marked fluctuations, in the asymptotic regime (where error decreases proportionally to 1 / N) the error reduction is approximately 1.7× from theory, and potentially greater in practice. This corresponds to requiring 2.89× fewer active particles with a coupled calculation to achieve the same error as with an uncoupled calculation. Whether this approach is worthwhile in practice will depend on the computational expense of the thermal-hydraulics solve and data exchange.

4 Conclusion

This work investigated the inline multi-physics coupling between Monte Carlo neutronics and thermal-hydraulics applied to a 3D PWR fuel pin. The conclusions are that this approach should be applied both during inactive and active cycles; otherwise there is a significant error penalty from fixing the thermal-hydraulics as a result of relatively noisy Monte Carlo tallies. Coupling during active cycles was found to be robust to the choice of neutron settings and update frequency for the PWR pin problem. More surprisingly, it was observed that inline coupling can give markedly lower errors than Picard iteration for the same neutron population per generation. This was shown to be due to inline coupling reducing inter-generational correlation between neutrons, and thus feeding less noise to the thermal-hydraulics solver.

The extent to which correlation reduction may occur across different problems has not been investigated here, but will be the subject of future work. Conceivably, this will depend on the physical strength of the feedback and on the size of the problem. It might also depend on whether the reactivity feedback is positive or negative. The dependence of statistical convergence on feedback strength is similar to what Kochunas et al. (2017) describe in a deterministic context: the strength of feedback can accelerate the convergence of a solve beyond that of mere power iteration. Developing a stochastic generalisation of such an analysis would be worthwhile to ascertain generally how multi-physics affects statistical convergence.

Future work could also consider larger systems. Thermal-hydraulic coupling strongly affects the axial power profile across a core, but the radial impact of coupling is likely to be weaker. Realistic reactor problems will have a larger radial extent but the same axial extent. Therefore, it is unclear whether the degree of correlation reduction when performing inline coupling on such systems will be comparable to that observed in this work.

The findings of this work seem likely to apply to enforcing xenon equilibrium as well (Griesheimer, 2010; Isotalo et al., 2013). Indeed, this is a commonly applied ‘inline’ coupling approach, used by both Serpent (Leppänen et al., 2025) and MC21 (Griesheimer et al., 2024). While we believe both implementations continue xenon updates during active iterations (and Serpent also does so for critical boron), this work provides additional justification for this choice.

We plan to use the developed coupling numerics in SCONE on a variety of problems in the near future. As opposed to the simple coupling script and thermal-hydraulics used here, this will be done with more mature coupling frameworks and higher-fidelity solvers.

Correction note

A correction has been made to this article. Details can be found at: 10.3389/fnuen.2026.2001814.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Mark D. DeHart, Abilene Christian University, United States

Reviewed by: Wen Yin, Harbin Engineering University, China

Khurrum Saleem Chaudri, Khalifa University, United Arab Emirates

Asmaa AboElnour, Egyptian Atomic Energy Authority, Egypt

Footnotes

Li (2017) discusses the inapplicability of this formula when acceleration is applied. This is because the expression assumes the ACC is independent of the cycle and only depends on the lag. Acceleration can, in principle, interfere with this–inline feedback may have the same effect. As performed by Li (2017), we investigated computing the ACCs over several subsets of the 20,000 active cycles and found little variation in their average. Hence we believe that the ACCs presented are valid for inline feedback.