SUPER-HISTORY METHODS FOR ADJOINT-WEIGHTED TALLIES IN MONTE CARLO TIME EIGENVALUE CALCULATIONS

Time eigenvalues emerge in several key applications related to neutron transport problems, including reactor start-up and reactivity measurements. In this context, experimental validation and uncertainty quantification would demand to assess the variation of the dominant time eigenvalue in response to a variation of nuclear data. Recently, we proposed the use of a Generalized Iterated Fission Probability method (G-IFP) to compute adjoint-weighted tallies, such as kinetic parameters, perturbations and sensitivity coefficients, for Monte Carlo time (or alpha) eigenvalue calculations. With the massive use of parallel Monte Carlo calculations, it would be therefore useful to trade the memory burden of the G-IFP method (which is comparable to that of the standard IFP method for k-eigenvalue problems) for computation time and to rely on history-based schemes for such adjoint-weighted tallies. For this purpose, we investigate the use of the super-history method as applied to estimating adjoint-weighted tallies within the α-k power iteration, based on previous work on k-eigenvalue problems. Verification of the algorithms is performed on some simple preliminary tests where analytic solutions exist. In addition, the performances of the proposed method are assessed by comparing the super-history and the G-IFP methods for the same sets of benchmark problems.


INTRODUCTION
Several key applications in reactor physics, encompassing reactor period measurements in nearly critical configurations, reactor noise analysis techniques such as the Rossi alpha or Feynman alpha in steady-state configurations, or material control and accountability in critical assemblies, require the assessment of the asymptotic evolution of the neutron population [1,2]. This is usually achieved by determining the dominant time (or α) eigenvalue associated to the Boltzmann operator for the problem under analysis [1]. In this context, experimental validation, uncertainty quantification and bias estimation through data assimilation may further require computing the perturbation and sensitivity coefficients of the dominant α eigenvalue with respect to nuclear data. This topic has recently drawn attention both in deterministic [3,4] and Monte Carlo methods [5,6]. For this purpose, we have recently implemented and successfully tested time-eigenvalue perturbation and sensitivity analysis techniques [7] using the Generalized Iterated Fission Probability (G-IFP) method [8]. Similarly to the standard Iterated Fission Probability (IFP) in k-eigenvalue calculations, the major drawback of the G-IFP implementation in time-eigenvalue problems is the strong increase in memory footprint. In order to overcome this limitation, in this work we introduce a novel adjoint-weighted tally capability based on the super-history method. This method is inspired from previous findings for k-eigenvalue problems [10]. The proposed algorithms will be tested for timeeigenvalue sensitivity coefficients and compared to the G-IFP method.

THE SUPER-HISTORY METHOD IN THE α-k POWER ITERATION
The α-k power iteration, as described in e.g. [11], is basically a modified version of the standard power iteration where neutrons are tracked along successive cycles from birth to death by absorption, leakage or α-reactions (either time-absorption or time-production, depending on the sign of α). Individual neutrons may generate fission or alpha-copy neutrons that will be promoted to the next cycle. An improved version of the α-k method that ensures numerical stability for negative α has been introduced in [12,13]. In this simulation scheme, a cycle matches a neutron generation (alpha-copies also contribute to the next generation). The value of α is updated at every cycle in order to make the fictitious k of the underlying power iteration converge to 1, and the neutron population is finally normalized at the beginning of every cycle.
For k-eigenvalue problems (i.e., for the standard power iteration), the super-history method consists in letting the neutron population evolve over several inner generations without normalization [14]: for the super-history method, a single power iteration cycle contains several consecutive inner generations, and the normalization of the neutron population is done at the beginning of every cycle. The number of inner generations per cycle is denoted by L and is a user-defined parameter.
Since the iteration scheme of the α-k is very similar to the power iteration, it is in principle possible to adapt the super-history algorithm to the α-k method (for an illustration, see Fig. 1 (left)): the neutron population would evolve over several inner generations, the value of α being fixed until the following cycle is reached. We assign a generation identifier g ≤ L to each neutron pertaining to the same generation. For a given cycle c, the starting neutrons for the current cycle have a generation identifier g = 0. Then, the direct progeny of a given neutron, i.e., the neutrons promoted to the following generation either by fission or alpha-copy, will be tagged with a generation identifier increased by one. Finally, neutrons will be promoted to the next cycle when their generation identifiers reach g = L.
If the super-history method is naively applied to the α-k iteration scheme, numerical tests have shown that the neutron population during the convergence phase (i.e., the inactive cycles of the the α-k scheme) can grow very large or decay within a cycle due to the fact that the cycle-wise α value may be far from the true value: in this case, the apparent multiplication factor may be very far from one and thus lead to undesired growth or complete die-away of the neutron population. In order to circumvent this problem, we choose to trigger the super-history method once the convergence of the α-k is attained: the super-history becomes then numerically stable, as the apparent multiplication factor k of the system is very close to one and the neutron population is less likely to grow or decay from one cycle to the next.

ADJOINT-WEIGHTED TALLIES FOR SUPER-HISTORY SCHEMES
First-order perturbations, sensitivity coefficients or kinetics parameters for time-eigenvalue problems involve adjoint-weighted quantities, similarly as for their k-eigenvalue counterparts [7,8]. For example, the sensitivity coefficients of the α eigenvalue with respect to a given parameter x p can be expressed as where ∂ x p applied to a generic operator A denotes the scaled derivative B p is the prompt part of the Boltzmann operator, F d the delayed operator, λ j the delayed neutron precursor decay constant of family j, and is the speed of the particle. In addition, ϕ α and ϕ † α represent the direct and adjoint fundamental α-eigenmodes, respectively. Bilinear forms involving both ϕ α and ϕ † α can be successfully computed by using the G-IFP method [8], which has been adapted to time-eigenvalue problems by extending the IFP algorithm proposed for k-eigenvalue problems [15,16]. The key idea is that the fundamental adjoint function ϕ † α can be determined by estimating the neutron importance in the system, which is proportional to the neutron progeny after a large number of generations. As a consequence, ϕ † α can be estimated by recording the descendants after M latent generations for an ancestor injected into the system at coordinates r, Ω, E. A neutron injected in the system will yield on average neutrons at a latent generation M. Here, k 0 denotes the cycle-wise estimation of the apparent multiplication factor of the system (close to one as we are in a α-eigenvalue iteration), which is constant throughout the cycle, and ω i the weight of the neutrons at generation M. Similarly as for IFP, the G-IFP method is known to be possibly hindered by a large memory footprint, which is proportional to the number of neutrons per cycle, to the number of latent G-IFP generations * and to the number of required adjoint-weighted tallies. This memory overhead may become overwhelming when a large number of sensitivity coefficients are required [7].
In order to cope with the memory burden of the IFP method, for k-eigenvalue problems the superhistory method has been successfully applied to compute adjoint-weighted tallies [10] and sourceterm perturbations [17] in the context of k-eigenvalue sensitivities and perturbations: this technique is inherently history-based and thus allows strongly reducing the memory footprint (at the expense of increased computer time). Based on these findings, we have chosen to adapt the super-history method to α-eigenvalue problems. The super-history method simulates a neutron and its progeny over L generations before simulating the next neutron: as such, it easily allows computing the random contribution (associated to the ancestor neutron) together with its importance at the L-th distant generation (associated to the progeny) for the required adjoint-weighted tallies. As opposed to the G-IFP, we are no longer required to keep in memory individual random contributions over several cycles, and the contribution of the current simulated neutron to the adjoint-weighted tallies can be computed without storing further information. A super-history-based scheme for adjointweighted tallies is illustrated in Fig. 1 (right). Random contributions are computed only for the neutron whose generation identifier is zero (the ancestor). Then, once its whole progeny has been simulated, the ancestor random contribution is weighted by the total weight of the progeny (whose generation identifier is L): the number of latent generations for adjoint-weighting is set to L − 1.
The specific tallies needed for perturbations, sensitivity coefficients or kinetics parameters have been thoroughly described in [7,8].

VERIFICATION TESTS
Verification tests for the super-history scheme in α-eigenvalue calculations have been performed on two simple benchmark configurations where analytical solutions are available. The first is a twoenergy group transport model (with two precursor families) in an infinite multiplying medium. The second is a continuous-energy neutron thermalization problem in an infinite medium.

Two-energy group transport in an infinite medium
We consider a homogeneous infinite medium with two energy groups and two precursor families: the full model description and the exact solutions can be found in [8]. The dominant α eigenvalues computed by applying the super-history method with L = 5 inner generations within the α-k power iteration for several reactivity levels are displayed in Tab. 1. A good agreement with respect to analytical solutions is found (the last considered case lies outside the 3-σ confidence interval, but the α-k method is known to have large cycle-to-cycle correlations, which possibly leads to an underestimation of the error bars). These findings show that the super-history method is unbiased for direct time-eigenvalue calculations.
* Similarly as for the standard IFP method [9], we have implemented an overlapping-cycles G-IFP method for the adjoint-weighted tallies in time-eigenvalue problems [7].  As for adjoint-weighted tallies, the Monte Carlo results using the super-history method (with L = 5 inner generations) are shown in Tab. 2 for the sensitivity coefficients S α,x p corresponding to a deep sub-critical configuration including delayed neutrons. The simulation results are in very good agreement with reference solutions: relative discrepancies lie below 0.2%, and the analytic solution is in the 3-σ confidence interval.

Thermalization problem
In order to probe continuous-energy transport, we have further considered a homogeneous infinite system made of an ideal gas of non-fissile nuclei. Scattering is assumed to be at thermal equilibrium, with a so-called amnesia kernel (capable of thermalizing the neutrons in a single collision), with a Maxwellian spectrum f s (E) ∝ √ E exp(−E/θ), θ being the average post-collision energy. The corresponding α-eigenvalue equation reads Choosing a constant scattering cross section and a 1/ capture cross section allows analytically solving this configuration for the dominant ϕ α (E) and ϕ † α (E), and thus obtaining the sensitivity coefficients [7]. The sensitivity profile S α,Σ c (E) of α with respect to the capture cross section is displayed in Fig. 2: for comparison, we have computed S α,Σ c (E) using the super-history (with L = 5 inner generations) and the G-IFP method (with M = 5 latent generations) and contrasted both to the analytical solution. The sensitivity profile to the capture cross section has been computed by taking 1000 energy groups g for the intervals [E g−1 , E g ], i.e., 1000 tallies S α,Σ c (g) to be estimated within the same calculation. A very good agreement is found between the Monte Carlo calculations and the reference solution. Larger fluctuations are observed for E θ, due to poor statistics in the tails of the post-scattering distribution.

PERFORMANCE ANALYSIS
In order to compare the performances of the super-history and G-IFP methods as applied to computing adjoint-weighted tallies, we have considered the continuous-energy transport problem described above. In particular, we were interested in assessing standard deviations (for several tallies), single-processor computation time and memory overhead. All calculations have been performed on a Intel(R) Xeon(R) CPU E5-2630 v4 @ 2.20GHz using a CentOS 7.6 operating system. Time has been estimated on a non-parallel execution. Memory storage has been estimated through the resident set size of the "ps" Linux command. The performances, expressed in terms of memory footprint and computation time, of the G-IFP and super-history methods as a function of the number of inner cycles or equivalently latent generations (L and M, respectively, according to the notation above) are illustrated in Fig. 3: the memory burden of the G-IFP grows linearly with the number of latent generations M, as expected, since an overlapping-cycle scheme is used for the G-IFP. On the contrary, the memory occupation of the super-history method does not depend on the latent generations L and stays constant. The opposite behavior is observed for computation time: the super-history method displays a linear increase with respect to the number of latent generations, as expected, whereas the G-IFP is roughly constant.

CONCLUSIONS
The super-history method has been successfully combined with the α-k power iteration in order to compute adjoint-weighted tallies, such as the sensitivity coefficients of the time-eigenvalue with respect to nuclear data. However, during the convergence phase (i.e., the inactive cycles of the α-k algorithm), the super-history method might possibly lead to abnormal termination due to a large discrepancy between the current value of α and the target value. Thus, it is advised to activate the super-history method only after the convergence of the α-k is attained. The major advantage of the super-history method as applied to computing adjoint-weighted tallies is that the neutron progeny is simulated within each cycle: contrary to the G-IFP method, the super-history method does not need to keep track of every progenitor for several latent cycles, and the random contribution to adjoint-weighted tallies can be computed separately for each neutron history. These considerations allow considerably reducing the memory occupation, at the expense of an increased computer time. Future work will consist in testing this methodology on more realistic continuous-energy benchmark configurations, such as those issued from the ICSBEP database.