Coupled Probabilistic and Point Kinetics Modelling of Fast Pulses in Nuclear Systems

This paper describes a probabilistic method of modelling point nuclear systems with low numbers of neutrons including the effects of delayed neutron precursors and its coupling with standard point kinetics equations. This coupling allows the simulation of the non-deterministic progression of a system transitioning from subcritical to supercritical and the resulting power peak. Through analysis of large numbers of realisations various statistical parameters of such transients can be obtained. The method of simulation presented here successfully replicates the survival and extinction probabilities predicted by the Backwards Master Equation and experimental and analytic results from the literature regarding the Godiva reactor and extends the examination of that reactor. In particular the effect of delayed neutrons on the simulated response of Godiva is highlighted. With its implementation in a parallel computer code, the model is able to simulate at a reasonable speed a range of systems where low neutron populations are important.


Introduction
In many nuclear systems the instance where the neutron population is small can be of specific interest from operational or safety viewpoints. When the neutron population is very low there is more time for the reactivity to change before the neutron population is able to respond enough to produce a negative feedback as the neutron population must first grow sufficiently large that it is unlikely to become extinct and then it must grow several orders of magnitude to produce a significant feedback. This means the case of a low neutron population is one scenario under which even a system with strong negative reactivity feedbacks may find itself in a situation of a significant neutron population and supercritical reactivity with all the resulting safety implications. Such cases include the start-up of power reactors, accident scenarios in which the reactivity changes over time, such as the Y12 accident (Patton et al., 1958;Zamacinski et al., 2014) and pulse reactors such as Godiva (Hansen, 1960) and Caliban (Humbert and Méchitoua, 2004;Authier et al., 2014;Williams, 2015).
When the neutron population is sufficiently low, the behaviour of the neutron population is inherently non-deterministic as each neutron has a probability of meeting a number of different fates resulting in the removal of the neutron from the population and the potential creation of new neutrons and/or delayed neutron precursors. The resultant behaviour of a system of neutrons will be different on each occasion. An early formulation of the behaviour of such systems may be found in Bell (1963Bell ( , 1965 which examine the probability of a certain number of neutrons and neutron precursors existing at a given time. However, such treatments quickly become computationally demanding in all but the most simplified cases so examination of derived parameters such as the survival probability through the lumped backward Master equation is often employed (Cooling et al., 2016). Such treatments tend to assume a lumped reactor with point variables and delayed neutron precursors are frequently neglected. Treatments which include the effects of delayed neutrons are beginning to appear (Hayes and Allen, 2005;Ha and Kim, 2011;Williams and Pázsit, 2015) but this is still not common.
At low neutron populations a probabilistic method with similarities to the Monte Carlo method as discussed in Booth (2010) and Gang (2011). However, it is the coupling with point kinetics models which lends the present work its novelty. The point kinetics method is a well-known approximation to the behaviour of a reactor where the spatial variation of the neutron flux is not taken into account when calculating the rate of change of the total neutron population (Blue et al., 1964;Lapenta et al., 2001;Cooling et al., 2013). A major assumption of the method is that the number of neutrons and precursors is high enough that the average behaviour calculated is representative of a real system. Even when the neutron population is large, each neutron will still have a finite probability of spawning new neutrons through fission when it is removed from the system which is independent of the fates met by other neutrons in the system. However, with enough neutrons, the relative standard deviation of the ratio of the number of neutrons in one generation to the next becomes small and so the approximation is valid.
The goal of the work presented in this paper is to first develop a model of the change in neutron population which describes the probabilistic behaviour of independent neutrons at low neutron populations and switches to using the point kinetics approximation at higher neutron populations. This allows a single deterministic run to be analogous to a single instance of a neutron population within a given configuration of fissile material. Using this model, a large number of realisations are created and the relevant data from each saved to calculate the desired outputs of the ensemble of realisations.

Low Population Probabilistic Modelling
When the neutron population is low the model treats each neutron individually and at each timestep calculates the probability of a number of reactions occurring to a neutron, precursor, and source in a manner similar to that employed in Gang (2011). The information supplied to the program is the p n (the probability of a fission producing n prompt neutrons), β m (the probability of a given neutron produced in fission being via a delayed neutron precursor in delayed neutron precursor group m), λ m (the decay constants of delayed neutron precursors in group m), Λ (the generation time of neutron) and Σ c , Σ f and Σ esc (the cross-sections for capture, fission and escape). Note that, for most of the examples in this paper these parameters are independent of time or the state of the reactor but these could be arbitrary functions of such variables. For the purposes of this paper we also make the simplifying assumption that the system is close enough to critical that the prompt neutron lifetime is effectively equivalent to the generation time.
The cross-sections represent the relative probabilities of a neutron finally experiencing each fate (capture in non-fissile material, causing a fission, and escaping the system). As such we may write: where p c , p f , p esc represent the absolute probability of a neutron finally meeting each fate (assuming the state of the system doesn't change in such a way that the Σ values are altered).
Thus, we may write the probabilities of a neutron being absorbed, of undergoing a fission, and of escaping the system in a small time interval dt as: and the probability of a neutron surviving the time interval dt without undergoing any interaction as: We proceed to evaluate the probabilities of a neutron causing an individual fission which produces a given number of prompt neutrons and delayed neutron precursors. If a neutron undergoes fission we make the assumption that at most one delayed neutron precursor is created. Data relating to correlations between the number of neutrons produced in a fission and the probability of a fission producing a delayed neutron precursor is not available. In the absence of this data, we also make the assumption that the probability a fission producing a delayed neutron precursor of precursor group m is independent of the number of prompt neutrons produced in a fission. Finally, we assume that β = m β m << 1.
Before proceeding we remind the reader that p f is the probability of a neutron causing a fission in a small time period dt, p n is the probability of a fission producing n prompt neutrons and β m is the fraction of the total number of neutrons produced by fission which will be via the production of delayed neutron precursor in group m. Thus, the probability of neutron causing a fission produces n prompt neutrons in the short time period dt is p f p n . We may also find the probability of a fission also producing a delayed neutron precursor of group m as beingνβ m whereν is the mean number of neutrons produced per fission (whether prompt neutrons or via a delayed neutron precursor) which is defined by:ν =νβ + n p n n, which, after rearranging, gives:ν = n p n n 1 − β .
First we consider the case of the neutron causing a fission including the production of a delayed neutron precursor, we may write: where p f,n,m is the probability of a neutron causing a fission in a short period of time of duration dt which releases n prompt neutrons and a delayed neutron precursor of group m. p f,n is the probability that the neutron causes the release of n neutrons with no delayed neutron precursors in a short period of time of duration dt. p f p n is the probability that a neutron causes a fission which causes the release of n prompt neutrons (regardless of whether any delayed neutron precursors are produced) in this time period. Thus we may sum the probabilities of the neutron causing a fission which results in the creation of all different combinations of prompt neutrons and delayed neutron precursors which involve the creation of n neutrons in total and set it equal to p f p n . This allows us to write that: We may combine Equations (11) and (12) to write: Note thatνβ < 1 as as a consequence of the earlier assumption that β << 1.
We may also write the probability of a delayed neutron precursor of group m decaying in a short period of time of duration dt as: We allow for several independent sources. The ith source releases S i,n neutrons in each decay and has a decay rate S i,s decays per second. The probability of a source producing S i,n neutrons in a short period of time of duration dt is thus: To construct the state of the system at time t + dt from the state at time t the state at t + dt begins by containing no neutrons or delayed neutron precursors. Then, for each prompt neutron in the state at time t a random number is generated and what happens to it in that timestep is decided according to the probabilities p surv , p f,n,m and p f,n . If the survival case is obtained the state at t + dt gains a single neutron. If the case corresponding to a fission which produces n prompt neutrons and a delayed neutron precursor in group m is obtained then n prompt neutrons are added along with a single delayed neutron precursor in group m and if the case of fission without the production of a delayed neutron precursor is selected then n prompt neutrons are added and. For each source a random number between zero and one is generated and if it is below p s,i , S i,n prompt neutrons are added. For each delayed neutron precursor in group m at time t a random number between zero and one is generated. If it is less than p decay,m a prompt neutron is added to the state at time t + dt and, if not, a delayed neutron precursor is added to group m to the state at time t + dt. Through this process the state of the system at time t + dt is constructed from the state of the system at time t.
dt should be chosen such that the probabilities of multiple events occurring to a single neutron/precursor in a timestep are small so that the code will not commonly overlook cases where, for instance, a delayed neutron precursor decays and then the resultant neutron initiates a fission all within a timestep dt. We define a tolerable probability of such a chain of events occurring with a single timestep as p multi . There are three criteria which may determine p multi . First, if there are prompt neutrons in the system at time t: where p surv is defined in terms of dt in Equation (8). If delayed neutron precursors are present then, for each m: where p decay,m is defined in terms of of dt in Equation (14). If a source is present then, for each i: where p s,i is defined in terms of of dt in Equation (15). dt is picked such that it is the largest it can be whilst fulfilling all of whichever of Equations (16), (17) and (18) are appropriate given the number of prompt neutrons, delayed neutron precursors and source at time t. Throughout this paper p multi wil be taken to be 10 −4 which leads to timesteps falling roughly in the region between 10 −7 s and 10 −4 s In this manner the timestep of the model of the reactor used in the case of a low neutron population is adaptive in terms of the timestep which greatly aids the efficiency of the simulation.

High Population Deterministic Model
The model described in Section 2.1 is not efficient once a neutron population gets too high as each neutron and delayed neutron precursor must be treated separately. Instead the neutron and precursor populations are modified according to a simple set of point kinetics equations: where N (t) is the number of prompt neutrons at time t, k eff (t) is the effective multiplication factor at time t and C m (t) is the number of delayed neutron precursors in group m at time t. The following equations are used to calcualte k eff (t): whereν is the average number of neutrons released per fission (including any neutrons released via delayed neutron precursors) and is defined in Equation (10). Note that Equations (21) and (10) are also relevant when the probabilistic model is being used. The equations are solved by a subroutine based upon an algorithm detailed by Shampine and Gordon (1975). The algorithm is designed to solve Ordinary Differential Equations (ODEs) and uses adaptive timestepping. It will not be discussed in detail here as it is not the focus of this paper but it has proved accurate and robust.

Simulating a Realisation
Sections 2.1 and 2.2 describe the two methods used to progress a simulation in time. In the simulations performed here, a realisation will begin with the construction of the initial conditions of the system and using this as the state of the system at t=0. At each time the state of the reactor is examined and either the process described in Section 2.1 or that described in Section 2.2 is used to advance the simulation. The criteria used to chose which of the processes is used is the combined number of prompt neutrons and delayed neutron precursors in the system N tot (t) = N (t) + m C m (t). If this is larger than the threshold value N thresh then the deterministic point kinetics method is used but, if it is lower, the probabilistic method is used. In this fashion a simulation of a realisation can switch method multiple times as t increases in order to use the most appropriate method at any given moment. The choice of a suitable value for N thresh is discussed in Section 5.2.

Multiple Realisations
Section 2 describes how a single realisation may be simulated. However, the results of this code will be non-deterministic due to the probabilistic method used to model the neutron and delayed neutron precursor populations when the combined population is low within each realisation. As a result, a large number of realisations are run so that statistics may be recorded. Effectively, this is a Monte Carlo sampling of probabilistic neutron histories which shares fundamental similarities to Monte Carlo codes such as MCNP (X-5 Monte Carlo Team, 2003), even though the manner in which neutron histories are simulated is significantly different.
As each realisation is entirely independent of any other this process is easy to run in a parallel mode on several cores at once. This is achieved using the MPI parallel coding interface. At the end of each run the values relevant to a particular statistic are recorded on tallies local to each rank. This means it is not necessary to store all the results of each realisation which significantly reduces the memory required. After all realisations have been run the relevant statistics are collected from each rank, collated and the output plotted.

The Example System
In this paper we will examine an abstract system which is not designed to reflect any specific real system, only to demonstrate this method of simulating neutron populations. The values of Σ s , Σ f and Σ c will vary from case to case to vary the criticality of the system as required. The parameters Λ, β i and λ i are taken from the Aqueous Homogeneous Reactor (AHR) described in Cooling et al. (2013) and the values of p n are taken from Zucker and Holden (1986). These parameters are summarised in Table  1 and are common to all simulations in this paper. Again, it is stressed that this example system is only being used to provide parameters which are broadly representative of reality to demonstrate this method of simulation and there is no claim that this system is any way special or even representative of a real or potential system.

Probabilistic Neutronics Sample Results
We begin by showing two typical results of the probabilistic method. The first is a subcritical case with a single source with a strength of 1000 incidents per second with each incident releasing a single prompt neutron (such that S 1,n = 1 and S 1,s = 1000). Both Σ c and Σ esc are equal to 0.3 whilst Σ f =0.4. Through Equation (21) we find this equates to a k eff of 0.970. The results are shown in Figure 1 and are typical of a subcritical system. The prompt neutron population varies on a timescale similar to that of the neutron lifetime and the reciprocal of the source strength. Incidents where the neutron population grows and shrinks can be seen but, given the system is subcritical, it never diverges. The number of delayed neutron precursors varies significantly more slowly as the rate of production is much lower and they are more long lived. In the realisation presented it happens that only two delayed neutron precursors are created: one in group 2 and one in group 4. The next case we observe is that of a supercritical system which uses the same conditions as those used to produce Figure 1 but with Σ c and Σ esc are equal to 0.27 whilst Σ f =0.46 (resulting in a value of k eff of 1.116). Early in the simulation the behaviour of the system is very similar to that found in Figure 1. Whilst the neutron population could have remained low for this finite time period, the population grows in this case and, after a period of time, tends towards an exponential rate of growth with a fairly steady time constant. The numbers of delayed neutrons also begin to increase at a similar pace although few or no decays are seen as the lifetime of the precursors is longer than the simulated period. The ratio of delayed neutron precursor populations to each other varies over time but is roughly proportional to the size of the corresponding β i found in Table 1.

The Choice of N thresh
When choosing a value of n thresh (which was introduced in Section 2.3) the key question is how many neutrons need to be in the system so that the non-deterministic physical system resembles the deterministic point kinetics questions introduced in Section 2.2 with acceptable accuracy. We now examine at what size of neutron population this occurs.
We examine this by creating a number of cases where the number of neutrons varies across a large range and recording the observed reciprocal period of the neutron population as simulated by the probabilistic method. By repeating the simulation a large number of times the mean and the standard deviation of the reciprocal period can be obtained as a function of the total number of neutrons and delayed neutron precursors. The relevant equations are as follows: where N tot (t) is the total number of prompt neutrons and delayed neutron precursors at time t and 1 T (t → t + ∆ t ) is the measured reciprocal period as measured in the time interval t to t + ∆ t . To examine this three simulations were run. The supercritical case from Section 5.1 will be used. Although in each case the time step used for advancing the model dt was held constant at 10 5 s, the value of ∆ t used calculating the observed inverse period was 10 −5 s for the first simulation, 10 −4 s for the second and 10 −3 s for the third. In each simulation, 100,000 realisations were conducted each lasting 1s. As the simulation in Section 5.1 shows, N tot (t) is very likely to diverge in that time and so each realisation will sample a number of points over the range of values of N tot (t). The results are shown in Figure 3.
In all cases the mean value stays approximately constant, showing that the proportional rate of increase of the number of neutrons and precursors was independent of both the number of neutrons in the system and the period over which a measurement occurs. In all cases both the number of times the ensemble of realisations observed the reactor having a certain number of prompt neutrons and precursors decreased as the number of prompt neutrons and precursors increased. This is because, as the expected rate of change is proportional to the number of prompt neutrons and precursors and so a given realisation of a system is more likely to skip a specific number of prompt neutrons and precursors as the number of prompt neutrons and precursors increases.
The standard deviation of the measured reciprocal period also drops significantly as the number of prompt neutrons and precursors increases. This is because the denominator of Equation (23) grows larger with N tot faster than the numerator. The standard deviation also drops more quickly and to a lower level as the value of ∆ t increases. This is due to the fact that the size of ∆ t relative to Λ (1.8202 × 10 −4 s in this case) determines the number of interactions which happen per prompt neutron over the measurement period ∆ t . The larger the number of interactions which occur the more predictable the value of 1 T (t → t + ∆ t ). These results allow us to estimate the number of prompt neutrons and delayed neutron precursors required for the deterministic modelling described in Section 2.2 to be appropriate. The condition is that the standard deviation in the reciprocal period over the temporal resolution of the measurement being taken be significantly lower than the mean. For the remainder of the simulations in this paper the finest temporal resolution used when the deterministic model may be used due to a high number of prompt neutrons and precursor will be around 10 −3 s . As a result we select a value of N thresh of 2000.

Survival Probability
One basic simulation which may be performed is the observation of a single initial neutron and the observation of the probability of the resulting neutron chain existing as a function of time in the absence of any sources. For critical and subcritical cases this will always tend to zero as t → ∞ and will tend to a finite non-zero value for supercritical cases as all chains will eventually die out or attain a high enough population that the probability of it decreasing to zero tends to zero. The two cases studied here will be the subcritical and supercritical cases described in Section 5.1 but with one initial prompt neutron and no sources present. As the reactivity is constant in this case, these results may be compared with results simulated using the Backward Master Equation which is utilised for a similar problem in Cooling et al. (2016). That version of that method is limited to considering prompt neutrons only and matches very well with the corresponding results obtained by using the method described in this paper. Figure 4 shows the results for the subcritical case both with and without the presence of delayed neutrons. The survival probability tends to zero as t → ∞ as all neutrons chains will eventually die out in a subcritical or critical system. The majority of this decrease occurs over the same timescale as the neutron lifetime (1.8202 × 10 −4 s) but the survival probability begins decreasing much more slowly when it reaches about 0.1s. This is because there is a finite chance the neutron chain has produced one or more delayed neutron precursors which decay over the timescale of 1 λ i (typically 0.1s to 100s). It is over these timescales that the final descent of the survival probability to zero occurs. The results share the major features of results of simulations of similar systems presented in Williams and Pázsit (2015). In the case where no delayed neutrons are included (achieved by temporarily setting β i = 0) this secondary descent is absent and, instead, all neutron chains decay on the timescale of the neutron lifetime.  Figure 5 shows the decay of the survival probability in the supercritical case described in Section 5.1 both with and without delayed neutron precursors. In the case without delayed neutron precursors, the k eff value as calculated by Equation (21) is actually slightly lower (1.108) than the case with delayed neutrons (1.116) the survival probability as t → ∞ is also slightly lower (0.1094 compared to 0.1127). Like the subcritical case, the survival probability initially drops with a timescale equal to the prompt lifetime. In the case without delayed neutrons the survival probability, again, tends to its final value immediately whilst, in the case with delayed neutrons, this progression is slowed to the timescale of the lifetime of the delayed neutron precursors as a number of precursors will be produced.

Extinction Probability with the Presence of a Source
Another basic result which may be examined is the probability that no neutrons are present at a particular time in the presence of a source. Two types of "extinction probability" will be considered. In the first the neutron chain will be considered extinct if no prompt neutrons or delayed neutron precursors are present. In the second the chain will be considered extinct when no prompt neutrons are present (regardless of whether delayed neutron precursors are present or not). At t = 0 no neutrons or precursors are present. Again, the subcritical and supercritical cases from Section 5.1 will be used with the same source with a strength of 1000 incidents/s with each incident releasing 1 neutron. We will also examine the case where no delayed precursors can be created (by setting β i to zero in the same fashion as Section 5.3). In this final case no delayed neutron precursors are produced and so the two definitions of extinction probability are identical.  Figure 6a shows the two different definitions of extinction probability and the case without any delayed neutrons for the subcritical case. In both cases the extinction probability begins to decrease on a timescale similar to 1 S 1,s as the first neutron chains are initiated.
In the subcritical case the definition of extinction which requires there to be no delayed neutrons present for the chain to be extinct tends to zero as delayed neutron precursors decay relatively slowly and are produced in a large enough fraction of fissions that there are soon significant numbers of delayed neutron precursors in each realisation (see Figure 7a). For the case where only the lack of prompt neutrons is required for the system to be considered extinct or the case where no delayed neutrons are modelled the behaviour is similar and the extinction probability tends to a value around 0.48. This is consistent with prompt neutrons being created by fission and the source and being removed from the system through capture at a similar rate. In the supercritical case the differing definitions of "extinct" or the inclusion or exclusion of delayed neutron precursors makes little difference and the extinction probability tends to zero. This is because, in every realisation the population of prompt neutrons eventually rises high enough such that it will never return to zero. This can be seen in Figure 7b where the number of delayed neutron precursors continues to grow at a significant rate.
Both the survival probability and extinction probability are special cases of the more general ability of the model to asses the probability of the system having different numbers of prompt neutrons and/or delayed neutron precursors. It is easy to extract from the model the probability of the system having at least, no more than or an exact number of prompt neutrons and/or delayed neutron precursors in a specific group or groups at a given time or as a function of time through examination of each realisation and keeping a tally. This would allow the construction of the type of probability distribution described in Bell (1963) which forms the basis of the Forward equation.

Power Peaks
This model is also capabale of modelling a power peak. To do this the relevant cross-sections are re-defined as follows: where t ramp is the duration of the reactivity ramp which causes the peak and is equal to 1s in this case. γ E,esc is the rate at which the escape probability increases with the total energy released and is equal to 10 −9 J −1 . This is designed to qualitatively approximate negative reactivity feedbacks which are found in many nuclear systems where, over short timescales, the reactivity may be approximated to decrease linearly with the energy deposited. It is noted that these feedbacks themselves may also contain some stochastic effects (such as those associated with the transitions between different boiling regimes in a water-cooled reactor), as these stochastic effects are not the focus of this paper we make the simplifying assumption that the feedbacks of the systems studied are deterministic. E(t) is the energy released by fissions in the system. It is defined as: where χ f = 3.204 × 10 −11 J is the energy released per fission and F (t) is the fission rate. Note that this fission rate can either be a series of delta-functions in the instance where the probabilistic model is being used or a continuous function of time where the deterministic model is being used. In this case it defined as: A singleton source with a strength of 1000n/s is included and there are no neutrons or delayed neutron precursors present at t = 0. All 6 delayed neutron groups described in Table 1 are included. The result of these definitions is that the system will be initially subcritical before the falling capture cross-section causes the system to become supercritical, allowing a neutron chain to grow large. As fissions release energy the escape cross-section will increase and the system will become subcritical and the number of neutrons and the power will decrease; thus, a power peak is created. As the system is a conceptual system used as a proof of concept for this model, the details of what the cross-sections and their variations might physically represent are not of any importance: the overall qualitative behaviour of the k eff in response to the progression of time and the energy released is what is important.

Individual Realisations
We first examine a number of individual realisations in Figure 8. In each case the prompt neutron population is initially fluctuating at a low value in a similar fashion to Figure 1 as the system is subcritical. The value of k eff first goes above one at 0.30s and the system becomes prompt supercritical at t=0.36s. At this point any neutron chain has the potential of increasing significantly in number but, as shown in Figure  7b, this is not guaranteed to happen immediately. In the realisations shown in Figure  8, Realisation 1 increases and reaches its peak first and Realisation 3 increases last. It can be observed that the later a chain starts increasing significantly in number, the higher its peak power and the later it is attained. At the peak power the total energy released increases significantly and the k eff drops sharply, causing the prompt population to decrease. However the definition of Σ c in Equation (24) means that the k eff subsequently rises again until t = 1s. The higher k eff and the decay of delayed neutrons created in the power peak cause a gentle rise in power in this region which is more pronounced the earlier the power peak occurred. The non-zero power causes the k eff to reduce over time and so the prompt neutron population decreases at a rate largely determined by the rate of decay of delayed neutrons. Also displayed in Figure 8 are the results obtained when only the deterministic model is employed (in practice this is achieved by setting N thresh = −1). This result displays a much smoother behaviour whilst the system is subcritical and the population is small. Once the population begins to rise the behaviour is very similar to the coupled realisations. The deterministic case happens to rise slightly earlier than all three realisations but there is no guarantee this will be the case for any coupled realisation.

Ensemble of Realisations
More information on the behaviour of such a system can be obtained through the examination of a large number of realisations. Details of the time-progression of several statistical quantities are displayed in Figure 9. Before the system becomes critical at 0.30s the neutron and delayed neutron precursor populations are low and vary significantly from realisation to realisation due to the non-deterministic behaviour in this physical regime. This causes the standard deviation in the prompt neutron population and the total amount of energy released to be comparable in size to the relevant means. The k eff rises but does not see a significant standard deviation as little enough energy has been released in any realisation that its value remains largely unaffected by the feedback term in Equation (26).
When the system becomes supercritical the mean prompt neutron and delayed neutron precursor populations increase significantly as the neutron population in some realisations begin to diverge. This also causes the standard deviation of the k eff to rise as, in some realisations, significant amounts of energy have now been released whilst, in others, almost none has been released. The standard deviation of the prompt neutron population also increases sharply and remains the same order of magnitude as the mean. This is because, as shown by Figure 8a, in each realisation the prompt neutron population rises to a very high level for a short period of time. This means that, across the ensembles of realisations at this time, a few have very high prompt neutron populations whilst all the others have significantly lower populations either as they are yet to have their power peak or because they already have and have proceeded to the lower plateau after the peak.
Eventually, by around 0.9s almost all chains have experienced their power peak and, although the k eff still rises until t=1s due to Equation (24), it remains subcritical for the remainder of the simulation. The prompt neutron and delayed neutron precursor populations begin to decrease relatively slowly as delayed neutron precursors formed in the power peak decay, releasing prompt neutrons. At this time the standard deviation of the prompt neutron population, the total energy released and the k eff drops so it is significantly lower than the mean. This is because, regardless of when the power peak occurs a similar total amount of energy has been released at approximately the same time, meaning there are approximately the same number of delayed neutron precursors to decay in an environment with a similar k eff . The deterministic results appear qualitatively similar to the mean results. The peak power of the deterministic case is actually a little larger and is obtained a little earlier than the peak of the mean power. This is because the peak value of the average neutron population is an average of several sharp peaks and appears when these peaks are most tightly grouped but is less sharply peaked than any of the individual peaks.  We may gain more information by looking at the ensemble of results of the distribution of peak powers and total energy released over the course of the 100s of the simulation, which are displayed in Figure 10. The peak power shows a distribution with a mean of 3.62×10 9 W and a standard deviation of 4.89×10 8 W. A small minority of the realisation showed a peak power of over 5×10 9 W and the largest power observed was 7.51×10 9 W. However, this is not a hard upper limit to the actual distribution as the distribution shows that there are relatively few realisations which attain this power (hence the jagged appearance of the distribution). This means the scenario which theoretically attains the highest possible peak power was simply not simulated in the 10,000 realisations performed. Nonetheless, the distribution here clearly shows the distribution which will be displayed in the vast majority of instances.
The total power released is much more tightly grouped. The long tail of low probability events with higher values is still present, but it does not extend as far compared to the mean as the peak power. The mean energy released is 9.96×10 7 J whilst the standard deviation is only 1.90×10 6 J. This means that the amount of energy produced by the system in this transient is fairly consistent.
The values of the deterministic results are also included in Figure 10. As can be seen, the deterministic peak power and total energy released are 3.27×10 9 W and 9.84×10 7 J respectively and are slightly lower than the mean values produced by the coupled model by 9.7% and 0.1% respectively.

Comparing Power Peaks
We may simulate variations of the system introduced at the start of Section 5.5 and compare the distributions of peak power and total energy produced. The two parameters we will chose to vary here are the length of time over which the reactivity ramp occurs t ramp (see Equation (5.5)) and the strength S 1,s of the single neutron emission source. The results are shown in Figure 11. tramp=1s, S1,s=1000/s tramp=2s, S1,s=1000/s tramp=10s, S1,s=1000/s tramp=1s, S1,s=500/s tramp=1s, S1,s=10000/s 1.0e+08 1.1e+08 1.2e+08 Probability Density (1/J) Energy Released (J) at t=100s tramp=1s, S1,s=1000/s tramp=2s, S1,s=1000/s tramp=10s, S1,s=1000/s tramp=1s, S1,s=500/s tramp=1s, S1,s=10000/s (b) Energy Released at t=100s For the same value of t ramp , the source intensity does not affect the lowest peak power or lowest total energy released in a realisation. This is because both of these events occur when a neutron chain is initiated just as the reactor becomes prompt supercritical (meaning the transient occurs when the reactivity is as low as possible) and this is possible for any non-zero source strength. However, for higher source strengths the distributions of the peak power and total energy released are more tightly grouped to these lowest values, resulting in lower means and standard deviations for these values. This is because the higher the source strength the more neutrons are likely to be released when the reactivity is only just prompt supercritical and the more neutron chains have a chance of growing and forming the power peak when the reactivity is low. This results in a lower peak power and total energy released as the amount of energy needed to be released to cause the system to be subcritical again. Additionally, the lower the reactivity when the neutron chain diverges the slower the increase in power as a function of time will be and so the lower the peak power will be.
For the same source intensity, the higher the value of t ramp the lower the peak power. This is for two reasons. Firstly, in a fashion similar to the higher source strength, the higher the value of t ramp the slower the reactivity increases, meaning there are more neutron chains initiated when the reactivity is relatively low, increasing the probability that a neutron chain will diverge and form the power peak when the reactivity is low. The second effect occurs because the rise in power occurs over an appreciable amount of time. For instance, for realisation 3 in Figure 8 the power begins to increase at t = 0.60s and the peak power occurs at t = 0.80s. For a smaller value of t ramp the reactivity will increase more while the power is rising, meaning the power will need to rise higher and more energy will need to be released to cause the reactor to become subcritical.

Godiva
Godiva is an unreflected sphere of Highly Enriched Uranium (HEU) which was dominated by fast fission and was operated in a pulsed operation mode (Peterson, 1953). Hansen (1960) presents a substantial amount of theoretical and practical work surrounding the Godiva reactor which has become somewhat of a benchmark in the analysis of systems with a low neutron population. In this work we attempt to recreate some of the results present in Hansen's work both as an act of validation for this work and as an attempt to present new evidence and explanations to a historic problem. This section attempts to recreate two key experiments discussed in Hansen (1960) which involve a ramp insertion and a delayed supercritical step insertion.

Ramp Insertion
Hansen (1960) describes a scenario under which the reactivity of Godiva is steadily increased in the presence of a weak neutron source, resulting in a power peak at some non-deterministic time and power. We wish to model only the first peak and, as a result, when the total energy released E(t) exceeds the value E thresh (indicating the power peak has begun) the representation of Godiva will change to allow the power peak to end naturally and prevent further neutron chains being started. We adopt a value of 1kJ for E thresh . This change is partially for performance as fewer neutrons will need to be modelled and partially to prevent other neutron chains releasing significant amounts of energy after the initial power.
Hansen made several approximations to the physics in the modelling of Godiva. As a result, two approximations to the Godiva reactor will be presented here: one which attempts to replicate Hansen's model of Godiva and another which attempts to model Godiva as closely as possible. The two main differences are the treatment of delayed neutrons and the treatment of the source.
There are several parameters which are common to both representations of Godiva. Hansen states the prompt neutron lifetime of Godiva was 6 × 10 −9 s and this value will be adopted as Λ in this work. Hansen doesn't state the relevant values of p n for Godiva but, as Godiva was comprised of HEU we approximate it as being that of pure 235 U which is given in Zucker and Holden (1986) and are used in the example system examined in Section 4 and are the same values of p n used in that section which are detailed in Table 1. Hansen stated Godiva had a well-known quenching constant of 5 × 10 −20 per neutron absorbed. Hansen (1960) describes the case where the reactivity increased at a rate between 10 −2 s −1 and 10 −1 s −1 so we choose to use a value of 3×10 −2 s −2 in order to be in the middle of this range. We may combine this with the quenching coefficient of 5 × 10 −20 per neutron absorbed to produce the following expressions for cross-sections (recalling that it is only the relative size of Σ esc + Σ c to Σ f that is important: Σ c (t = 0) = 1.4105 for simulations without delayed neutrons 1.4262 for simulations with delayed neutrons (31) where This formulation allows the reactivity to stop increasing after the initial power peak, meaning the same neutron chain cannot reignite after the power peak as the system remains subcritical. Note that there are two different values for Σ c where which one is used depends on if the model used contains delayed neutron precursors or not. The reason for this is because, as shown in Equation (10), the inclusion of delayed neutron precursors increases the mean number of neutrons produced per fission which, in turn through Equation (21), results in a higher k eff because delayed neutrons are being produced in addition to prompt neutrons. Thus, having two values for Σ c ensures that the system has is exactly critical at t = 0 in both cases. 1.14 6 27.3×10 −5 3.01 Table 2: Delayed neutron precursor group data for the refined model of Godiva. The total value of β is 650×10 −5 . These are the values for pure 235 U given in Lamarsh and Baratta (2001) To create the more refined representation of Godiva six delayed neutron precursor groups are included with the values for each presented in Table 2. Hansen (1960) states the source strength within Godiva as 90 neutrons/s. When attempting to replicate Hansen's representation of Godiva this will be treated as a single source with S 1,s =90s −1 and S 1,n = 1. However, when trying to represent Godiva more exactly we use the fact that we know the fission source in Godiva was largely spontaneous fissions to represent the sources with different values of S i,n as a number of different sources as in Table 3. For the time being, in both cases the source will be switched off once E(t) > E thresh which prevents other neutron chains from starting after the main peak. From here on we refer to the representation of Godiva with a fission source and the presence of six delayed neutron groups as the "modified representation". i S i,n S i,s (s −1 ) 1 1 6.441 2 2 12.550 3 3 11.345 4 4 4.727 5 5 0.978 6 6 0.0977 7 7 0.00530  The results showing the distribution of the energy produced under the different assumptions are shown in Figure 12. A major feature of this figure is that the modified representation produces a significantly larger amount of energy compared to the representation designed to replicate the representation in Hansen. This is primarily due to the the inclusion of delayed neutrons. Without delayed neutrons it is possible that, almost immediately a neutron will be produced and lead to a chain which causes a power peak. Such a peak will be small because little reactivity has been added by the ramp. In the case with delayed neutrons, the system must become prompt supercritical in order to for the power to increase significantly. Once the system drops back to being delayed subcritical the prompt neutron population drops more slowly as a large number of delayed neutron precursors are decaying. This causes the power peak to be much wider and so, particularly for the less energetic realisations, much more energy is released.
Another key observation is that deterministic value for a particular representation falls at the very lowest values of the total amount of energy released as simulated by the coupled model which indicates the deterministic model severely underestimates the mean energy released in the transient. This is one of the key results of Hansen's work. Hansen derives an expression for the ratio of the mean energy released in a non-deterministic model to the mean released in a deterministic model: whereĒ is the mean energy release predicted by by the stochastic system, E k is the energy release predicted by the deterministic system, χ 2 is the second factorial moment of the neutron multiplicity distribution, S is the source strength, a is the reactivity ramp rate and b is the quenching parameter. Hansen obtains a value for this ratio of between 235 and 252. The corresponding value for this work comes from a comparison of the two Hansen representations in Figure 12 which provides a value for this ratio of 230 which is in good agreement with Hansen's work. However, we may then use the data in Figure 12 to perform the same analysis for the modified representation of Godiva (which includes the fission source and delayed neutrons) for which we obtain a value of 4.21. This indicates that the effect of delayed neutrons significantly reduces the impact of the non-deterministic behaviour early in the simulation, which deserves further study.
The reason this ratio changes so much with the introduction of delayed neutrons is as follows. In the case where there are no delayed neutrons modelled, as soon as k eff exceeds 1 the number of prompt neutrons may begin to increase under the coupled model and will increase under the deterministic model. The number of prompt neutrons will increase rapidly until the energy released is sufficient to counter the positive reactivity insertion. This may occur when the system has a reactivity of, for instance, 10pcm above critical. This requires a small release of energy to counter. This results in the deterministic case and the least energetic of the coupled simulations to have a similar low energy produced over the course of the power peak. For the coupled case it is also possible that the prompt neutron population may not start to rise until the reactivity is at a few hundred pcm, which would result in a much larger amount of energy released. This description applies to the lines marked "Hansen Representation" in Figure 12 and explains why the "Hansen Presentation, Coupled" case produces a large range of energies produced.
With delayed neutrons however, the system must reach the state of being prompt supercritical (that is the system is 650pcm above critical in this case). At this point the number of prompt neutrons may begin to increase under the coupled model and will increase under the deterministic model. The number of prompt neutrons will this increase until the energy released causes the system to drop back to being delayed supercritical. However, many delayed neutron precursors were created in this prompt peak and they will cause the power to remain at a non-negligible level until the amount of energy released is enough to counter the ∼650pcm positive insertion from the ramp. This description applies to the lines marked "Modified Representation" in Figure 12. As a result the smallest amount of energy released as predicted by the coupled model and the amount of energy released predicted by the deterministic model increase significantly compared to the case where no delayed neutron precursors were present. For the illustrative numbers used in this and the previous paragraphs this increase is a factor of ∼70. In fact we see the value for the total energy released produced by the deterministic codes increase from 79.5kJ for the Hansen representation to 5.09MJ (a factor of 64.1) and the lowest energy observed to be released from in coupled case produced by the coupled model increase from 70.4kJ for the Hansen representation to 5.08MJ (a factor of 72.2). This effect explains why the low energy tail of the distribution present for the Hansen representation is not present for the Modified representation in Figure  12. However, it is only the low energy tail of the distribution which increases by this factor: the mean changes from 18.3MJ to 21.5MJ (a factor of 1.17). This is the reason that the ratios of the energy produced are so different in the cases when delayed neutrons are neglected and included: the introduction of delayed neutrons increases the energy released in the deterministic case much more than it increases the mean energy released in the coupled case.
Another variation in this scenario which may be examined is where neither the source nor ramp insertion of reactivity terminates when the first power peak occurs. When this is done the results are as in Figure 13. The difference in behaviour between the different methods of simulations in marked. In the deterministic version of the case which replicates Hansen's model the power oscillates very rapidly due to the low prompt neutron lifetime. For this reason Figure 13a) plots the average power over the past 1ms rather than the instantaneous power -it is not practical to plot outputs at a sufficiently large number of times to resolve the oscillation fully. In this instance, as soon as the reactor is supercritical the power rises almost instantaneously and the reactivity and power fall. As a result the cumulative energy released appears to rise smoothly and the k eff remains very close to one.
The coupled model simulation of Hansen's representation of the system sees a small number of larger distinct peaks. In this case the k eff can be seen to rise steadily until a neutron produced by a source happens to instigate a neutron chain which grows to a sharp power peak rather than dies out. This makes the system subcritical and the k eff rises steadily until another power peak occurs. As a result the k eff takes the appearance of an irregular sawtooth and the cumulative energy released rises in sharp steps.
The modified representation follows a similar pattern whether the deterministic or coupled method is used. Initially the k eff rises steadily until there is an initial power peak. In the deterministic case this is when the system becomes prompt supercritical and in the coupled case this is when a neutron released by a source when the system is prompt supercritical happens to start a chain which grows instead of dying out. Once this has happened however, both see a decrease in power. The power does not drop as much as when Hansen's representation is used as the delayed neutron precursors created in the power peak continue to decay and prevent the population falling far. In the coupled case the power falls to a trough before rising again as more negative reactivity was inserted in the initial power peak. However, both cases tend towards the same value in terms of both current power and cumulative power released. The power slowly drops in this region as the excess delayed neutron precursors created in the initial power peak slowly decay.
The distribution of the total energy released in Figure 13d is much tighter than those of Figure 12. This is because the total reactivity insertion over 10s is the same in all cases, unlike the previous case where the ramp reactivity stopped when the energy released reached a certain level or, as Hansen modelled the system, at the first peak power. As a result the total amount of energy released is very close to the amount of energy required to provide enough negative feedback to counter the positive feedback inserted.
For the deterministic functions the amount of energy released is, again, a delta function. The amounts of energy released in the deterministic cases (for the Hansen and the modified representations) are very similar. The amount of energy released in the simulation performed with the coupled model of the modified system is not a delta function but a very narrow distribution with a standard deviation of 10.5kJ and a mean of 82.780MJ that lies almost exactly on the value produced by the modified deterministic model (82.825MJ). This is because the system in both cases, after an initial power peak, the delayed neutrons cause the system to tend smoothly to a state where the total amount of energy released exactly provides sufficient negative reactivity to exactly counter the positive reactivity inserted into the system. It is only the case where the coupled model is used with Hansen's original model of Godiva (which neglected delayed neutron precursors) that sees a significant distribution in the total amount of energy released and even this is fairly tightly grouped (it has a standard deviation of 6.67MJ compared to a mean of 79.572MJ which is very close to the deterministic value of 79.752MJ). This distribution occurs because, for a given realisation at t=10s he system may have just undergone a power peak and so have produced more energy than required to make the system subcritical or the system may currently be supercritical as not enough energy has been released to make it subcritical (Figure 13c illustrates this well).
We may again compute the key result from the work of Hansen (1960): the ratio of the mean energy released in the stochastic case to the energy released in the deterministic case. For both cases (Hansen's representation and the modified representation with delayed neutrons) this value is almost exactly one.     A summary of the ratios of the average energy released predicted by a non-deterministic model to the energy released predicted by a deterministic model is displayed in Table  4. These results imply that the reason Hansen obtained large values for this ratio were the assumptions he made regarding the physics of the system and the value he chose to simulate. Specifically, Hansen neglected the effects of delayed neutrons, only the energy released to the first peak was used the value of interest and both the ramp reactivity insertion and source were considered to be switched off when the peak power is obtained. If the effects of delayed neutrons are accounted for this value drops significantly and if the value of interest is the energy released following a fixed period of time (which contains the first peak) instead of to the first peak only and if the source and reactivity ramp remain after the first power peak then this ratio drops to one. Examination of Figure 13 clearly shows that different sets of assumption can still have a very marked effect on the behaviour and the progression of the system with time even if the specific ratio discussed here is close to one.

Step Insertion
Hansen (1960) describes another set of experiments performed with Godiva: those of a step insertion of reactivity such that, after the step k eff − 1 = 0.0047. The result is a delayed supercritical insertion with a slowly increasing power. The experiment was repeated 30 times and, on each occasion, the amount of time taken for the prompt neutron population to reach 4200 neutrons was observed. Hansen presents a model containing one delayed neutron precursor group which approximates the mean and standard distribution of this time well. We now attempt to replicate this experiment with this model. Unfortunately Hansen does not state the decay constant of the one delayed neutron precursor group he used in his model so it is not possible to recreate his representation of Godiva in this case. Instead, we rely on the modified model of Godiva used in Section 6.1 with the following definitions of the cross-sections: Again we use a generation time of 6 × 10 −9 s. For each realisation the time t when n prompt is first at least 4200 neutrons is recorded. The distribution for this value is presented in Figure 14. The mean and standard deviation obtained are 29.8s and 5.2s respectively, compared to the experimental values presented in Hansen of 31.8s and 4.6s and the results of Hansen's model which were 32.5s and 3.3s. All of these values are fairly close to each other indicating good agreement and helping to validate this model.

Model Performance
The simulations presented in this paper were all run on a 20 core desktop. To run a simulation of 10,000 realisations such as that in Figure 5.5 typically took around one hour (although this could vary significantly with the problem being modelled). The amount of memory required is minimal as each rank only needs to store the information regarding the realisation it is currently working on and the relevant tallies and statistics for the desired observations before the information from all ranks is combined after all realisations are complete. As there is almost no communication between ranks it is anticipated that the parallel efficiency would remain very high even when simulations are performed on computers with a very large numbers of cores.  There are different regimes under which the model may operate and the rate at which the model progresses through simulated time will be different in each case. A rough sketch of these regimes is shown in Figure 15. This diagram is very approximate and is used to demonstrate broad behaviour. It neglects changes in k eff and source strength over time and neglects the effects of delayed neutrons.
The first requirement for the code to run slowly is the frequent need for small timesteps (see Section 2.1). This is due to the frequent presence of prompt neutrons or the presence of a high intensity source. The presence of delayed neutron precursors does not often cause very small timesteps because even the shortest-lived precursors have a lifetime of ∼ 0.1s, which is relatively long compared to the lifetime of prompt neutrons or 1 S i,s for a large S i,s . For the system to be inefficient this must also be coupled with a high combined number of prompt neutrons and delayed neutron precursors as these must be simulated at each time step (bearing in mind that this combined population may not exceed N thresh before the more efficient deterministic model takes over).
In regime 1 of Figure 15 the system is supercritical and the neutron population will be high (note that the system will not exist in this state for long before the population becomes very large and overflows in the code, so care must be taken if this regime is part of a simulation). In this case the deterministic model described in Section 2.2 will be employed which is very efficient. In regime 2 the system is subcritical but the source strength is sufficiently high that the equilibrium population is above N thresh and the deterministic model is again employed.
In regime 3 the prompt neutron population is frequently zero and, when it is nonzero, will not frequently be high. This means that the probabilistic model described in Section 2.1 is employed but the time step will be long much of the time and, when it is not, there are not many neutrons to simulate and so the model is efficient in this regime. In regime 4 the probabilistic model is again employed but, this time, the neutron population is frequently non-zero. However, the neutron population does not rise very high. This means that, whilst the time step is short, there are not many neutrons to simulate and so the code runs at a reasonable speed. In regime 5 the probabilistic model is still employed but there may frequently be many hundreds or thousands of neutrons which need to be simulated. This is the regime in which the code runs slowest as the time step is short and there are many neutrons. This is progressively worse the higher the source and the k eff until the point when the neutron population exceeds N thresh . Thus, by lowering N thresh this regime can be minimised, although this comes at the cost of losing the simulated non-deterministic behaviour. This is why studies like that in Section 5.2 are useful for informing this compromise. In regime 6 the source strength is high enough that the time step is forced to be very short which will make the code slow to progress through simulated time.
A low value of Λ may also degrade the efficiency of the code if there are frequently neutrons present and the system is subcritical. Although each neutron chain lasts a smaller amount of time as the generation time is smaller, a sufficiently high source could mean the time step is frequently very small indeed. This limits the practical application of the code for fast fission systems.
It is worth stressing that behaviour of this nature, although they may reduce the efficiency of the code, do not mean the code will not work: it will just take longer to run. A user should be aware of this behaviour when using this model so they may understand its practical limitations and chose parameters for the simulation of a given scenario effectively.

Conclusions
This method has proved to be a very feasible way of modelling a system with a low neutron population and its transition into a critical state with a much higher population. The method is simple to implement and runs reasonably quickly in a large subset of relevant scenarios due in large part to its straightforward implementation in a parallel fashion. The fact that each realisation is performed individually makes extracting meta-data an easy process. This paper has also demonstrated how the low neutron population model may be integrated with traditional deterministic models.
The model shows the impact of delayed neutrons on standard tests of populations with low neutron populations such as survival probability and extinction probability. The largest impact is on the rate at which the system tends towards its final value for a particular measure although the final value itself is generally unchanged. However, due to the long lifetimes of delayed neutron precursors, the final value of the extinction probability is changed significantly if it is considered that no delayed neutron precursors are present for a system to be extinct. The presence of delayed neutron precursors was also shown to increase the minimum energy released in a peak formed by a ramp insertion when the coupled model is applied and also the energy released when the deterministic model is applied.
The coupling of the two simulation methodologies has also allowed the effect obtained when a supercritical system with a low population takes a non-deterministic amount of time to diverge in terms of neutron population and how this can alter the peak power and total energy released if the reactivity is changing over time. The effects of different source strengths and different rates of change of reactivity was also investigated with higher source strengths and slower rates of changes or reactivity being observed to produce more predictable power excursions which release less energy and have a lower peak power.
The model has successfully replicated results predicted by other methods such as the survival probability of a single neutron chain, the extinction probability in a system with a neutron source and experimental data and analytically by Hansen (1960) regarding the Godiva reactor. By modifying the Hansen's representation of Godiva this work has demonstrated how some results are drastically altered by including delayed neutrons and how considering slightly different values of interest can produce very different results.