TOWARDS ZERO-VARIANCE SCHEMES FOR KINETIC MONTE-CARLO SIMULATIONS

The solution of the time-dependent transport problem for neutrons and precursors in a nuclear reactor is hard to treat in a naive Monte-Carlo framework because of the largely different time scales associated to the prompt-fission chains and to the decay of precursors. The increasing computer power and the development of variance-reduction techniques specific for reactor kinetics have recently unlocked the possibility to calculate reference solutions to the time-dependent transport problem. However, the application of timedependent Monte Carlo to large systems (i.e., a full reactor core) is still stifled by the enormous computational requirements. In this paper, we formulate the construction of an optimal Monte-Carlo strategy (in the sense that it results in a zero-variance estimator) for a specific observable in time-dependent transport, in analogy with the existing schemes for stationary problems. As far as we are aware, zero-variance Monte-Carlo schemes for neutron-precursor kinetics have never been proposed before. We verify our construction with numerical calculations for a benchmark transport problem.


INTRODUCTION
The time-dependent form of the neutron transport equation in the presence of delayed neutrons is generally difficult to solve with Monte-Carlo methods because the system is stiff, i.e. the time scale of the prompt fission chains is much smaller than the time scale of precursor decay (by a factor of about 10 −4 in a typical pressurized-water reactor). It is only in the last ten years that these problems have become tractable, by virtue of the groundbreaking work of Sjenitzer and Hoogenboom [1,2] and the subsequent emergence of a new class of variance-reduction techniques (kinetic Monte Carlo). However, the computation cost of kinetic Monte-Carlo simulations is still prohibitively high, which limits the size of the systems that can be treated in a reasonable time (a cluster of fuel assemblies at most [3]). The extension of kinetic Monte Carlo to larger systems is instrumental for the validation of kinetic deterministic solvers, with or without the simulation of feedbacks (such as thermal-hydraulics or thermomechanics), but likely requires more efficient variance-reduction techniques than those available at the present time.
This paper introduces zero-variance schemes for kinetic Monte Carlo. Our construction is intimately related to the formalism for zero-variance schemes in stationary problems, which have been extensively reviewed by Lux and Koblinger [4] and more recently revisited by Hoogenboom [5].
Our work might inspire approximate zero-variance schemes, such as those relying on deterministic solvers for guidance, which have the potential to unlock enormous gains in computing time for kinetic calculations.

FRAMEWORK AND NOTATION
The natural framework for Monte Carlo is the integral formulation of the neutron transport equations. First, let P = (t, r,Ω, E) denote the coordinates of a generic point in the extended phase space (including time). The integral form of the transport equations reads Here ψ(P ) and χ(P ) respectively represent the collision and emission density in P , and Q(P ) is the source term. The flight operator T reads Tg(P ) = T (P , P )g(P ) dP , where the integration over dP is shorthand for integration over all the coordinates of the extended phase space: The flight kernel reads the cross sections being assumed to be independent of time. Similarly, the collision operator C reads where the kernel C(P , P ) consists of a prompt and a delayed term: with Notation for the flight and collision operators is customary; the quantities f s , χ p and χ j d respectively stand for the scattering kernel, the prompt neutron spectrum and the delayed neutron spectrum; the summation extends over the precursor families (for the sake of simplicity we consider a single fissile isotope). We will also need the adjoint transport equations, which we write as The adjoint operators involve integration over the kernels of the direct operators, Eqs. (1) and (2): T † g(P ) = T (P, P )g(P ) dP , C † g(P ) = C(P, P )g(P ) dP .
We recall that the adjoint source η ψ (P ) associated to the adjoint collision density ψ † (P ) can be interpreted as a detector response function, because R = η ψ (P )ψ(P ) dP is the expected response.
We remark in passing that, up to the additional integration in dt, the transport equations are formally identical to their stationary counterparts [5]. Therefore, a number of useful results for the stationary case are easily transposed to the time-dependent case.
We restrict our attention to flight-based estimators, i.e. estimators that depend only the initial and final point of the neutron flight * . Such estimators are defined by a function f (P, P ) that specifies the quantity to score when a flight from P to P is sampled. Let M i (P ) be the i-th moment of the score generated by a particle starting from P with unit weight. The estimator is required to be partially unbiased: This condition guarantees that f will yield an unbiased estimate of the contribution associated to any particular flight. Thus, if f is partially unbiased, it is also unbiased, in the sense that it yields an unbiased estimate of the response associated with a particle starting with unit weight from any point P in phase space. The unbiasedness condition is expressed as or, equivalently, which are easily verified to be corollaries of Eqs. (5), (6) and (7).

ZERO-VARIANCE SCHEMES FOR KINETICS WITH DELAYED NEUTRONS
Zero-variance schemes for neutron transport have been extensively studied by Lux and Koblinger [4] and later revisited by Hoogenboom [5]. Most of the attention has been devoted to the study of the stationary transport equation, with a few exceptions [6]. To the best of our knowledge, this paper presents the first formulation of a zero-variance scheme for the solution of the timedependent neutron transport equation including delayed neutrons.
The variance-reduction techniques that have been recently proposed in the literature [1,3,7] rely on an explicit representation of the precursor population. In this work, instead, we choose to work with the condensed representation of delayed fission, Eq. (4). Moreover, we shall limit our analysis to the case of branchless sampling of the collision [1,4]. It is actually possible to construct zerovariance games for schemes that include an explicit representation of the precursor population. Furthermore, although branching schemes are usually less efficient than branchless schemes, the construction of a zero-variance scheme with branching collisions is possible. For the sake of conciseness, however, these discussions have to be omitted and reserved for a future publication.

Construction of the Scheme
We follow the strategy of Ref. 4 in this section. Our goal is to construct a zero-variance game for the estimator f . LetT (P , P ) andĈ(P , P ) be the hitherto unknown kernels of the flight and collision operatorsT andĈ of the sought zero-variance game. We wish to relate them to the original kernels T (P , P ) and C(P , P ), Eqs. (1) and (2). First, we impose that the zero-variance kernels be normalized to unity: In addition, we need to specify how the statistical weight of the neutron must be adjusted in the zero-variance game. The rules for the weight corrections are fixed by the unbiasedness condition, Eq. (8), and stipulate that a flight P → P sampled fromT (P , P ) must be associated with a weight correction w T (P , P ) = T (P , P )/T (P , P ), and a collision P → P sampled fromĈ(P , P ) must be associated with a weight correction w C (P , P ) = C(P , P )/Ĉ(P , P ).
Finally, we impose the zero-variance condition whereM i (P ) is the i-th moment of the score due to a particle starting with unit weight from a generic point P in the random walk defined byT andĈ. By virtue of the unbiasedness condition, Eq. (8), the first moment satisfiesM 1 (P ) = M 1 (P ) = χ † (P ). The second moment satisfies the backward transport equation M 2 (P ) = T (P, P ) f 2 (P, P ) w 2 T (P, P ) dP + 2 T (P, P ) f (P, P ) w 2 T (P, P ) Ĉ (P , P ) M 1 (P ) w C (P , P ) dP dP + T (P, P ) w 2 T (P, P ) Ĉ (P , P ) w 2 C (P , P )M 2 (P ) dP dP . (12) Proceedings of the PHYSOR 2020, Cambridge, United Kingdom In deriving the form of Eq. (12), we have assumed that the collisions of the zero-variance game are sampled with a branchless strategy [1]. The construction of the zero-variance kernelsT (P , P ) andĈ(P , P ) is similar to the one carried out in Ref. 4. We omit it for the sake of brevity and we limit ourselves to stating the final results: C(P , P ) = C(P , P ) χ † (P ) It can be proved [4] that the random walk defined by Eqs. (13) and (14) for flights and collisions, respectively, results in a score of R, as expected, with zero variance, provided that the starting particles are sampled from the modified source distributionQ(P ) = Q(P )χ † (P )/R.

NUMERICAL RESULTS
As a verification, we test our zero-variance scheme in an extremely simplified, yet non-trivial, model: we assume one-speed transport in an infinite medium, so that time is the only relevant phase-space coordinate. Zero-variance schemes for more complex models will be presented in a future publication.

One-speed, Infinite-Medium Model
Since time is the only non-trivial coordinate, we write t in place of P where necessary and we omit the other coordinates. The explicit forms of the kernels (Eqs. (1)-(4)) for the one-speed, infinite-medium model are where H represents the Heaviside function. We limit ourselves to one precursor family. The target response is taken to be the total flux integrated over the time interval t ∈ [0, T ], which corresponds to the adjoint source . Let α ± be the α-eigenvalues associated to the Boltzmann equation for our system, is the prompt neutron decay constant and α + is the inverse reactor period. The solution of the adjoint transport equations, Eqs. (5) and (6), reads up to a normalization constant. The form of Eq. (15) is valid in the general case, i.e. when α + and α − are both different from zero. In the particular case where one of the eigenvalues vanishes, the solution has a slightly different functional form, which we omit.
The definition of the source term Q(t) requires a short discussion. If the source is intended to represent a single neutron entering the system at t = 0, then obviously Q(t) = δ(t). However, if the system is critical at t = 0 and the source term represents the stationary state of the system, then neutrons need to be injected at t = 0 and at later times. Physically, the neutrons entering the system at t > 0 are delayed neutrons produced by the decay of precursors present in the system at t = 0. In this case, the source term reads where ζ = λ/vν d Σ f is the ratio between the number of neutrons and precursors present in the system at t = 0 [1].

Zero-variance Scheme for the Collision Estimator
We present the zero-variance scheme for the collision estimator, f (t , t) = η ψ (t) = Σ −1 t . It can be shown that the zero-variance flight kernel, Eq. (13), takes the form of a non-homogeneous exponential distribution, viz.
where the modified cross sectionΣ is given bŷ ∂t .
The weight correction associated with the flight, Eq. (9), is given by We choose to sample the zero-variance flight kernel by rejection. SinceΣ(t) ≥ Σ t , it is easy to verify that a majorant for the flight kernel is given by which has an exponential dependence on t and is therefore easily sampled from. Note that t is constrained to lie in the [t , T ] interval.
For the zero-variance collision kernel, it is convenient to introduce the precursor importance χ † c (t), which is the solution of the adjoint equation which can be interpreted as the importance of delayed neutrons emitted by a given fission event.
With this proviso, Eq. (14) takes the form and the associated weight correction (Eq. (10)) reads which bears a striking similarity to the weight correction for branchless collisions for the analog kernel [1]. In a branchless interpretation of the collision kernel, the probabilities to select the prompt or delayed channels are respectively given by Note that there is no practical difference between scattering and prompt fission in the one-speed model. If the delayed channel is selected, the time at which the delayed neutron appears is sampled by rejection from the majorant (with t ∈ [t , T ]) As customary, the score s that is accumulated at each collision is s = w/Σ t , where w is the current weight of the neutron. Such zero-variance random walk would in principle live forever, since the particle is never allowed to undergo capture and cannot escape from the domain. In practice, a Russian roulette is needed in order to keep the simulation time finite. When the neutron weight drops below the Russian roulette threshold w RR , it is killed with probability w RR ; if the neutron survives the Russian roulette, its weight is set to 1. The ideal zero-variance behaviour is attained in the limit of a vanishing threshold w RR .

Monte-Carlo simulation
A numerical demonstration of the fact that the scheme realizes a zero-variance estimate of the total neutron flux is presented in Table 1. We fix the values of the physical parameters (cross sections, etc.) and we perform calculations with 10 6 independent histories for different values of the Russian-roulette threshold. Table 1 shows that, as the roulette threshold w RR decreases, the estimated mean value approaches the analytical solution of the problem. At the same time, the estimated standard error decreases approximately like the square root of the threshold. These findings suggest that, barring numerical difficulties, the zero-variance limit would be attained in Table 1: Estimated mean and standard error for the total neutron flux integrated over 0 ≤ t ≤ T , for different values of the Russian-roulette threshold w RR and in different regimes (subcritical, critical and supercritical), with a sample of 10 6 neutrons.

CONCLUSIONS
We have sketched the theoretical construction of a zero-variance scheme for kinetic Monte Carlo including delayed neutrons. The formalism is a straightforward generalization of the zero-variance schemes for stationary calculations. We have verified the correctness of our scheme on a simple one-speed, infinite-medium transport problem: in the limit of vanishing Russian-roulette threshold, the estimated Monte-Carlo mean approaches the analytical solution and the estimated standard error approaches zero. Generalizations to other families of estimators (e.g. track-length) and more realistic transport problems will be presented in a forthcoming publication.