CALLISTO-SPK: A Stochastic Point Kinetics code for performing low source nuclear power plant start-up and power ascension calculations

This paper presents the theory and application of a code called CALLISTO which is used for performing NPP start-up and power ascension calculations. The CALLISTO code is designed to calculate various values relating to the neutron population of a nuclear system which contains a low number of neutrons. These variables include the moments of the PDF of the neutron population, the maturity time and the source multiplier. The code itself is based upon the mathematics presented in another paper and utilises representations of the neutron population which are independent of both space and angle but allows for the specification of an arbitrary number of energy groups. Five examples of the use of the code are presented. Comparison is performed against results found in the literature and the degree of agreement is discussed. In general the agreement is found to be good and, where it is not, plausible explanations for discrepancies are presented. The final two cases presented examine the effect of the number of neutron groups included and finds that, for the systems simulated, there is no significant difference in the key results of the code. 2017 The Authors. Published by Elsevier Ltd. This is an openaccess article under the CCBY license (http:// creativecommons.org/licenses/by/4.0/).


Introduction
In a nuclear system, each neutron will have a finite probability of causing a fission, being absorbed, escaping the system and so on. It is not possible to know in advance or simulate what fate will befall any given neutron. In addition, the exact time at which a source releases neutrons is not known in advance even if its average intensity is known. In most cases where nuclear systems are being discussed, such as in nuclear power stations under normal operational conditions, there are very large numbers of neutrons present and the change in the total number of neutrons over time may be well approximated by analysing the expected behaviour of this large population of neutrons. This is because the law of large numbers means the overall behaviour in such a case tends towards the mean behaviour.
However, in systems with smaller populations of neutrons, the number of neutrons will not be predictable and repeating identical initial conditions may lead to different results. This is relevant in the case of reactor start-up where it is important to ensure the probability of an undesired stochastic transient is sufficiently low for the start-up process to be classed as safe. There are several methods which have been proposed and, generally, each tends to have strengths and weaknesses.
The CALLISTO Stochastic Point Kinetics code has been constructed based upon the work presented in Williams and Eaton (2017), which provides an in-depth analysis of the physics and mathematics represented here. This code is able to simulate an arbitrary number of prompt neutron energy groups and delayed neutron precursor groups with no spatial or angular dependence. It contains multiple modules designed to calculate the moments and the generating function and its derivatives of the population density function of the number of neutrons in neutron energy groups of interest, the maturity time of the system, the k inf and reactivity of the system and the source multiplier.
This paper summarises the mathematical model contained in CALLISTO code before presenting various example results produced using the code in order to validate or verify the code or to examine the physics simulated by the code.
Mathematical symbols not defined locally in the text are defined in Appendix B.

CALLISTO calculations
Within CALLISTO variables, such as the cross-sections, probabilities of a fission producing different numbers of neutrons and the delayed neutron fractions are defined by the user and may each https://doi.org/10.1016/j.anucene.2017.11.022 0306-4549/Ó 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). be functions of time. For a given time-dependent description of a system, CALLISTO may perform any one of a number of different calculations, such as calculating the k inf or reactivity of the system. This section describes some of the calculations which may be conducted by CALLISTO. This section aims only to summarise the general method as this has been previously discussed in more detail in Williams and Eaton (2017). If more detail on the methods employed or their justification are desired, that paper should be consulted.

k inf and reactivity
Both the k inf and reactivity values of the system in question at a given time are calculated within CALLISTO as a function of time by solving the eigenvalue equation using the power iteration method. The value of k inf rather than k eff is relevant for the systems presented as the systems are considered infinite in extent.

Moments of the probability density function of the number of neutrons in the system
To calculate the mean and standard deviation of the number of neutrons in the system at a given time t the system of equations presented in Appendix A.2 are solved. This solution is performed backwards in time with the variable s representing the time variable which is being advanced backwards. When s ¼ 0 the solution is complete. For instance, the variable N S ðt; UjsÞ represents the mean number of neutrons present in the system in energy range U at a time t due to neutrons released by a source in the time period between time s and time t. As such, by reducing s from t to zero the total number of neutrons present in the system at time t due to neutrons produced by sources since time 0 can be obtained.
This approach means that the calculation of the moments of the PDF at each value of t is independent of the calculation of the moments of the PDF at every other value of t.

Generating functions and derivatives of the probability density function of the number of neutrons in the system
The generating function of the PDF of the number of neutrons in the system in energy range U at a time t due to neutrons released between a time s and a time t is defined by: where z is the generating function variable and Pðn; t; UjsÞ is the probability of there being n neutrons in the system in energy range U at a time t due to neutrons released between a time s and a time t.
To calculate this value and its first and second derivatives with respect to z at a given time t the system of equations presented in Appendix A.1 are solved. This solution is performed backwards in time with the variable s representing the time variable which is being advanced backwards. When s ¼ 0 the solution is complete.
This approach means that the calculation of the generation function and derivatives at each value of t is independent of the calculation of the generating function and derivatives at every other value of t

Maturity time
The maturity time is the time at which the rate of change of the RSD of the number of neutrons in the system becomes small as the RSD approaches an asymptotic value. This is also the time at which the RSD in the number of neutrons in the system first becomes close to the RSD of the number of delayed neutron precursors in the system, as demonstrated in Fig. 1. What is deemed ''small" or ''close" is not uniquely defined meaning the maturity time is not uniquely defined and is to some degree subjective.
The maturity time is found by solving the equations presented in Appendix A.2 with the appropriate final conditions for different values of t until the time t mat is found such that the following equation is satisfied: where mat is the maturity time convergence criterion, which is set to a value of 1 Â 10 À5 for the calculations performed in this paper. RSD neutron ðtÞ is the RSD of the number of neutrons in the system at time t and RSD precursor ðtÞ is the RSD of the number of delayed neutron precursors in the system at time t. The choice of mat effectively selects what is considered ''close" in terms of the difference in the RSDs of the neutrons and delayed neutron precursors. We may examine Fig. 1 in a little more detail. When we consider the state of the system after a very short period of time dt the only event that may have occurred is the release of a single neutron from a source disintegration. This means the PDF of the number of neutrons in the system is given by: Pð1Þ ¼ Sdt; As a result, both the mean and the variance of the number of neutrons in the system is equal to Sdt, meaning the RSD is equal to 1 ffiffiffiffiffi Sdt p . A similar argument may be applied to the precursor population, but here the mean and variance will be proportional to dt 2 as both a source emission and a fission must occur to produce the first precursor. It follows that the precursor RSD will be proportional to 1 dt . In the middle region of this figure it is helpful to recall that, as the system is delayed super-critical (i.e 0$ < q < 1$), any chain of prompt neutrons will die out on a timescale equal to a fairly small number of prompt neutron lifetimes. As a result, the neutron population is made up of short-lived chains of prompt neutrons caused by the release of a neutron from a source or the decay of a delayed neutron precursor (with the latter becoming progressively more important as the number of delayed neutron precursors increases). While one of these chains persists it may cause a large number of Fig. 1. The RSDs of the neutron and delayed neutron precursor populations for a sample system which becomes critical at t $ 120 s and prompt super-critical at t $ 350 s (with the system having no neutrons or precursors present at t ¼ 0 and source which begins releasing neutrons at t ¼ 0). Note that the neutron population refers to prompt neutrons created in a fission, delayed neutrons produced by the decay of delayed neutron precursors or neutrons released in source disintegration.
neutrons to be present in the system (of the order of the tens or hundreds of neutrons at its peak).
Near the start of the middle region, the mean neutron population will be low. This corresponds to there being a small probability of a non-zero number of prompt neutron chains being present in the system. This means that there is a high probability that there will be no neutrons present and a small probability of there being a larger number of neutrons present (of the order of the tens or hundreds of neutrons that a single chain can cause). For a given specified history of delayed neutron precursors the number of chains of neutrons present will be represented by a Poisson distribution with a mean of: where s chain is the duration of a single chain of neutrons. In reality this will vary from chain to chain but we will assume it constant for simplicity here. The assumption of a Poisson distribution is valid if the source intensity, number of delayed neutron precursors and reactivity do not change on timescales smaller than s chain . Thus, the RSD of the number of chains in a given realisation of the system is given by 1 ffiffiffiffiffiffiffiffi n chain p and this is proportional to the RSD of the number of neutrons. This represents the RSD in the number of neutrons in a given realisation but, as different realisations will have different numbers of delayed neutron precursors, this is only one contribution to the RSDs of the ensemble of realisations. The remainder is due to the RSD in the numbers of delayed neutrons precursors across realisations. As time increases the number of delayed neutron precursors will increase, causing the probability of there being zero neutron chains present to drop and the probability of multiple chains being present to increase. This causes the RSD to decrease as the mean number of neutrons chains present increases.
In qualitative terms, the RSD of the delayed neutron precursors may be expected to be lower as the lifetime of delayed neutron precursors is much longer (0.1 s-100 s). This causes a much smaller RSD in the number of delayed neutron precursors as their population is not varying so rapidly or to such extremes as the population of neutrons.
At later times the RSD of both the neutrons and precursors tends towards the same steady value. This is because the component of the RSD of the neutron population relating to uncertainty due to the different number of neutrons chains for a given number of precursors becomes small as n chain increases. This means the only component of the RSD remaining is that of the number of delayed neutron precursors themselves.
To examine this mathematically, we use the fact that, for a highly multiplying medium, the PDF takes the form of a gamma PDF (Radkowsky, 1964, p. where g ¼ hni 2 r 2 . We observe that, as t ! 1; g ! g 1 (a constant). Thus we may write (7) as: Similarly, the PDF of the precursor population in precursor The variance of the PDF is: Setting x ¼ r 2 n hnðtÞi 2 , we find that Eq. (9) becomes: Thus, the RSD rn hni is the same for all self-similar PDFs including the neutron population and the delayed neutron precursor populations. As such, the maturity time may be seen as the time at which these populations may be considered self-similar.

Source multiplier
The source multiplier is the factor by which the source must be multiplied such that the probability of the number of neutrons present in the system at the maturity time (found with the original source) is less than or equal to the mean number of neutrons present at this time in the system with the original number of neutrons is equal to a given probability Q. To illustrate this, consider the probability that the neutron population is less than n Ã : If n Ã is the mean number of neutrons in the system at time t mat with the unmodified source then the source multiplier may be thought of as the factor the source must be multiplied such that the probability calculated in Eq. (11) is equal to some prespecified value Q. To put it mathematically: This value may be used in reactor design in order to account for the effects of a low number of neutrons in a system on safety. To demonstrate this, consider a case where a reactor is being switched on with a linear increase in reactivity in the presence of a source of a particular intensity. In this case a power peak may be formed as the neutron population increases until its increase is prevented by changes in reactivity brought about by the control system or negative feedbacks in the system due to the high neutron population. The timing of the first power peak is non-deterministic due to the low population of neutrons and, in this case, the later the power peak the larger it will be. This effect is discussed at length in Cooling et al. (2016).
In the above case, it may found through deterministic analysis that, during a reactor start-up, a particular combination of reactivity profile and source intensity produces a mean number of neutrons as a function of time which is on the limit of what is considered safe in terms of the power of the system. However, in a given realisation of the system, the actual number of neutrons present is not necessarily equal to the number predicted in the deterministic model. We proceed by considering the number of neutrons present in a specific realisation of the system at the maturity time (at which time the neutron population is expected to be large enough to be considered deterministic), which is before the power peak. If the number of neutrons in this realisation is lower than the mean number of neutrons at that time, the power peak will be later and larger than the deterministic case and, thus, would be unsafe. As such, it is desirable to limit the probability of such a realisation occurring to being below a prescribed value. To achieve this, the source which produced the just-safe power peak in the deterministic analysis can be increased by a factor equal to the calculated source multiplier relating to the desired probability. Increasing the source intensity changes the PDFs such that the peak power will occur at an earlier time, thus it reduces the probability of a dangerous peak power at a later time.
To calculate this value within CALLISTO, the maturity time is first calculated as described in Section 2.4. Next, the mean number of neutrons present at that time is calculated as described in Section 2.2. Next, a number of different values of the generating function variable z are sampled until the one which corresponds to the desired probability Q at the maturity time is found, utilising the equations in Appendix A.1 and A.3. The number of neutrons this corresponds to, n prob , is found and the mean number of neutrons is divided by this number to find the source multiplier.

Hurwitz curves
Hurwitz Curves are named after their appearance in Hurwitz et al. (1963). These curves display the source multiplier as a function of the related probability Q (see Section 2.5).
To calculate the relevant values, the maturity time is first calculated as described in Section 2.4. Next, the mean number of neutrons present at that time is calculated as described in Section 2.2. Next, a number of different values of the generating function variable z are sampled and the related probabilities Q and the number of neutrons this corresponds to n prob at the maturity time are found, utilising the Equations in Appendix A.1 and A.3. For each value of z the relevant source multiplier is found by dividing the mean number of neutrons by the related value of n prob . Finally, the source multiplier is plotted as a function of Q.

Results
Within this section, different calculations are performed by CALLISTO. In some cases these will be compared to results in the literature or results produced by other codes. In other cases the results explore the capabilities of the code or the implications regarding the underlying physics and models. Hurwitz et al. (1963) provide a number of Hurwitz curves for several reactivity ramp insertions. These curves describe the source multiplier as a function of the corresponding value of Q. Two of those curves are represented here. The transients which produced these curves are characterised by the ratios R k 1 (the ratio of the reactivity ramp rate to the decay rate of the single delayed neutron precursor group) and S k 1 (the ratio of the source rate to the decay rate of the single delayed neutron precursor group). The first case selected has a relatively low reactivity ramp rate, whilst the second has a relatively fast reactivity ramp rate. Both have the same source intensity.

Hurwitz curves
The data for these two cases are given in Table 1. Some data is duplicated directly from Hurwitz et al. (1963) whilst other values, which are not given in that paper, are approximations. The results are not expected to be particularly sensitive to the values which have been approximated.
The results of these two cases are displayed in Fig. 2. The corresponding results from Hurwitz et al. (1963) are also shown for comparison. The agreement for Case 1 is very good and Case 2 is fair. One possible explanation for the deviations in Case 2 is that Hurwitz's model assumed that prompt neutrons have a negligible lifetime whilst CALLISTO does not make that assumption. This is more important in Case 2 as the transient is faster and so the reactivity is changing on a timescale that is closer to that of the duration of a burst of prompt neutrons, rendering the behaviour of prompt neutrons more important.

Multi-group neutron populations
This case is based upon Myers (1995) which presents a number of stochastic calculations of an infinite extent of 93% enriched uranium utilising four neutron energy groups. The neutronics parameters of the system are invariant in time. The calculation neglects the effect of delayed neutron precursors. A summary of the data used is presented in Table 2. Some of this data was not present in Myers (1995) in numerical form but, instead, was present in bar charts from which the data was read manually. This may result in some small differences between the data used here and the data used in the calculations performed by Myers.
The value of k inf , as expected for a static infinite slab of high enrichment uranium, is very high and is also independent of time. The exact value is 2.22408987 as calculated by CALLISTO. A corresponding value does not appear to be given by Myers for comparison.
The mean and standard deviation of the neutron population in energy group 1 as a function of time as calculated by CALLISTO are shown in Fig. 3. The source intensity in question is 100n/s and the time simulated is 1 Â 10 À7 so, in the majority of cases, Table 1 The variables used to define the data for the cases presented in Section 3.1 and selected derived quantities. R and S refer to the reactivity ramp rate and source intensity respectively.
Hurwitz curves generated by Case 1 and Case 2 of the systems described in Section 3.1. The data series entitled ''Hurwitz" were found by reading data from graphs within Hurwitz et al. (1963) and so there is some error attributed to these data series from this process. no neutrons will have been released. However, due to the high reactivity, in the instance that a neutron has been released, the population will grow very quickly. This means the probability distribution of the number of neutrons in the system will contain a sharp peak at zero neutrons and a long tail extending to a large number of neutrons. For this reason, the mean number of neutrons is low and the standard deviation is much higher than the mean.
The mean and standard deviation of the number of neutrons present in the system as a function of time due to a single neutron injected into energy group g at t ¼ 0 is presented in Fig. 4. This again shows how a single neutron injected in the system will give rise to a chain of neutrons which rapidly increases in number. Myers (1995) contains many plots of different variables as a function of time and many of these are comparable to data produced here. However, the data presented is in terms of a ''Non-Dimensional Time t l " which is not simply converted into a dimensional time as a value for l does not appear to be given by Myers. Additionally, no instances where a source is present in the infinite medium are discussed. However, some comparisons may still be made between Myers' results and simulations performed by CALLISTO.
One comparison which may be made is that of the ratio between the mean number of neutrons present in energy group 1 at a given time (once the system has attained a state of exponential growth) for neutrons injected into different energy groups at t ¼ 0. Another ratio that may be compared is the standard deviation divided by the mean for the number of neutrons in energy group 1 for a neutron of a particular energy group inserted at t ¼ 0 The results of these comparisons may be found in Table 3. It should be noted that the ratios from Myers are obtained by manually reading data from a logarithmic axis in a document reproduced from microfilm so there is a considerable uncertainty in the estimate here of the value obtained by Myers.
Overall, the results of Table 3 show similar overall trends but the exact agreement is patchy, with some results being very close and others differing significantly. As noted, obtaining results from Myers (1995) was an imprecise endeavour and this plausibly makes up for a considerable amount of this variation, although it is difficult to quantify this effect.
The results of the third calculation provide the generating function of both the total number of neutrons present in the system due to the source and also the number of neutrons present in the system as a function of time following an insertion of a single neutron of energy group g at t ¼ 0. These results are shown in Fig. 5. Before considering these results, it is useful to recall the definition of the generating function relating to the number of neutrons present at time t due to source operation since time s and its derivatives: where P S ðn; t; UjsÞ is the probability of there being n neutrons in energy range U in the system at a time t given a source that has been present since time s. It is also useful to recall the generating function of the number of neutrons present at time t due to a single neutron injected at time s in energy group g: e G g ðz; t; UjsÞ ¼ 1 À X 1 n¼0 z n P g ðn; t; UjsÞ; where P g ðn; t; UjsÞ is the probability of there being n neutrons in energy range U in the system at a time t resulting from a single neutron injected in energy group g at time s. For this case, we have Table 2 The variables used to define the data for the case presented in Section 3.2.  Fig. 3. The mean and standard deviation of the number of neutron present in the system described in Section 3.2 as a function of time. selected a value for the generating function variable z of 0 and so these simplify to: G S ð0; t; UjsÞ ¼ Pð0; t; UjsÞ; ð17Þ G 0 S ð0; t; UjsÞ ¼ Pð1; t; UjsÞ; ð18Þ G 00 S ð0; t; UjsÞ ¼ 2Pð2; t; UjsÞ; ð19Þ e G g ð0; t; UjsÞ ¼ 1 À P g ð0; t; UjsÞ: In Fig. 5a, the value of G S ð0; t; Uj0Þ reduces from 1 to 0.99825 over the 2 Â 10 À7 s over which these values are calculated, implying that, at t ¼ 2 Â 10 À7 s, Pð0Þ ¼ 0:99825. This matches up well with the fact that source intensity is 1 Â 10 4 n/s. One would expect the probability of zero neutrons being present in a system with time-independent neutronics parameters and with a very short prompt neutron lifetime to be given by: where P E is the extinction probability for a neutron injected from the source. Solving this equation for t ¼ 2 Â 10 À7 s gives an estimate of 0.124 for the extinction probability. We will go on to see that this is indeed a good approximation. Fig. 5b also gives information regarding the survival probability.
Recalling that e G g ð0; t; Uj0Þ ¼ 1 À P g ð0; t; Uj0Þ, we see that the asymptotic value of e G g ð0; t; Uj0Þ should give the survival probability for a neutron injected at t ¼ 0 into group g. Similar data is also available in Myers (1995) and this comparison is tabulated in Table 4. Again, issues relating to the difficulty of extracting data from the graphs present in Myers (1995) causes a significant uncertainty in the values in this column.
The agreement between Myers' results and CALLISTO's result in Table 4 is poor. However, note that the extinction probability estimated in the discussion of Fig. 5a of 0.124 (recalling the source injects neutrons only in group 1) is close to that of the CALLISTO estimate in Table 4 of 0.140.
We may form another estimate of the extinction probability which is actually a lower bound on it by considering the crosssections given in Table 2. We may do this by considering the probability that a neutron in a given energy group will result in a fission. This may be achieved with the following equations: Fig. 4. The mean and standard deviation of the number of neutrons present in the system described in Section 3.2 as a function of time due to a single neutron injected into energy group g of the system at t ¼ 0.

Table 3
Ratios of different moments of the number of neutrons in energy group 1 during the exponential growth phase following the injection of a neutron into the system described in Section 3.2. ni the mean number of neutrons present in the system after a neutron is injected into energy group i at t ¼ 0 and ri represents the standard deviation in the number of neutrons present in the system after a neutron is injected into energy group i at t ¼ 0. Note that difficulty in reading the graphs in which the data is presented in Myers (1995) contributes a significant potential error to results in this column.

Ratio
CALLISTO Myers (1995) 5. The generating functions of the system described in Section 3.2 as a function of time due to either the source or due to a single neutron injected into energy group g of the system at t ¼ 0. All results use z ¼ 0 where R tot;g is the total cross-section of energy group g and p f ;g is the probability that a neutron in energy group g eventually causes a fission. Due to the fact that only down-scatter is present we may easily solve this set of equations to obtain: This provides two pieces of information. Firstly, it provides a lower limit on the extinction probability for each group. If the neutron that is present at t ¼ 0 does not cause a fission then it cannot sponsor a persisting chain. This leads to the relation: where P g ðEÞ is the extinction probability of a neutron injected in group g. As can be seen in Table 4 the results produced by CALLISTO obey this inequality whilst the results obtained from Myers (1995) do not for cases where the neutron is injected in groups 1 or 2. Given the high probability that a neutron will cause a fission, the number of neutrons produced per fission it is reasonable to expect that the true extinction probability will be similar to this lower limit. Again, the CALLISTO simulations are consistent with this expected result. The second piece of information that may be gained from the values of p f ;g relates to the ratios between them. Once a neutron has caused a fission it will produce a number of neutrons as determined by the values of p f ;g;m . We define P E;fission;g to be the probability that a fission caused by a neutron originally injected in group g (noting that it may be in a different energy group when it actually causes a fission) does not produce a chain of neutrons that persists. We thus form the relation: P E;g ¼ 1 À p f ;g ð1 À P E;fission;g Þ: Although P E;fission;g is not easily available we may rearrange and take the ratio of this expression for neutrons of different energy groups to obtain the expression: We tabulate these ratios for both the results produced by CAL-LISTO and Myers (1995) in Table 5.
We note here that, as the fission spectrum is identical for each fission regardless of the energy group of the neutron that caused it, the difference between the values of P E;fission;g must arise from the different fission multiplicities for fissions caused by neutrons of each energy which derive from the values of p f ;g;m in Table 2. A neutron injected in group g is most likely to cause a fission in that group as opposed to any other so it follows that the more likely it is that a fission releases a larger number of neutrons the more likely the chain is to survive. Faster fissions release more neutrons (for this system a fission caused by a neutron in group 1 will release an average of approximately 2.98 neutrons and a fission caused by a neutron in group 4 will release an average of approximately 2.41 neutrons). As a result we expect P E;fission;g to increase for values of g which correspond to lower energy groups. This behaviour is observed in the CALLISTO results but not the Myers results in Table 5.
The analysis of this system using CALLISTO has shown that the results are self-consistent and conform to what may be expected from physical reasoning. The comparison with Myers (1995) has produced mixed results but it has been demonstrated that, where disagreement exists, there is significant evidence to favour the results of CALLISTO. The reason for this discrepancy is unclear. The most obvious explanations are that the system studied by Table 4 The extinction probabilities for neutrons of group g injected into the system described in Section 3.2. Note that difficulty in reading the graphs in which the data is presented in Myers (1995) contributes a significant potential error to results in this column.

Energy Group
CALLISTO Myers (1995) 1 0 :140 0:03 2 0 :175 0:08 3 0 :210 0:16 4 0 :288 0:23 Table 5 The ratios of the extinction probabilities for a chain following fission caused by a neutron injected in group g of the system described in Section 3.2. The values marked ''N/A" are values for which it as not able to obtain a positive value from Myers (1995 Table 6 The variables used to define the data for the case described in Section 3.3.
Variable Value Myers has not been precisely replicated here, the results of Myers could not be read with sufficient accuracy or that there is an error in the code used to produce the results found by Myers. The first two options seem particularly plausible given the poor quality of the available copy of that manuscript.

Delayed ramp
This case is a single-energy case designed primarily to demonstrate CALLISTO's ability to calculate the source multiplier. It contains a ramp reactivity insertion such that the reactivity is as given by Eq. (31): Á $ for 100 s 6 t < 500 s 2$ for t P 500 s The full specification of the relevant parameters may be found in Table 6.
A plot of reactivity against time produced is shown in Fig. 6. The maturity time of the system was calculated to be 378.194 s and the source multiplier of the system for Q ¼ 1 Â 10 À8 was calculated to be 2:41636 Â 10 5 . For comparison, a separate code was used to solve the same equations that are solved by CALLISTO. This code was the code used to produce the data presented in Williams and Eaton (2017). For this case, that code calculated a maturity time of 379.2 s and a source multiplier of 2:40 Â 10 5 . This close agreement helps to verify CALLISTO.
The slight differences are likely to be primarily due to small differences in approach. For example, Williams' code takes its input data in a slightly different form (it takes values of v fgn as opposed to p fgm for example) and this may lead to slight differences in the parameters propagated through to the ODEs as the data will have gone through different degrees of rounding and truncation.

Sixteen group system
This case is intended to explore the multi-group capabilities of CALLISTO. The physically modelled system is that of an infinite mixture of 20% enriched uranium of hydrogen when H:U ratio of 400. The cross-sections of this system are modelled using the sixteen group cross-sections from Hansen et al. (1961) whilst the delayed neutron group data is taken from Wilson and England (2002). CALLISTO calculates that this system has a k inf value of 1.02889, which corresponds to a reactivity of approximately 3.85 $. Into this system a source is placed which has a decay rate of 1000 disintegration/s with each disintegration producing exactly one neutron in the highest energy group.
This system may also be collapsed down into a one-group representation using the following equation: where R onegroup is either the absorption or fission macroscopic cross-section for the one-group case and R g is the corresponding cross-section for energy group g in the sixteen group cross-section data-set. / g is the neutron flux in energy group g in the sixteen group representation of the system in the neutron flux eigenvector of the neutron transport equation solved for the eigenvalue relating to the growth mode associated with the k inf of the system. This is calculated by CALLISTO when the k inf of the system is calculated. The equivalent one-group system is thus also used as an input to CALLISTO and CALLISTO calculates the reactivity of this system to be 3.83$, which is close to that calculated for the sixteen group system. CALLISTO may also calculate the maturity time for both representations of the system. For the sixteen group case this is calculated to be 1:189663 Â 10 À2 s and for the one group system this is calculated to be 1:16218 Â 10 À2 s, showing fairly close agreement.
The mean and standard deviation of the neutron population as a function of time may also be compared as in Fig. 7. It can be seen here that the agreement is fairly good, although the growth in both the mean and the standard deviation is larger in the one group case than the sixteen group case, despite the sixteen group case having a slightly higher reactivity, as noted previously. This is likely due to the one group case having a different PDF of the generation time of a neutron as neutrons do not need to be slowed down before they have a significant chance of causing fission in the one group case. This also explains the slightly shorter maturity time in the one group case.
In fact, we see that the mean and standard deviation of the number of neutrons at the respective maturity times are very similar in the sixteen and one group cases. The mean is 8:485 Â 10 5 and the standard deviation is 6:245 Â 10 6 at t = 1.189663Â10 À2 s for the sixteen group case and the mean is 8:442 Â 10 5 and the standard deviation is 6:311 Â 10 6 at t ¼ 1:16218 Â 10 À2 s for the one group case).
Next we may compare the calculated source multiplier for a target probability of Q = 0.1. In the sixteen group case this is 9:39905 Â 10 6 and in the one group case this is 9:53204 Â 10 6 . The fact that these values agree fairly closely helps to validate the multi-group treatment within CALLISTO and to increase Fig. 6. The reactivity as a function of time of the systems described in Section 3.3 Fig. 7. The mean and standard deviation in the number of neutrons present in the system as calculated for the sixteen group and one group representations of the system described in Section 3.4 confidence that correctly collapsing a complex energy group structure to a one group approximation can still produce useful values.
A final comparison that may be made is that of the running time of CALLISTO for these two different representations of the system. These data are presented in Table 7. As can be seen, the one group calculations are significantly faster. This is primarily because the inclusion of even a single fast group dramatically increases the stiffness of the equation set as a whole, significantly slowing the solution of the ODEs.
As a result, a user may wish to utilise a one group representation of a system in order to decrease the solution time, as the loss in accuracy does not appear to be very large in the simplification to a one group system. However, the system chosen here was deliberately chosen to be solvable in a reasonable amount of time, even with all sixteen energy groups represented. Specifically, the specification of the system ensured that the maturity time was very small, which reduces the amount of time that must be simulated. This has, in turn, placed a lower limit on the value of Q which may be used when calculating the source multiplier. It is not practical to repeat this comparison for a system with a longer maturity time as the sixteen group approximation would take a prohibitively long time to solve. This limits the utility of this result, although Section 3.5 attempts to extend the confidence of the conclusions found here over a larger sets of regimes.

Two group system
Section 3.4 represents a realistic representation of a system whose material properties are derived from a specific material composition in an infinite configuration. However, the inclusion of high-energy groups greatly increases the stiffness of the problem and prohibited simulating the system for a long simulated time. This, in turn, prohibited choosing a low value for Q. In this section a different, less realistic representation of a system without any high-energy groups is presented with the aim of comparing a two group and one group representation of the system regarding the value of the source multiplier obtained. The details of the system being simulated are presented in Table 8. This data is also collapsed into a one-group representation (with the resultant group having an energy of 0.025 eV) using Eq. (32) and the results are compared. Both systems have a reactivity of approximately 0.5$.
Using the two group representation the maturity time is calculated to be 59.35 s and the source multiplier is calculated to be 33.40. For the one group case the maturity time was calculated to be 60.37 s and the source multiplier was calculated to be 33.81 with Q ¼ 1 Â 10 À8 . The agreement of these two values provides a further piece of evidence that the collapse of a multigroup system to a single group can provide an acceptably close agreement at a greatly reduced computational cost where the maturity time is larger and the value of Q is small. The running times are recorded in Table 9.

Conclusions
This paper has presented the CALLISTO code which is based upon the mathematics formulated in another paper (Williams and Eaton, 2017). The code has been applied to a series of five low neutron source verification test cases. In three of these systems the results have been verified or validated against other methods or papers. In two cases the agreement was good, whilst the third case showed significant differences. However, these differences could be at least partially explained and it was demonstrated that CALLISTO appeared to be producing plausible and self-consistent results.
In the fourth system it was demonstrated that, for this system, the reduction in the number of energy groups from sixteen to one did not produce large changes in the results produced by CALLISTO in terms of the reactivity, maturity time, neutron population or source multiplier. In the fifth system a reduction in the number of energy groups was reduced from two to one and it was found that this reduction did not significantly change the maturity time or source multiplier. This provides supporting evidence that a single energy group is sufficient to simulate the effects of a low neutron population in these types of calculations. Both the fourth and fifth cases are consistent with the fact that one-group models of fast burst systems such as GODIVA and CALIBAN are found to provide good comparisons with the available experimental data.
The work performed here could be extended by the introduction of spatial and angular dependencies for the neutron population (Williams and Eaton, 2018). This would allow for the direct simulation of more complex systems. Table 7 The time taken to perform different calculations for the two representations of the system described in Section 3.4. These calculations were all performed on the same computer and include overheads such as reading inputs and creating outputs. These overheads should only make a significant contribution to the one group calculations.

Calculation
Sixteen group One group Mean and Standard Deviation (Fig. 7 3.488824/s Table 9 The time taken to perform different calculations for the two representations of the system described in Section 3.5. These calculations were all performed on the same computer and include overheads such as reading inputs and creating outputs. These overheads should only make a significant contribution to the one group calculations.
: ðA:41Þ To calculate the mean and second moment of the number of neutrons N S ðt; UjsÞ and hNðN À 1Þi s ðt; UjsÞ in energy range U caused Table B.10 Summary of variables and their meanings (Roman letters). Note that a ''0" following a variable name refers to the differential of it with respect to z.

Variable
Description Definition Characteristic energy of neutron for group g Problem-specific F fg ðsÞ Probability of a fission neutron being in group g Problem-specific F dig ðsÞ Probability of a neutron released by decay of a precursor in precursor group i being in group g Problem-specific F sig ðsÞ Probability of a neutron produced in a disintegration of the ith source being of energy group g Problem-specific G Number of neutron energy groups Problem-specific G fg ðz; U; tjsÞ Generating function for a neutron beginning at time s P 1 n¼0 z n Pðn; t; Ujg; sÞ e Gg ðz; U; tjsÞ Generating function for a neutron beginning at time s in group g for neutrons in the energy range U at time t in group g for neutrons in the energy range U at time t 1 À Gg ðz; U; tjg; sÞ G S ðz; U; tjsÞ Generating function for the neutrons released by a source between time s and time t for the neutrons in the energy range U at time t P 1 n¼0 z n P S ðn; t; Ujg; sÞ mn Mass of a neutron 1:675 Â 10 À27 kg The total number of sources. Problem-specific p fgm Probability of a fission initiated by a neutron of energy group g producing m neutrons Problem-specific Pðn; t; Ujg; sÞ Probability of a single neutron injected at time s producing n neutrons in energy range U at time t Variable P S ðn; t; Ujg; sÞ Probability of the neutrons released by a source between time s and time t producing n neutrons in energy range U at time t Variable p sim ðsÞ Probability of a source disintegration of the ith source producing m neutrons Problem-specific Q ðn prob ; t; UjsÞ Probability that the neutrons released by sources in a time period between time s and time t gives rises to less than n prob neutrons in energy range U Pn prob n¼0 P S ðn; t; UjsÞ S i ðsÞ The disintegration rate of the ith source Problem-specific tmat The maturity time: the time at which the RSD of the number of neutrons in the system rSðt;UjsÞ

RagðsÞ
Macroscopic absorption cross-section for group g User defined

Rs;g!g0 ðsÞ
Macroscopic scattering cross-section from group g to g 0 User defined

RsgðsÞ
Macroscopic scattering cross-section from group g P G g 0 ¼1 Rs;g!g0 ðtÞ v fgn ðsÞ The nth neutron multiplicity (prompt only) of fission caused by neutrons in energy group g Pm max;f m¼n m!
ðmÀiÞ! p f m v sin ðsÞ The nth neutron multiplicity of the ith source P mmax;si m¼i m! ðmÀiÞ! p sm