A theory of low source start-up based on the Pál-Bell equations

Article history: Received 18 August 2016 Received in revised form 6 December 2016 Accepted 10 December 2016 The safe start-up of a nuclear reactor depends upon the presence of a steady neutron source in the core. This source, however, does not always have to be physically inserted into the reactor because there exist in the core natural neutron sources from spontaneous fission, cosmic rays, photo neutrons, fission products, etc. Nevertheless, so that the source magnitude is well defined, it is generally thought judicious to have a specially constructed source of the (a, n) type present. From an operational point of view, it is vital to assess the strength of the natural sources to see if they will be sufficient in magnitude to ensure safe stochastic startup without the addition of an extraneous source. The most important case for source evaluation is that of a reactor starting up with fresh, unirradiated fuel because then the natural background sources will be at a minimum. It is the purpose of this paper to examine the criteria necessary to ensure that the source strength is high enough to reduce the probability of any undesirable stochastic transient occurring to a specified value, e.g. 10 ; it may also be considered as an update of the classic work of Hurwitz and co-workers (1963). To carry out the calculations, we use the Pál-Bell backward formalism (Pázsit and Pál, 2008) and apply it to the point model in order to make comparisons with the earlier work of Hurwitz. We also extend the study to include space and energy dependence which are found to have a not insignificant influence on the results. The usefulness of the Gamma distribution is explored and its accuracy assessed. Tables and figures are given to illustrate the conclusions. 2016 The Authors. Published by Elsevier Ltd. This is an openaccess article under the CCBY license (http:// creativecommons.org/licenses/by/4.0/).


Introduction
To understand how a reactor behaves as the control rods are raised, or reactivity is increased in some other way, it is necessary to solve the equations of reactor kinetics. A full space-energy-time dependent study should be carried out if possible, but initially a point model approach is generally adequate for guidance. In conventional reactor startup theory it is generally assumed that the source strength and the neutron density are of sufficient magnitude to reduce any statistical fluctuations arising from the underlying random processes to negligible proportions. For a sufficiently low source strength, however, this may not be the case and there is a possibility of large fluctuations. This may be understood better if we recognise the fact that the concept of criticality does not depend on neutrons. A system may be supercritical by a great margin and nothing may happen. However, as soon as a neutron, the carrier of the chain reaction, is introduced then multiplication can proceed very rapidly. It is this possibility that we must study. A practical example is when a control rod is withdrawn in the presence of a low density neutron field. If there are few neutrons, then the count-rate in a detector will be small and so the operator, or automatic control system, might assume that the rod had not been withdrawn far enough. The rod is then withdrawn further, but at this point the density, which was initially low, may quite rapidly (within a prompt neutron lifetime) become relatively large and, as there is also by then a larger amount of reactivity present (due to the continued withdrawal of the control rod), the doubling time may be very short indeed. It is necessary therefore to specify the source strength such that the probability of the neutron level not exceeding some prescribed value which, on deterministic calculations is safe, is acceptably small (Hurwitz et al., 1963a,b;MacMillan and Storm, 1963). Starting a reactor with no independent source, or a source of unkown strength, is known as a 'blind start' and is to be avoided for the reasons explained above (Shaw, 1969).
We would, in practice, define the source strength and associated reactivity insertion rate by deterministic means such that the neutron level does not exceed a specified value during startup over a given period of time. The probability that this source strength, due to inherent fluctuations, will lead to the neutron level not exceeding that specified value is then calculated using the stochastic theory of neutron transport, as opposed to the deterministic method. Thus both deterministic and stochastic methods are employed. In practice, a scram level is set at which value the reactor is tripped and the control rods inserted to shut the system down. However, even after a very rapid negative reactivity insertion, the power may continue to increase with time above the scram level and will then pass through a maximum, before dropping to a very low level. It should also be recalled that there is a delay time between the trip and the actual motion of the control rod due to inertia. It is the value of the maximum that is the safety level, i.e. the transient following the trip must not cause plant damage either due to the peak power or energy deposited. If the shim rate (i.e. reactivity added per unit time) and source strength are considered to be safe using the deterministic method, we will also need to know the associated probability that the source does not lead to an excessively large power transient arising from the stochastic fluctuations; if it does, then the source strength must be increased until the criterion for safety is met. To put it another way, during a period of low neutron intensity, the reactivity may well increase to a large value, approaching prompt critical, so that when the neutron population does increase, a very severe transient will occur in a short time which is beyond the capability of the control system. Such calculations involve the use of the stochastic equations of neutron dynamics and are the main goal of this work.
In summary, the work involves the calculation of the probability distribution function for the neutrons which is carried out by a tried and tested method, using the saddlepoint approach of Hurwitz et al. (1963a,b) to invert the associated generating function. Also the accuracy of the saddlepoint method is assessed by comparison with an exact inversion formula. We include delayed neutrons, both one and six groups, and come to the conclusion that six groups are essential for an accurate probability evaluation. We also examine the influence of energy dependence of the neutrons and the effect of spatial variation, e.g. the position of the source in the core. Results are presented in graphical and tabular form. In contrast to the work of Hurwitz et al. (1963a,b) which uses the forward form of the equation of probability balance, we shall use the backward form known as the Pál-Bell equation (Pál, 1958a;Bell, 1965); also in contrast to Hurwitz et al., we study the neutron density behaviour rather than that of the precursors. Additionally, we have found that in some circumstances the Gamma pdf is a reasonable approximation to the actual pdf and can be used for guidance regarding the influence of a low source without having to calculate the exact pdf.

Start-up behaviour based on deterministic reactor kinetics
Before commencing our stochastic studies it is useful to describe the time dependent behaviour of the reactor by means of the conventional equations of reactor kinetics. These equations regard the neutrons as having a well-defined mean and no fluctuations. We begin with the point reactor kinetics equations in the one speed, I delayed neutron group approximation, as follows (Keepin, 1965 NðtÞ ð 2Þ NðtÞ is the mean neutron number and C i ðtÞ the mean number of delayed neutron precursors of the ith group. The initial conditions are Nð0Þ ¼ C i ð0Þ ¼ 0. Note that N; C and S denote the total number of neutrons, precursors and source strength in the system and not the number per unit volume. The other parameters are defined below, S(t) = independent source strength, neutrons/s k i = average delayed neutron precursor decay constant of ith group, sec À1 b i = delayed neutron fraction of ith group a i b m = mean number of neutrons emitted per fission (prompt + delayed) k a ¼ vR a ¼ vðR c þ R f Þ k c þ k f = absorption probability per unit time per neutron; this is sometimes referred to in Markov processes as the absorption transition rate R is the macroscopic cross section The parameters in the above equations are suitably averaged over energy and homogenised. These kinetic equations may be solved for a variety of reactivity variations, however as the most likely change of reactivity in a typical power reactor is by withdrawal or insertion of control rods, we assume that reactivity changes are made by changing the capture cross section. Scenarios in which the power is allowed to rise to a specified fiducial level and then shut down have been studied. In particular, we are interested in defining a value of startup source strength S m and reactivity rate R m such that the power increases smoothly to the operating level without any chance of an energy release that would cause reactor damage. The combination (S m , R m ) will later be used in the stochastic formalism to estimate the probability that a 'rogue transient' will occur and also to calculate, by that same stochastic method, by how much one should increase the source strength, or reduce the reactivity rate, to ensure an acceptable degree of safety. The deterministic calculations are not reported here as they are well-known, but we mention them to indicate how the stochastic procedures enter the problem.

Stochastic behaviour
As we noted above, when the source strength is low there will be large fluctuations in the neutron density. To assess the influence of these fluctuations on the start-up behaviour we must write down the generating function equations for the probabilities. Thus, let us define Pðn; c 1 ; c 2 ; . . . c I ; tÞ as the probability distribution function (pdf) that at time t there are n neutrons present and c 1 precursors of type 1, c 2 precursors of type 2. . .c I precursors of type I. The generating function is defined as Fðx; y 1 ; y 2 ; . . . y I ; tÞ ¼ x n y c 1 1 y c 2 2 . . . y c I I Pðn; c 1 ; c 2 ; . . . c I ; tÞ The equation obeyed by F is (Pázsit and Pál, 2008;Williams, 1974;Courant and Wallace, 1947) @Fðx; y; tÞ @t ¼ S d ðtÞ X K k¼1 q k x k À 1 ! Fðx; y; tÞ þ X I i¼1 k i ðx À y i Þ @Fðx; y; tÞ @y i þ k c ðtÞð1 À xÞ þ k f ðtÞðf T ðx; yÞ À xÞ Â Ã @Fðx; y; tÞ @x ð4Þ where y ¼ ðy 1 ; y 2 ; . . . y I Þ. Eq. (4) is known as the forward form of the generating function equation. The initial conditions normally used are that, at t = 0, there are no neutrons or precursors present. Thus The subsidiary generating function f T ðx; yÞ is defined as . . . X 1 c I ¼0 x n y c 1 1 y c 2 2 . . . y c I I pðn; c 1 ; c 2 ; . . . c I ; tÞ ð 6Þ where pðn; c 1 ; c 2 ; . . . c I Þ is the probability that n prompt neutrons are emitted in a fission event, together with c 1 precursors of type 1, c 2 precursors of type 2. . . c I precursors of type I. For our purposes we write (Pázsit and Pál, 2008;Williams and Pazsit, 2015) pðn; c 1 ; c 2 ; . . . c I Þ ¼ p n p c 1 p c 2 . . . p c I i.e. the emission of prompt neutrons and precursors are independent events. This may not be strictly true, but there is no experimental evidence to suggest that this is not a reasonable assumption. The generating function then becomes  (Terrell, 1957, Diven et al., 1956. A more general way of writing the generating function is due to Bell (1965). Thus we write f(x) as where N is the maximum number of prompt neutrons emitted in a fission event and m refers to prompt neutrons. Now let us use the binomial theorem to get ð1 À ð1 À xÞÞ m ¼ X m n¼0 ðÀ1Þ n m n ð1 À xÞ n  This formulation is consistent with the definition of f ðx; yÞ given by Pázsit and Pál (2008) in the sense that they give f x ð1; 1Þ ¼ mð1 À bÞ; f xx ð1; 1Þ ¼ ð1 À bÞ 2 m 2 À ð1 À bÞ m; f xy i ð1; 1Þ ¼ bð1 À bÞ m 2 Values of v n for prompt neutrons are given by Holden and Zucker (1986) for various fissile elements. Bell also introduces In Eq. (4) there is a source term is the number of disintegrations per second and q k is probability that the source emits k neutrons per disintegration. This term allows for the possibility that more than one neutron is emitted in a disintegration, as for example in correlated sources from spontaneous fission or spallation sources or a sum of such sources. One can therefore define a new subsidiary generating function as The average number of neutrons emitted per second is then 1 and if only one neutron is emitted per disintegration then S d is the number of neutrons emitted per second. The general term f q ðxÞ may be described as a compound Poisson process. Finally, we give in Table 1 below the values of v n from Holden and Zucker (1988). It should be noted that in principle it is possible to extend the forward equation to energy and space dependent processes (Stacey, 2001, Matthes, 1962, 1966 but a more convenient method will be described in the next section. To obtain the pdf it is necessary to solve the equation for the generating function F and then invert it to obtain the actual pdf Pðn; c 1 ; c 2 ; . . . c I ; tÞ. Eq. (4) is a first order partial differential one and it requires some numerical effort both to solve it and to invert the resulting solution. Moreover, as we will need to consider energy and spatial variations, the forward form of the stochastic balance equation is not ideal for extension to this situation as has been shown by several authors and explained in some detail in Williams (2008). For this reason a complementary approach is adopted based on the so-called backward equation of probability balance and encapsulated in the Pál-Bell equations. In their simplest, one speed, point model form these may be written as (Pázsit and Pál, 2008) Table 1 Multiplicities for fission at 2200 m/s (Holden and Zucker, 1988). where Gðz; tjsÞ is defined as the generating function for Pðn; tjsÞ which is the probability that if a single neutron is born at time s it will give rise to n neutrons at time t later. G di ðz; tjsÞ is an analogous quantity for the delayed neutrons and G S ðz; tjsÞ defines the pdf that results when a source of strength S d disintegrations per second is present. These three coupled non-linear equations may seem more complicated than the single, first order partial one for the forward equation, but in practice that is not so because the above are first order, ordinary differential equations which may be readily solved by standard numerical methods. It is clear also that Eqs. (7)- (9) are the characteristic equations of the partial differential equation (4). We use Eqs. (7)-(9) throughout, together with modifications to include space and energy dependence. The final conditions on G; G di and G S , i.e. when s ¼ t, are given by Gðz; tjtÞ ¼ z, which corresponds to one neutron present, G di ðz; tjtÞ ¼ 1, which corresponds to no precursors present. G S ðz; tjtÞ ¼ 1 which assumes that the number of neutrons in the subcritical state before startup is negligible; Hurwitz et al. also made this assumption. Should the latter assumption be invalid, then it would be necessary to calculate the associated steady state form of G S . Details of how this is done are given in Appendix D, but experience shows it to be unecessary.

Harris' approach
In the early 1960's when stochastic startup problems were first being seriously studied, Harris (1964) and Hurwitz et al. (1963a,b) were both involved in the development of the theory. Although the general concepts of Harris and Hurwitz were similar, their numerical approaches differed markedly. Indeed, as far as we can tell, Harris has not given in the open literature any explicit values for the modified source strength required to suppress any rogue transients but has simply noted that a safety probability (SP) can be defined by where t h is a time at which the reactivity has achieved a large value and the neutron and precursor levels are NOT low, i.e. the sum term is very small. He goes further to simplify this criterion to Pðn; c; t h Þ ð 11aÞ As the neutron and precursor probabilities are strongly correlated, Harris argues that if the system is super-prompt critical one may write and if it is sub-prompt critical Pðn; c; t h Þ ð 11cÞ A number of semiqualitative figures are constructed of neutron density vs precursor density and regions of 'acceptability' defined. Harris also argues that the Gamma probability distribution is an acceptable form for the true distribution. No numerical values are given and while one follows the argument, this does not seem a very profitable route to follow at the present time. However, in order to check Harris' proposal of using the Gamma pdf, we will assess its accuracy later in the paper.

The method of Hurwitz et al.
An important characteristic of a safe start-up is that the power level must rise to a detectable level before the reactivity has risen to an undesirably large value (Hurwitz et al., 1963a,b). A primary objective of the statistical calculation is therefore to demonstrate that the probability that the power remains low, while the reactivity is appreciably above delayed critical, is vanishingly small for the classes of start-ups which may occur in practice. In other words we want the power to be relatively high as soon as possible. A calculational procedure for attaining this objective is to choose a fission rate F at which statistical fluctuations are definitely unimportant (e.g. at the maturity time t mat , which is the time at which the system becomes dominantly deterministic) and then to determine the earliest time t at which, with probability 1 À e, the actual fission rate exceeds this chosen rate. e is a very small number, e.g.
$ 10 À8 , which reflects the maximum uncertainty that one is willing to accept in describing start-up. We see the beginnings of some degree of subjectivity here, a characteristic common to all risk assessments. The subsequent behaviour of the system for t > t mat is calculated by the deterministic equations of reactor kinetics.
The reasons for examining the nature of the weak source effect is because, if the neutron density is low and the system is supercritical, there is the possibility of single persistent chains developing from one source neutron. Such chains can lead to large energy bursts at random times which are not predictable by the normal equations of reactor kinetics (Hansen, 1960). Only by ensuring that the number of neutrons present is large at the prescribed time (or not too small) can single persistent chains be avoided. That is done by having a sufficiently strong source which leads to the generation of many overlapping chains which are, by definition, predictable. This situation can easily arise in reactor startup as control rods are being withdrawn and so our goal is to obtain a relationship between the source strength and the probability that the power remains low in the presence of a high reactivity. In principle if we can obtain the pdf, Pðn; tÞ, of the neutron population in the reactor, then we can define Pðn; tÞ ð 12Þ as the probability that the neutron population is less than n ⁄ . Now it turns out that a characteristic of the pdf is that, as time proceeds, Pðn; tÞ ! P n nðtÞ as t ! t mat ð13Þ i.e. at some time t mat , the pdf becomes essentially a function of n= nðtÞ only and hence The time t mat at which this occurs is called the maturity time; it is not a precise time but can be defined when the relative standard deviation, rðtÞ= nðtÞ, reaches a constant value within a certain tolerance, e.g. 10 À3 . The operational requirement is that Q ðn Ã ; t mat Þ be a specified value (e.g. 10 À8 ). To put it otherwise, t mat is the time before which (to a given probability) no undesirable stochastic transients will occur. For t > t mat we are of course in the deterministic region and the normal rules apply. The value of the power at t mat is generally low and of the order of 10 W and is therefore a safe condition with no significant feedback. We wish to know what source strength will ensure that Q ðn Ã ðt mat Þ= nðt mat ÞÞ is equal to the specified value, say Q Ã . The neutron population n Ã ðt mat Þ is now regarded as the value of n associated with a just safe source as obtained from the deterministic calculation that provided the just safe pair ðR m ; S m Þ. Once t mat has been reached, we may use the deterministic equations of reactor kinetics to describe the neutron behaviour and n Ã ðt mat Þ will be the initial condition. Note that to achieve Q Ã we may alter the source strength or the rate of reactivity insertion but in practice it is usually more convenient to change the source strength S. We now assume that nðt mat Þ is associated with a new source S which is to be adjusted so that Q Ã is achieved. Thus we can write n Ã ðt mat Þ / S m and nðt mat Þ / S ð15Þ This argument can be made clearer if we note that the desired value of n Ã ðt mat Þ is determined by the value chosen for Q from the relation Pðn; t mat Þ Our calculations will give us the value of n Ã ðt mat Þ for the just safe source S m corresponding to Q and the associated mean density nðt mat Þ. What we now need to calculate is the value of the new source strength S which will make n Ã ðt mat Þ ¼ nðt mat Þ. The ratio nðt mat Þ=n Ã ðt mat Þ is read from what we term 'the Hurwitz curve' and will be calculated below. It is also a useful exercise to use the corrected source S as the initial 'just safe' source. If this is done then the resulting factor nðt mat Þ=n Ã ðt mat Þ is close to unity, which is a comforting result. The fact that the revised ratio is not unity we believe to be due to the fact that, even in the so-called deterministic regime, there remain fluctuations albeit of a very small maginitude.
In physical terms, Q Ã ¼ Pðn < n Ã ðt mat ÞÞ is made small to ensure that the neutron density is not low below n ⁄ and hence avoid the first persistent chain surges about which we wrote above. We believe that this argument for specifying how to calculate the new source strength is rather clearer than that given by Harris as described above. It is also quantitative and we will describe below how the maturity time and its associated parameters are calculated. Hurwitz et al. (1963a,b) use the saddlepoint method (Wilf, 1990 and Appendix F) to find an approximate inversion of the generating function G S ðz; tjsÞ which leads to Q . Our work is based upon that technique which we show is very accurate, and very fast to evaluate numerically. Finally, we note that another explanation of the meaning of Q Ã may be deduced if we write it in the continuous form Q ðn Ã ; t mat Þ % We now interpret this as the probability that there will be an adverse event triggered by a first persistent chain between times t ¼ 0 and t ¼ t mat . Or, alternately, because the region ð0; t mat Þ is the stochastic regime, then 1-Q is the probability that the system can traverse the stochastic regime without producing a stochastic surge due to a persistent chain. Of course, this says nothing about the magnitude of the surge should it happen. For that we must turn to the work of Hansen (1960). Hurwitz et al. (1963a,b) (designated HMSS) derived an equation analogous to (4) to solve the low source startup problem. We will briefly explain their method and why it is now outdated. HMSS published two definitive papers on the subject, in one of which they used the zero prompt neutron lifetime approximation with one group of delayed neutrons. To describe their approach we find it clearer if we use their equation which is given, in the notation used by HMSS, as @Fðu; v; tÞ @t ¼ Sðu À 1ÞFðu; v; tÞ þ kðu À vÞ @Fðu; v; tÞ @v þ 1 Casting this in the notation of Eq. (4), we find @Fðu; v; tÞ @t ¼ Sðu À 1ÞFðu; v; tÞ þ kðu À vÞ @Fðu; v; tÞ @v þ ½k c ð1 À uÞ þ k f ðgðuÞð1 À mbð1 À bÞð1 À vÞÞ À uÞ @Fðu; v; tÞ @u where the generating function for prompt neutron emission gðuÞ is given by  Table 1 and only go up to N = 5. The form of Hurwitz's equation is identical to our Eq. (4) if the term bð1 À bÞ $ b in the square bracket, which introduces a negligible difference of Oðb 2 Þ. We will use the HMSS notation in the above equation to explain their approach. The first assumption is that the reactivity does not exceed prompt critical which allows (approximately) the neglect of the prompt neutron lifetime s. When the reactivity is close to delayed critical, the insertion of a single source or delayed neutron will cause a chain of prompt fissions which will persist for the order of magnitude of about a hundred generations (a generation lasts about s seconds).
The prompt chain during this time is of order 10 À3 s. During this time, the reactivity changes very little and the population of delayed neutron precursors will change almost exclusively by virtue of the fact that additional precursors may be produced in the course of the prompt fission chain. It is therefore legitimate to make the assumption that the prompt fission chain occurs in an extremely short time. This is equivalent to the assumption that the prompt neutron lifetime is zero. To proceed, HMSS argue that it is only necessary to know the statistical distribution of the number of precursors produced in a single prompt fission chain which they show is equivalent to Pð0; m; 1Þ, i.e. the extinction probability of neutrons given m precursors. This is obtained from the above equation for F where k p is time independent and the source and the delayed neutron terms are ignored (over this short time span), viz: @Fðu; v; tÞ and hence by inversion, still using the saddlepoint method, to obtain P d ðm; tÞ ¼ P n Pðn; m; tÞ. It must be clear by now that the methods used by HMSS are both approximate and very tedious numerically. It is also restricted to one group of delayed neutrons which, it will be seen, is a severe limitation. It turns out that although these methods are ingenious, and in some cases accurate, use of the backward formalism removes the need to make any approximations other than that of the saddlepoint method, which in fact is shown to be very accurate. Although the HMSS method is not used here, we have cited it for completeness and historical purposes. As far as can be seen from the literature, HMSS were the first to deal in any depth with this matter.

The low source problem via the Gamma distribution
It is well known (Bell, 1963;Longmire et al., 1950;Harris, 1964) that a very reasonable approximation to the actual pdf is the Gamma distribution and below we go through the necessary steps using the Gamma pdf and calculate Q for various cases. Later, we will develop a method for solving the exact backward generating function equations which will enable the accuracy of both the Gamma pdf and other approximations to be assessed. Harris (1964) gives a convincing proof that, in a highly multiplying system, the pdf will always tend to that of a Gamma distribution. The Gamma distribution has two adjustable parameters and may be written as where the mean and variance are l ¼ ab and r 2 ¼ ab 2 . If now we calculate the exact mean and variance of the pdf we find l ¼ nðtÞ and r 2 ðtÞ ¼ n 2 ðtÞ=gðtÞ. Thus we may write P(x) as P S ðn; tÞ ¼ gðtÞ nðtÞCðgðtÞÞ where gðtÞ ¼ n 2 ðtÞ=r 2 n ðtÞ, which will be shown to tend to a constant value as t ! 1 or more practically as t ! t mat ; note that a misprint occurs in Harris's work and he has written gðtÞ ¼ nðtÞ=r n ðtÞ rather than the correct value above. This does not seem to have affected any of Harris's numerical work and so we assume that it is a genuine typographical error. Using the expression for gðtÞ enables us to write the limiting form of the pdf in terms of n= nðtÞ, viz: P S ðn= nðtÞÞdn ! g mat dn nðtÞCðg mat Þ g mat n nðtÞ g mat À1 Â exp À g mat n nðtÞ ; for large t where g mat ¼ gðt mat Þ. In practice the statement 't ! 1' really means as 't ! t mat ' as will become clear when we present numerical results below. The pdf is said to mature when it is essentially a function of n= nðtÞ only. This is a vital property of all pdf's describing multiplying processes and even the exact pdfs have this property. In fact the property arises from a knowledge of the first and second moments of the pdf, i.e.
nðtÞ and r 2 ðtÞ , equations for which can be easily derived (Harris, 1964;Pázsit and Pál, 2008;Williams 2016). We also note that the relative standard deviation ffiffiffiffiffiffiffiffi ffi gðtÞ p Þ is a good measure of the magnitude of the fluctuations and will be used in the criterion to define t mat . The corresponding cumulative distribution of the Gamma pdf is by definition Cðgðt mat ÞÞ Before exploring the specific details of the low source problem, let us examine the actual random nature of the problem by simulating the neutron density by sampling from P S ðn; tÞ as given by Eq. (19). We write the pdf as There is an IMSL library subroutine for this simulation, namely DRNGAM and DSCAL. For data we use S=k ¼ 5000 and R=k ¼ 0:2 and initial reactivity of qð0Þ ¼ 0, where R is the reactivity insertion rate in $/s, k is the mean value of the delayed neutron decay constants (=0.08519 s À1 ). We then find the result in Fig. 1 for one group of delayed neutrons. The full line denotes the mean density nðtÞ < nðtÞ >, whilst the fluctuating line denotes a simulation from the Gamma pdf at various times. It is clear that as time proceeds the fluctuations damp out and the curve tends closely to the mean value. However, even when the so-called deterministic (or mature) region is reached there still remain some relatively small fluctuations. These can be seen in the figure by noting that gðtÞ ! constant after about 55 s; a somewhat more physical parameter is the relative standard deviation r= n which is also shown in the figure. This reduces from around 7, at the beginning of the transient, to the asymptotic value of 0.36 in the deterministic region. It is interesting to note that the fluctuations in the early part of the transient have magnitudes that are well below the mean value. This can be understood from the associated pdf where g < 1 and hence biases small n values. It can be seen from Fig. 1 that as soon as g > 1, the fluctuations start to reduce markedly in magnitude and the transition to the deterministic regime is quite abrupt. The residual fluctuations for t > t mat , which are always present, are known as reactor noise and are useful diagnostic indicators regarding the physical condition of a power reactor (Williams, 1974).
In order to assess the accuracy of the Gamma pdf, Fig. 2 shows the Gamma pdf compared with the exact value as computed using the Abate-Whitt inversion formula at different times (60-85 s) as shown on the individual curves (Abate and Whitt, 1992, see Appendix A). It will be noted that over the most prominent part of the curves the Gamma pdf follows the exact form very closely. However in the 'wings' of the distribution, for very large and very small values, the Gamma pdf deviates markedly from the exact value although this is not evident from the figure. This feature will show up in the curves of the cumulative pdf to be discussed later. The error in the wings of the Gamma pdf is unfortunate because, in the low source problem, it is this region of small Q(n ⁄ , t) which is the most important. Notwithstanding these limitations, as we will see, there are some situations where the Gamma as 'pdf is very useful and in particular for scoping calculations to assess the influence of a low source (Williams, 2016).
The number of delayed neutron groups used in the calculation can prove crucial and Fig. 3 shows the Gamma pdf and the cumulative pdf for one and six groups of delayed neutrons; it is observed that there are significant differences. These differences are discussed in more detail below.

Calculation of cumulative pdf from the generating function equations
In order to calculate Q without appealing to the approximate Gamma pdf, it is necessary to solve numerically the Pál-Bell equations, as discussed above, for the generating function and then to invert that generating function. There are several techniques for inverting a generating function to recover the pdf and we will discuss some of the more useful of these.

Saddlepoint method
The basic problem is to evaluate the quantity where the generating function for P S ðn; tjsÞ is G S ðz; tjsÞ ¼ X 1 n¼0 z n P S ðn; tjsÞ Now it is easily shown from the definition of the generating function that and a prime on a symbol means differentiation with respect to z. We suppress the arguments of G S for notational simplicity. For the case when an analytical expression is not available, which applies to most situations, we have to deduce the saddlepoint quantities directly from the differential equations. Thus we need equations for G S ðzÞ; G 0 S ðzÞ and G 00 S ðzÞ where primes again denote differentiation with respect to z. This is explained further in Section 9 with numerical results.

The Abate-Whitt method
The saddlepoint method is approximate and to assess its accuracy it is necessary to compare it with an exact method. For the probability that there are no neutrons present, i.e. Pð0; tÞ we need only set z = 0 in the Pál-Bell equation and solve for Gð0; tjsÞ. For the first few values of Pðn; tÞ; n ¼ 1; 2; 3 . . . we may use the relation and apply the operation directly to the Pál-Bell equation. However for large values of n, this is impractical and we must use the method described by Abate and Whitt (1992) which is essentially exact and is therefore, in principle, better than the saddlepoint method. Details of how to use this method in a practical manner are given in Section 10, and Appendix A describes the mathematics behind the method. The major disadvantage of the Abate-Whitt method is the time of execution. This can be hours compared with less than a minute for the saddlepoint method. We have assessed the accuracy of the saddlepoint method for the case of one group of delayed neutrons by calculating the value of Q ðn; t mat Þ vs n= nðt mat Þ and comparing it with the exact method; numerical results will be given below. We are now in a position to discuss the general form of the Pál-Bell equation in space and energy and in a form applicable for use in practical low source studies.

Backward equation
The alternative form of equation for calculating probability balance, known as the backward equation, has been discussed for a point model above but it is far more general and will describe both energy and space dependent problems. In discussing the backward equation, we shall show it for phase space, i.e. position and velocity space. Later we show how energy and space dependence can be handled numerically and quantify their effect on parameters of interest. Following Pál (1958a,b) and Pázsit and Pál (2008) it is convenient to define the probability distribution function as which is the probability that at time t the system contains n neutrons in the volume of phase space R, on condition that at an earlier time s a single neutron was injected into the system at the point r 0 with velocity v 0 . As the probability in the low source problem is concerned with the total number of neutrons in the entire system, we will regard R as covering the whole reactor volume and all neutron velocities although, in some special cases, it could cover the volume of a detector or source, or an energy range embracing a resonance. The associated generating function is defined as It is shown in Pázsit and Pál (2008) that G satisfies the following generalised transport equation k s ðr 0 ; v 0 ; sÞ is the scattering probability per unit time per neutron (or transition rate) and (33) are defined below. Note that we use a slightly different notation from Pázsit and Pál (PP) in that they use Q to denote the transition rate whereas we use k. Also PP use Q a for their capture (non-fission) transition rate while we use the more traditional k c , with k a ¼ k c þ k f . Eq. (35) may also be written as a differential equation, viz: If the transition rates k x are independent of time, we have Gðz; t; Rjr 0 ; v 0 ; sÞ ¼ Gðz; t À s; Rjr 0 ; v 0 Þ and Eq. (36) simplifies to @G di ðz; t À s; Rjr 0 Þ @s In general k x is not independent of time and so we continue to use the more general formalism. F 0 ðvÞ is the fission spectrum of the prompt neutrons and F i ðvÞ that of the ith delayed group. The functions G di are of course closely associated with the delayed neutron precursors. The final conditions associated with Eq. (33) are In some texts (Pázsit and Pál, 2008), when k x is independent of s, the condition s = t is referred to as the initial condition. In this work we shall use 'final condition' to be consistent with Bell and because in many practical startup problems k x ¼ k x ðsÞ. The generating function Gðz; t; Rjr 0 ; v 0 ; sÞ, as defined above, is the 'single particle generating function' and is always time-dependent as it relates to the chain initiated by a single neutron. V r is a sub-region within the reactor and U v is an energy range; generally the complete range. The usefulness of G is that it does not change when different sources are used and, in a loose sense, is analogous to a Green's function. To relate G to the case where there is an independent source, we proceed as follows. Such a source can itself emit varying numbers of neutrons at each disintegration and this will also be dealt with in the same way as for the forward equation. Now if the source emits neutrons with a compound Poisson distribution with a varying multiplicity, we may write for the source generating function (Bartlett, 1955) with f q ðzÞ the generating function for the source emission as defined above. Note that while G is always time-dependent, G S will be asymptotically time-independent for a sub-critical system. The equation for G S can also be written in differential form as A more convenient notation for representing the above generating function equations has been devised by Bell (1965) in terms of 1 À G, but we defer discussion of that until later. For comparison with experiment it is G S that is used, G and G di being intermediate quantities. The final condition on G S is G S ðz; t; RjtÞ ¼ 1, i.e. there are no neutrons in region R when s = t. In a steady state subcritical reactor we may replace the value of s in Eq. (41) by À1. One of the major advantages of the backward formalism is highlighted by Eq. (42), namely that the bulk of the calculations can be carried out independently of the source. Thus if we had a complex source arrangement such as a spontaneous fission source, uniformly distributed over the core, and a number of different ða; nÞ sources at various positions, the only change in computational procedure would be in the integral on the right hand side of the equation for G S ðz; t; RjsÞ.
A number of functions occurring in the above equations have not been defined; here we do so. Firstly, there is f ðx; v 0 Þ and f i ðy; v 0 Þ which, respectively, refer to the number of neutrons emitted in prompt fission and by delayed neutron precur- where v n are the multiplicities defined by Bell (1965) modified for delayed neutrons, e.g. v 1 ¼ mð1 À bÞ. The function f q ðxÞ in Eq. (41) is analogous to f ðxÞ but refers to source neutrons which may be emission products of other elements.

Point model one speed equations
In the point model, as already noted, we shall not be concerned with the spatial and velocity dependence and so Eq. (33) reduces to Gðz; t; Rjr 0 ; v 0 ; sÞ ! Gðz; tjsÞ, with R covering the whole system and À @Gðz; tjsÞ The generating functions for the delayed neutron probability equations are which we may also write as To obtain the generating function when an independent source is present we use the equation for G S , namely À @G S ðz; tjsÞ @s ¼ S d ðsÞ½f q ðGðz; tjsÞÞ À 1G S ðz; tjsÞ ð 49Þ The initial value G S ðz; tjtÞ, i.e. when s = t, is equal to unity if we assume that there are no neutrons present initially. A formal reduction of the space and energy dependent equation to a point, one speed, diffusion theory model has been described by Bell (1965) and by Pázsit and Pál (2008). The effective point model data will be homogenised over space and energy in the usual manner. Bell (1965) has shown that it is often more convenient to write the backward equation in terms of the functionGðz; t; Rjr 0 ; v 0 ; sÞ ¼ 1 À Gðz; t; Rjr 0 ; v 0 ; sÞ. Thus Eqs. (33) et seq become @Gðz; t; Rjr 0 ; v 0 ; sÞ @s þTGðz; t; Rjr 0 ; v 0 ; sÞ

Bell's notation
where f ½G p and f i ½G di are as defined in Eq. (43). Also we havẽ Thus for the point model we may write the above, using Eq. (43), as @Gðz; tjsÞ @s These equations are consistent with those of Hurwitz et al. (1963a,b) and are related to the equations of the characteristics in the forward form of the probability equation. The source equation becomes, from (49), which will be useful in some limiting cases. Setting N = 2 leads to the quadratic approximation which can also prove very useful.

The diffusion approximation
To include the effect of spatial variation, it will be convenient in many practical problems to use diffusion theory in an homogenised core. The most logical way to arrive at the diffusion theory version of the Pál-Bell equations is to expand Gðz; t; Rjr 0 ; v 0 ; sÞ in spherical harmonics in the angular variable associated with v 0 and apply the usual approximations to the transport equation given by (33). However, Pázsit and Pál (2008) have pointed out that if this procedure is followed, the possibility of setting up a probability balance equation is lost, i.e. the expansion process is not consistent with a true probability balance. Pázsit and Pál argue that the main task is to replace, with an alternative process, the particle streaming operator, v 0 :r r 0 , which represents trajectories that consist of straight sections of random lengths, and whose last section terminates at the point where the scattered neutron is absorbed, scattered or multiplied. Instead, therefore, we need to determine a diffusion kernel which describes the probability density that a neutron existing at a particular point at a given time will appear, through diffusive motion, at another point after a given time. With this diffusion kernel, one can then use the conventional accounting procedure for setting up a probability balance. The details of this process are described in some detail in Pázsit and Pál (2008) and Pál (1962). Although this reasoning is physically appealing, the same outcome arises if a spherical harmonics expansion is used and such an approach also allows an extension to higher orders. The outcome of the above arguments is that, for diffusion theory, we may replace the transport operator, where Dðr 0 ; v 0 ; sÞ is the diffusion coefficient corrected for anisotropic scattering, and the multiplicity term is In the terms above, 0 ; sÞ and we have made the very reasonable assumption that the fission spectra are isotropic. Thus we may now write the energy-dependent diffusion theory version of the Pál-Bell equation as @Gðz; t; Rjr 0 ; v 0 ; sÞ @s þT D Gðz; t; Rjr 0 ; v 0 ; sÞ Explicitly, the multiplicity terms are which does not depend explicitly on velocity other than through the intial value v 0 . When a source is present we have from Eq. (42) À @G S ðz; t; RjsÞ @s ¼ The final conditions associated with the diffusion approximation are G di ðz; t; Rjr 0 ; tÞ ¼ 1 and G S ðz; t; RjtÞ ¼ 1. The condition at a vacuum boundary is Gðz; t; Rjr 0 ; v 0 ; sÞ ¼ 0 and, at an interface between two different media, we have continuity of Gðz; t; Rjr 0 ; v 0 ; sÞ and Dðr 0 ; v 0 ; sÞr r 0 Gðz; t; Rjr 0 ; v 0 ; sÞ.
In terms of the Bell notation G ¼ 1 ÀG, Eq. (62) becomes À @Gðz; t; Rjr 0 ; v 0 ; sÞ The one speed, diffusion approximation with no delayed neutrons is therefore À @Gðz; t; Rjr 0 ; sÞ @s À vr 0 :Dðr 0 ; sÞr 0G ðz; t; Rjr 0 ; sÞ When a source is present which emits one neutron per disintegration, we have @ @s G S ðz; t; RjsÞ ¼ G S ðz; t; RjsÞ Z R dr 0 Sðr 0 ; sÞGðz; t; Rjr 0 ; sÞ ð 66Þ where G S ðz; t; RjtÞ ¼ 1 and we have replaced the volume V r by R. In the buckling approximation, a reasonable estimate for leakage from a bare body would be to write Eq. (65) with k a ! k a þ DvB 2 leading to a quasi-point model. Note, however, that the use of the buckling approximation in the Pál-Bell equation is not correct because the latter is non-linear, although this is unlikely to cause any significant error for small buckling values. However, a more satisfactory procedure would be to solve Eq. (65) directly by numerical methods; this will then allow for the spatial variation of the source position. If delayed neutrons are included in the diffusion theory approximation, the one speed generating function equation becomes À @Gðz; t; Rjr 0 ; sÞ @s À vr 0 :Dðr 0 ; sÞr 0G ðz; t; Rjr 0 ; sÞ ¼ k f ðr 0 ; sÞ À k a ðr 0 ; sÞGðz; t; Rjr 0 ; sÞ À k f ðr 0 ; sÞ plus the equations forG di , viz @G di ðz; t; Rjr 0 ; sÞ @s ¼ k iGdi ðz; t; Rjr 0 ; sÞ À k iG ðz; t; Rjr 0 ; sÞ ð 68Þ The equation for G S ðz; t; RjsÞ remains as in Eq. (66). More will be said about these equations below.

Mean value and variance equations
In order to calculate the safety probability, it is necessary to know the border line between stochastic and deterministic behaviour. We have discussed above that this occurs at the maturity time which is when the pdf becomes a unique function of n= nðtÞ only. Such behaviour occurs when the relative standard deviation, , becomes essentially independent of time; thus we need the mean and variance. To calculate these quantities consider Eqs. (44) and (45) and from them derive the mean value equations for the neutron density and the precursor concentration and also the variance. Let us write Eq. (44) (suppressing s and t) as Let us write Eq. (45) as Differentiating this with respect to z, we find We now note that G 0 ð1Þ ¼ nðtjsÞ, the mean density, f 0 G ðGð1ÞÞ ¼ mð1 À bÞ and f 0 G di ðG di ð1ÞÞ ¼ mb i , which leads to À @ nðtjsÞ @s where also Gð1Þ ¼ The quantity G 0 di ð1Þ is not the precursor concentration but is related to it. Fortunately we do not need it explicitly for our purposes, but we denote it by c i . With an independent source present, we have the equation À @G S @s ¼ S d ðsÞðf q ðGÞ À 1ÞG S or the mean value Now setting z = 1, and using the properties of G given above, we find for the second moment À @hnðn À 1Þi @s Similarly the subsidiary equation for G 00 di ð1Þ becomes À @ @s When there is a source present we use À @G S ðz; tjsÞ @s ¼ S d ðsÞ½f q ðGðz; tjsÞÞ À 1G S ðz; tjsÞ ð 82Þ From which, with f q ðGÞ ¼ 1 À v ðqÞ 1 ð1 À GÞ þ 1 2 v ðqÞ 2 ð1 À GÞ 2 þ Á Á Á, we get À @ @s For the case of only one neutron emitted per disintegration, we have v ðqÞ 1 ¼ 1 and v ðqÞ 2 ¼ 0 and S d is the neutron source strength. For a spontaneous fission source we would use the appropriate values of v ðqÞ n . Also, by definition, G 0 S ð1Þ ¼ N S and G 00 S ð1Þ ¼ hN S ðN S À 1Þi. The space and energy dependent versions of these moment equations, based on Eqs. (61)-(63) are given in Appendix E.
We will give some numerical examples of the mean and variance in due course. The quantity i.e. the relative standard deviation, sometimes called the sensitivity, will be a useful measure of the fluctuations in the neutron population. It is also worth noting that, for the point model, an equivalent set of equations for the mean and variance can be obtained from the forward form of the probability balance equation. For our purposes Eqs.
where KðtÞ ¼ 1= mk f is the generation time. Then if q depends only on the capture cross section k c ðtÞ, we have Any general form of startup procedure can be introduced here given qðtÞ but, for a ramp insertion of reactivity, we can write k c ðtÞ ¼ k c0 ð1 À ctÞ, whence the reactivity rate is The initial conditions are usually chosen such that all quantities are zero at t = 0. However, if the system starts from a subcritical state then more appropriate initial conditions can be used. In this work we assume that reactivity changes arise from control rod movement. We noted the need to calculate the relative standard deviation, R std , as a function of time to obtain t mat and in that connection it is instructive to consider the analytical form taken by R std in a simple case, namely, that of a point model with no delayed neutrons.
Thus r 2 ðtÞ As t ! 1, Eq. (91) which leads to almost the same limit as for the jump reactivity as t ! 1 but at a much faster rate. It is this type of behaviour that we observe in the low source problem with the control rods providing the ramp rate. The difference in the rate of convergence to the limiting value goes as e Àqt=K in the jump case and as e Àct 2 =2s in the ramp case. If these are compared numerically, one notes that for short times less than s or K, the jump case decreases more rapidly than the ramp case, on the other hand at later times the ramp rate converges to its limit much more rapidly than the jump one. There is no obvious physical reason for this. A more detailed discussion of the behaviour of r 2 = n 2 with time for jump and ramp variation may be found in Clarke et al. (1968).
A set of equations for the moments which includes space and energy effects has been derived by Humbert (2003) from Eqs. (33) and (41) without delayed neutrons. These have been solved numerically in Cartesian and curvilinear co-ordinates using discrete ordinate and multigroup methods. Humbert has also obtained the survival probability Pð0; t; Rj0Þ. However, from the discussion in Humbert's paper it would seem that numerical results for the mean and variance were only obtained for the point model. The energy and space dependent survival probability has also been studied numerically by Baker et al. (2009) using S N and multigroup methods and applied to a hypothetical criticality accident scenario. Again delayed neutrons are neglected, making this and the work of Humbert, unsuitable for low source startup studies. However, it would be straightforward to extend their works to include this aspect of the problem. Stacey (1969) has carried out three group calculations for mean and variance in a point model.
Finally we note that actual control rod movement does not follow a ramp variation but rather the form (Lamarsh, 1983) being the velocity of the control rod and H the core height. Use of this variation will not change our general conclusions in any significant way but clearly, in a practical situation, it is best to use the prescribed variation in the form qðtÞ. There is one other startup procedure of interest which is called 'jog and wait' or 'pull and wait' (Hurwitz et al., 1963a,b). In this, the control rod is withdrawn by a small amount and then kept in that position for a specified time; this procedure is repeated until the desired power is achieved.

Stochastic aspects in start-up
We have already discussed in Section 4 how the pdf enters our calculations for the safety probability 1-Q, where Pðn; tÞ ð 96Þ is the probability that the neutron population is less than n ⁄ . We have also seen how at some specific time t mat , the pdf becomes essentially a function of n= nðtÞ only and hence It is important to point out that this similarity solution is also obeyed by the precursor population such that P n nðtÞ ; t ! P c i c i ðtÞ ; t as t ! t mat or R std;n ðt mat Þ ! R std;c i ðt mat Þ i.e. the relative standard deviations of neutrons and precursors all tend to the same value for large time. This feature is illustrated graphically for one and six delayed neutron group solutions in Fig. 4. The initial reactivity is zero. Clearly the neutron and precursor curves converge to a single value as t increases; but to different values according the number of groups. We also note that the precursor fluctuations are much less than those of the neutrons except near the maturity times. The limiting value of R n ¼ r= n is 0.359 for one group of delayed neutrons and for six groups it is 0.480. The corresponding maturity times for one and six groups are, respectively, 55.2 and 42.3 s. The limiting behaviour of r= n was noted by Stacey (1969) and by Clarke et al. (1968). In order to calculate Q it will be necessary to solve numerically the Pál-Bell equations for the generating function, as discussed above, and then to invert that generating function. As noted in section 6, there are several techniques for inverting a generating function to recover the pdf and we will discuss two of the more useful of these below and in Section 10.

Saddlepoint method
The basic problem is to evaluate the quantity Q ðn Ã ; tjsÞ ¼ X n Ã À1 n¼0 P S ðn; tjsÞ; n Ã > 0 We have seen from Section 6 that the saddlepoint method gives the following result with the auxiliary conditions G S and the derivatives of G S are evaluated at z = z 0 , with z 0 given as the root of the equation The saddlepoint method (which is an approximation) has been shown to give accurate results except when Q is near unity, and become progressively more accurate as n ⁄ increases. When, as is usually the case, G S ðz; tjsÞ is the solution of a differential equation, the coupled set of differential-algebraic equations above must be solved by special methods (Mattsson and Soderlind, 1993). One exception to the algebraic-differential requirement is if we insert the values of z 0 ; G S and G 0 S into Eq. (102) directly and find the resulting value of n ⁄ . This does not enable us to calculate Q as a function of time, but it does enable us to calculate it as a function of the ratio nðtÞ=n Ã ðtÞ where n Ã ðtÞ is treated, not as an integer, but as a continuous function of time. This approach will be employed to get numerical results in Section 11. As t ! t mat , the ratio nðtÞ=n Ã ðtÞ ! nðt mat Þ=n Ã ðt mat Þ which leads to the desired safety probability Q.
For the case when an analytical expression is not available, which applies to most situations, we have to deduce the saddlepoint quantities directly from the differential equations. Thus we need equations for G S ðzÞ; G 0 S ðzÞ and G 00 S ðzÞ where primes denote differentiation with respect to z. Let us consider the case of no delayed neutrons when we may write from Eq. (59) in the quadratic approximation @Gðz; tjsÞ @s ¼Gðz; tjsÞ½k a ðsÞ À mk f ðsÞ þ and from Eq. (56) for a Poisson source that emits one neutron per disintegration, @G S ðz; tjsÞ @s ¼ SðsÞGðz; tjsÞ G S ðz; tjsÞ ð 104Þ with initial conditionsGðz; tjtÞ ¼ 1 À z and G S ðz; tjtÞ ¼ 1. For the saddlepoint method we find that the final conditions arẽ Gðz 0 ; tjtÞ ¼ 1 À z 0 and G 0 S ðz 0 ; tjtÞ ¼ 0. Thus from Eq. (102) the final condition becomesGðz 0 ; tjtÞ ¼ 1 À z 0 ¼ 1=ð1 þ n Ã Þ. Differentiating Eqs. (103) and (104) where the initial conditions areG 00 ðz; tjtÞ ¼ 0 and G 00 S ðz; tjtÞ ¼ 0. These equations will also allow the time-dependent transition coefficients k x ðsÞ to describe time varying reactivity. We digress for a moment to discuss the numerical algorithm for solving backward in time equations. For example in Eq. (103) the time variable must go from the final time back to a given value of s. To make this procedure simpler we have set w ¼ t À s and so, for example, Eq. (103) becomes @Gðz; tjt À wÞ @w ¼Gðz; tjt À wÞ Â mk f ðt À wÞ À k a ðt À wÞ À where 0 6 w 6 t. The algorithm is therefore to fix the final time t and solve the equation in w from zero to t. This gives the function value at t; now change t and repeat the process. The procedure is somewhat tedious as one has to solve the differential equations more times than in a normal initial value problem. But it is clear that the usual process, where the transition rates are independent of time, cannot be used as we are faced with the value of k x ðt À wÞ which is neither at t nor at w but rather at some intermediate point. This algorithm works and can be modified to improve efficiency.

Extension to one group of delayed neutrons
The generalisation to one group of delayed neutrons in the quadratic approximation is from Eqs. (53), (54) These equations are solved numerically. However, it will be noted in the saddlepoint method from Eq. (102) that, given n ⁄ , one must use a rootfinder to obtain the value of z 0 . This would require the whole problem to go through a Newton-Raphson, or similar procedure, for each iteration; clearly a mammoth task and one which would be very time-consuming. An alternative way forward is to use an algebraic-differential equation solver with Eq. (102) as an algebraic side condition. In their treatment of the forward equation of probability balance, which leads to a first order partial differential equation, Hurwitz et al. treat n ⁄ as a time-dependent variable and find the value of n ⁄ which corresponds to a given z 0 as defined by the final conditions. In a similar way, we proceed by selecting a value of z 0 and calculating the corresponding value of n ⁄ from Eq. (102). Thus Eq. (102) becomes a part of the equation set. This method has distinct advantages over that of Hurwitz et al. as will be demonstrated below and it works best for the case when t $ t mat , i.e. when the pdf has matured. This matter will be discussed in more detail later. It is possible of course to solve for fixed n ⁄ and varying time using the exact method of inversion due to Abate and Whitt to be described below, but it is very time consuming.

Extension to six delayed neutron groups
Following the procedure discussed in the last section, we use the modified generating functionGðzÞ ¼ 1 À GðzÞ and derive the associated equations for the saddlepoint method for I groups of delayed neutrons. The generating function is written as @Gðz; tjsÞ @s Let us now reduce the derivatives to a convenient form.
First therefore we need d dz To proceed further we now need which by using the procedure described above leads to Note that the term above may be accurately approximated by setting the denominator term, 1 À mb jGdj , to unity, leading to a term of Oðb 2 Þ, viz: Clearly, in the quadratic approximation, The equation forG 00 becomes therefore @G 00 ðz; tjsÞ @s ¼ k aG 00 ðz; tjsÞ þ k f f ðGÞ If we use the approximation noted in Eq. (127) then Eq. (130) becomes @G 00 ðz; tjsÞ @s ¼ k aG 00 ðz; tjsÞ þ k f f ðGÞ A very useful approximation which aids the numerical work is to use the relation, for a i << 1, This is particularly helpful when the saddlepoint method is extended to spatial problems and introduces a negligible error.
We now need the equations forG di ;G The final conditions arẽ Gðz; tjtÞ ¼ 1 À z;G 0 ðz; tjtÞ ¼ À1;G 00 ðz; tjtÞ ¼ 0 ð136Þ The above equations are then coupled with the saddlepoint equations to get Q ðn Ã ; tÞ $ 1 ffiffiffiffiffiffiffiffiffiffiffi ffi 2pr 0 p G S ðz 0 ; tjsÞ z n Ã 0 ð1 À z 0 Þ ð139Þ G S and the derivatives of G S are evaluated at z = z 0 , with z 0 given as the root of the equation Now instead of finding z 0 from Eq. (140) we observe that z enters the equations explicitly only through the final condition on, i.e.,Gðz; tjtÞ ¼ 1 À z. Thus if we insert the value of z into Eq. (140) it leads from G S and G 0 S to n ⁄ , whence we find r 0 and Q ðn Ã ; tÞ. Numerical calculations in a later section will illustrate this matter. Essentially, the idea is to fix z 0 and find n ⁄ rather than fixing n ⁄ and finding z 0 . As we have already noted, this works very well for the mature pdf and agrees with 'exact' results from direct inversion of the generating function by the method of Abate-Whitt (1992). The method does not work for the explicit time dependent case because then n ⁄ changes with each z 0 whereas it (n ⁄ ) must be kept constant. The algebraic-differential method must be used for this case (Mattsson and Soderlind, 1993) but for our problem this is not necessary.
A detailed numerical comparison of the results shown graphically in Hurwitz et al. (1963a,b) has been made with this work for both one delayed group (as used by Hurwitz et al.) and for six groups. Excellent agreement is obtained for the one group case (bearing in mind that it was necessary to read the values off the graph and also that Hurwitz used the zero prompt lifetime approximation). The six group curves, which are superimposed on the graphs in the next section, show once again the shortcomings of the one group approximation. This result is consistent with the one and six group work of Bell et al. (1963) in their interpretation of the Godiva experiment. As examples we give below three cases from the original Hurwitz paper. See also MacMillan (1970).
9.4. The curves of Hurwitz et al. (1963a,b) revisited We give here the curves in Hurwitz et al. (1963a,b) but calculated by our backward equations using a start reactivity of À0.5$.  As far as the value of the prompt neutron lifetime is concerned, Hurwitz et al. (1963a,b) propose a value based on the product ks ¼ 8 Â 10 À6 and a value of k ¼ 0:2 s À1 ; this leads to s ¼ 40 ls.
The value of k seems to be rather high since the mean value of the six group data in Wilson and England (2002) Fig. 4 (3000 on the curve should be 3500). The one delayed neutron group results are thick lines and those for six groups are thin lines shown in Figs. 5-7; we refer to these as 'Hurwitz curves'. It is clear that the six group results differ considerably from those of the one group case and that any practical calculations must use them; this conclusion is consistent with the work of Bell et al. (1963) which used one and six groups of delayed neutrons to interpret the Godiva experiment. In 1970, the Hurwitz curves were revisited by MacMillan (1970) who obtained accurate values of the generating function using the forward equation but also used a rather crude method to invert the generating function to recover the pdf. A computer code was developed, called SSB, which used six groups of delayed neutrons. We have been unable to obtain this code but the description of the algorithm in MacMillan's paper does not bode well for accuracy. The code SSB is based on the earlier code NDT3 . The numbers in the boxes to the right of the figures denote the values for source strengths of 250k to 5000k neutrons per second for one and six delayed neutron groups.  To demonstrate the difference arising from taking the start reactivity at different values we show Fig. 8. It is clear for the range of Q of interest (Q > 10 À10 ) that À1.0$ is acceptable. We chose À0.5 $ to be consistent with Hurwitz et al.

Comparison with Godiva
A crude comparison with the Godiva experiment as described in Bell (1963) and Hurwitz et al. (1963a,b) is shown in Fig. 9. The one group case is well below that of the six group one and the experimental results shown in Fig. 1 of Bell et al. (1963) lie very close to the six group curve as does our theoretical result (it was necessary to read the values from the figure and so about 20% uncertainty is expected). Also the one group result, although a poor comparison with experiment, is close to the one group results obtained by Bell et al. and so the results are consistent. The description of the experimental events that lead to the results in Fig. 1 of Bell et al. state that ''Godiva was rapidly brought from a very subcritical state to 70 cents above delayed critical and the time at which the fission level reached some fiducial level was measured". There is a background source of 90 n/s. Later in the article it is suggested that the fiducial level is 2000 neutrons. We have taken our time of measurement equal to the maturity time as has been done by Bell et al., i.e. to the asymptotic distribution (as they term it). The comparison with experiment is reasonable and the expression 'very subcritical' is open to question; we took this as À15$ and the term 'rapidly' was taken as 20$/s. We also found that the value of Q ðn Ã ; t mat Þ is insensitive to the prompt neutron lifetime when its value is less than 45 ls as noted in Fig. 14 and Table 2. There should be error bars on the experimental points but these have not been given in any reports or articles seen by the authors.
10. Numerical solution directly from the differential equations using the Abate-Whitt method In general, it is unlikely that an analytical solution for the generating function will be available and so we must start from the differential equations for the generating functions, solve them for a range of z and then reconstruct the pdf. Thus returning to Eqs.
To obtain the probability distribution function from the above backward equations we go to Eq. (A24) in Appendix A. But now, instead of an explicit expression for Gðre i# j;k Þ, we have to use the solution of the differential equations. As a check on the method and its viability, we consider the case of the quadratic approximation with no delayed neutrons and time-dependent transition rates. The resulting equations are with t À s ¼ w and k x ðsÞ ¼ k x ðt À wÞ. dG 1 dw where G 1 ¼ G 1 ðj; k; wÞ and G 2 ¼ G 2 ðj; k; wÞ. The initial conditions are G 1 ðj; k; 0Þ ¼ 1 À r cos # j;k and G 2 ðj; k; 0Þ ¼ Àr sin # j;k . The symbol a ¼ mk f À k a . The source generating function G S is given, for v ðqÞ n ¼ d n;1 ,Ĥ 1 ¼ ÀG 1 andĤ 2 ¼ ÀG 2 , by dG S1 dw ¼ ÀSðG 1 G S1 À G 2 G S2 Þ and dG S2 dw We have run several cases in which an increasing number of multiplicities v n are included. The conclusion is that increasing the number to the full seven, compared with the quadratic approximation, for reactivities less than 3$ leads to differences in Q of the order of 0.01%. The time of execution in going from two to seven multiplicities is about a factor of 10. Thus, when feasible, it is acceptable and practical to use the quadratic approximation. The effect of increasing the number of delayed neutron groups is, however, likely to be much more important.

Results via the saddlepoint method
We have shown in Section 6 that the saddlepoint method leads to a value for the cumulative pdf in the form We have discussed how G S and its derivatives are to be obtained and we now wish to illustrate the method by numerical results. We commence by assessing the accuracy of the saddlepoint method. Hurwitz et al. (1963a,b) assessed the accuracy of their approach using the saddlepoint method by comparing it with some numerical inverse Laplace transforms. We will do it by comparison with the exact inverse of a generating function as devised by Abate and Whitt (1992). We choose as data S=k ¼ 5000; R=k ¼ 0:2$; K ¼ 45 ls and one group of delayed neutrons. This defines the ramp insertion of reactivity and is one of the cases considered by Hurwitz et al. We also assume two situations; one is when the initial reactivity just before the ramp is zero $, i.e. it is just at delayed critical. A more likely situation is when the system is highly subcritical initially and for this we assume that qð0Þ ¼ À5$. As we have seen above, the effect of sub-criticality below À1$ has little effect on the subsequent ramp behaviour or, to put it otherwise, a sub-criticality of À15$ can be equally well represented by a sub-criticality of À1$ for all the effect it has on the subsequent results. Hurwitz et al. also found this behaviour. Fig. 10 shows the value of Q ðn; t 0 Þ vs n= nðt 0 Þ for different values of t 0 from 50 s to 56 s; in fact 56 s corresponds to the time at maturity, so the figure shows the approach to maturity. The thick lines are from the exact inverse and the thin ones from the saddlepoint method. Numerical results show the values to be within 0.5% of each other and so we may conclude that for all practical purposes the saddlepoint method is a reliable tool for studying the low source problem. Moreover, it is readily adapted, via the backward equations, to describe energy and space dependent effects; the forward equation cannot do this. Fig. 11 shows a curve similar to that in Fig. 10 but with the initial reactivity at À5$. The final time in Fig. 11, which corresponds to the maturity of the pdf, is 350 s. Of course, from a safety point of view, the time during which the system is vulnerable to stochastic surges does not start until the system becomes critical; thus in  n/<n(t)> Fig. 10. Q(n, t) cumulative pdf as a function of n/nav. Exact vs sadddlepoint q(0) = 0.0$. the case of Fig. 11 the first 350-56 = 294 s is 'dead' time. The exact method using the Abate-Whitt algorithm required so many differential equations that the computing time amounted to several hours for the larger values of n ⁄ . However, the saddlepoint method has no such problems and this demonstrates its versatility and utility. In Fig. 11 the thick lines denote the exact result and the thin ones those of the saddlepoint method; the results are to within 0.7% of each other. We also show the mature value of Q in Fig. 11 for the Gamma pdf. This is reasonably accurate for Q > 10 À5 , but is not always so reliable and, as we have a very accurate value via the saddlepoint method over a wide range of Q values, the Gamma pdf will not be used. It may however prove useful as a guide when the problem is energy and space dependent. This matter will be discussed in another place. It is important to note that the accuracy given above for the saddlepoint method is based on a logarithmic scale. On a linear scale the accuracy is around 7% but, because we are only concerned with small values of Q ($ 10 À6 ), it is more realistic to define the error in terms of log 10 Q. In that case the fractional error is around 0.7%.
Because Figs 10 and 11 use a time-dependent form of the saddlepoint method as described in Section 10, we feel that a more detailed explanation of the procedure is needed. Thus, having calculated the average density hnðt 0 Þi at time t 0 < t mat , we must fix the final condition onGðz; t 0 jsÞ for s ¼ t 0 by specifying z 0 where 1 À z 0 ¼ 1=ð1 þ kÞ for k ¼ 1; 2; :0:0:10 k . We then solve the equations forGðz; t 0 jsÞ up to s ¼ t 0 then insert values ofGðz; t 0 jsÞ, etc and z 0 into G S ðz; t 0 jsÞ and hence calculate From this we obtain n Ã ðt 0 Þ=hnðt 0 Þi and we also have the corresponding value of Q ðn Ã ðt 0 Þ=hnðt 0 Þi; t 0 Þ. So each curve in Figs 10 and 11 are for a fixed t 0 and varying n Ã ðt 0 Þ=hnðt 0 Þi.
To further illustrate the limitations of the Gamma pdf we show Figs 12 and 13 (referred to as the Hurwitz curves) for two values of the reactivity insertion rate and an initial subcriticality of À1$. Fig. 12, which has a reactivity insertion rate of R ¼ 0:2k $ s À1 and a range of source strengths, demonstrates clearly that the Gamma pdf fails for all source strengths less than 5000k n s À1 . On the other hand, in Fig. 13 for which R ¼ 0:01k $ s À1 , the Gamma pdf is close to the saddlepoint method for source strengths greater than 1000k n s À1 . It is quite possible, therefore that, if we are concerned with very slow startup scenarios, the Gamma pdf may be helpful. This is especially true if space dependence is included for then it will only be necessary so solve the space dependent moment equations for the mean and variance. We also note that in practice the desired value of Àlog 10 (Q) is unlikely to be greater than 8 and so Fig. 13 indicates that the Gamma distribution is valid over a wide range of source strengths for reactivity insertion rates less than 0:01k $ s À1 .
To illustrate the influence of the generation time K on the Hurwitz curves, we show Fig. 14 for one and six groups of delayed neutrons obtained via the saddlepoint method. The smallest value of Q normally of interest is 10 À8 , or Àlog 10 ðQ Þ ¼ 8. In the range 0:45 < KðlsÞ < 450, the values of n=n Ã at Q ¼ 10 À8 are as shown in Table 2. It is clear that the value of Q, or conversely the revised source strength, is sensitive to K, but not overly so. Hurwitz et al. (1963a,b) use the zero prompt neutron lifetime approximation and we can now check on the accuracy of that approximation. It is readily seen from Table 2 that, as K decreases, the value of n=n Ã tends to a constant value; 32 for one group and 159 for 6 groups.
We have used K ¼ 45 ls in our calculations and this gives close agreement with the figures in the Hurwitz series of papers; Table 2 shows why this is so. Much more important is the sensitivity to the number of delayed neutron groups. For example, from Table 2, the one group model underestimates the magnitude of the revised source by around a factor of five, which is a serious error and confirms the need to use a complete set of delayed neutron groups. We have not tried using an 8 group set of delayed neutrons but do not expect it to change our conclusions in any significant manner.

Space dependence of the generating function
As we have seen in Section 8, Eq. (67), the one speed diffusion theory version of the equation for the generating function, may be written in the quadratic approximation and with no delayed neutrons for a homogeneous medium as À @Gðz; t; Rjr 0 ; sÞ @s À vDr 2 0G ðz; t; Rjr 0 ; sÞ The final condition isGðz; t; Rjr 0 ; tÞ ¼ ð1 À zÞDðr 0 ; V r Þ, where D ¼ 1 if r 0 2 V r ; else D ¼ 0: When a source is present, we have ðR ¼ V r Þ, G S ðz; t; RjsÞ ¼ exp À where G S ðz; t; RjtÞ ¼ 1. It is important to note that the probability Q ðn Ã ; tjsÞ does not depend explicitly on position but is influenced by it through the geometry of the system and the position of the source via Eq. (156).

Expansion in eigenfunctions
There is a number of different methods available to solve Eq. (154) numerically, but a simple one is to use an expansion in an orthonormal set of functions appropriate to the system geometry. Thus let us assume that we have a homogeneous bare reactor in which the time dependent flux may be written /ðr; tÞ ¼ X n w n ðrÞT n ðtÞ ð 157Þ where the eigenfunction w n ðrÞ satisfies the Helmholtz equation r 2 w n ðrÞ þ B 2 n w n ðrÞ ¼ 0 w n ðrÞ is zero on the extrapolated boundary and obeys the orthonormality condition ðw n ; w m Þ ¼ N m d n;m , the integer n ¼ ðn 1 ; n 2 ; n 3 Þ. Let us now expand the generating function as Gðz; t; Rjr 0 ; sÞ ¼ X n w n ðr 0 Þg n ðz; t; RjsÞ ð 159Þ Inserting this expression into Eq. (154), multiplying by w m ðrÞ and integrating over the reactor volume, we find (suppressing the arguments of g n ð. . .Þ) À @g m @s subject to the final condition g m ðz; t; RjtÞ ¼ ð1 À zÞ where V is the region in which the probability Q is desired; in our case this would normally be the whole system. The equation for G S ðz; t; RjsÞ then becomes @ @s G S ðz; t; RjsÞ ¼ G S ðz; t; RjsÞ X n g n ðz; t; RjsÞS n whereS n ¼ R R dr 0 Sðr 0 Þw n ðr 0 Þ and G S ðz; t; RjtÞ ¼ 1 Inversion of G S ðz; t; RjsÞ then yields the probability distribution. If the source is localised, for example resides in a very small volume which may be approximated by a delta function, we have Sðr 0 Þ ¼ S 0 dðr 0 À r p Þ andS n ¼ S 0 w n ðr p Þ. We now have a set of coupled, first order differential equations for the expansion coefficients g n . If the saddlepoint method is to be used then we have to find equations for g 0 n and g 00 n . Convergence of the series (159) is very important and Bell (1965) has suggested that the coefficients of higher spatial modes will decrease rapidly with time compared with the fundamental. In that case we can set n = 0 in Eq. (160) to get À @g 0 @s with For a spherical reactor of radius R, we have w n ðrÞ ¼ ffiffi 2 R q sinð npr R Þ=r, the ratio D 000 ¼ 0:9702 . . .. This result suggests that the multiplicity v 2 is effectively reduced to 0:9702v 2 , implying that space dependence very slightly reduces the fluctuations. Bell and Lee (1976), using one speed transport theory, obtained an essentially identical expression for D 000 . From an historical point of view, we note that Schroedinger (1945) using integral transport theory, also found a factor which depended on D 000 . Using only the fundamental mode may well be accurate in a small system like Godiva where we would not expect the pdf to vary with position. However, in a large power reactor it is likely that higher harmonics may be excited and the pdf may be influenced by the spatial instabilities, analogues of which occur in the deterministic regime (Bell and Glasstone, 1970). We will report on space dependence in more detail in a forthcoming publication but we can say that its inclusion is vital if the source is localised. An example of a slab reactor in which the source is at three positions in the core is shown in Fig. 15. The reactor parameters are, with an initial reactivity of zero and one group of delayed neutrons, width of core a ¼ 100 cm; ¼ 5000k ns À1 ; R ¼ 0:2k $ s À1 ; k ¼ 0:08519 s À1 ; K ¼ 45 ls The position a/2 is the core centre and a/8 very close to the edge. It is clear that the shape and position of the source plays an important role in determining the value of Q. One aspect that deserves future investigation is to see how close the exact spatial result is to a point model with the source 'spatially importance weighted'. There is some information of vital importance that may be obtained from Fig. 15; namely the factor by which we must multiply the source strength ð5000k ¼ 426 n=sÞ to achieve a given Q value. Table 3 shows the case for Q ¼ 10 À6 . It is interesting to note that a source located at the mid-point of the reactor has a very similar effect to one that arises from the point model. It might be asked why the uniform case and the point model are not the same, this is because in the space dependent model the shape of the generating function is included, but in the point model it is not.

Energy dependent Pál-Bell equation and its influence on the point model results
It is now necessary to investigate the importance of energy dependence and how to incorporate it into the calculation of the generating function and ultimately the safety factor 1-Q. The necessary equation has already been given by Eq. (64) but we write it more explicitly as À @Gðz; t; Rjr 0 ; v 0 ; sÞ @s ¼ Àðk a ðr 0 ; v 0 ; sÞ þ k s ðr 0 ; v 0 ; sÞÞGðz; t; Rjr 0 ; v 0 ; sÞ The final conditions arẽ G di ðz; t; Rjr 0 ; tÞ ¼ 0 and G S ðz; t; RjtÞ ¼ 1: It is useful to note, that if the capture and fission cross sections are both proportional to 1/v and vDðvÞ and Hðx; y; vÞ are independent of velocity, thenGðz; t; Rjr; E; sÞ ¼Gðz; t; Rjr; sÞ because the scattering terms cancel. In general this is not true and so we must consider the explicit effect of energy dependence (however the 1/v case may be used as a benchmark). In order to put Eq. (166) in a more convenient form, let us convert it from velocity to energy, leading to À 1 v 0 @Gðz; t; Rjr 0 ; E 0 ; sÞ @s ¼ ÀðR a ðr 0 ; E 0 ; sÞ þ R s ðr 0 ; E 0 ; sÞÞGðz; t; Rjr 0 ; E 0 ; sÞ þ r r 0 :Dðr 0 ; E 0 ; sÞr r 0G ðz; t; Rjr 0 ; E 0 ; sÞ þ Z dE 0 R s ðE 0 ! E 0 ; r 0 ; sÞGðz; t; Rjr 0 ; E 0 ; sÞ This equation may be converted to multigroup form as shown below, with g denoting the energy group ðE gÀ1 ; E g Þ. À 1 v g @G g ðz; t; Rjr 0 ; sÞ @s ¼ ÀðR a;g ðr 0 ; sÞ þ R s;g ðr 0 ; sÞÞG g ðz; t; Rjr 0 ; sÞ þ R f ;g ðr 0 ; sÞ þ r r 0 :D g ðr 0 ; sÞr r 0G g ðz; t; Rjr 0 ; sÞ þ X G g 0 ¼1 R s;g!g 0 ðr 0 ; sÞG g 0 ðz; t; Rjr 0 ; sÞ The details are given in Appendix B, but it is useful to note that the formalism is general enough to deal with a mixture of source types, e.g. a combination of single emission ða; nÞ sources of various neutron energies and a spontaneous fission source with prescribed emission characteristics, all at various positions in the core. A measure of the influence of energy dependence can be assessed by repeating our point model calculations in multigroup form. Thus we consider Eq. (172) without spatial dependence, viz: À @G g ðz; t; RjsÞ @s ¼ Àðk a;g ðsÞ þ k s;g ðsÞÞG g ðz; t; RjsÞ þ k f ;g ðsÞ with k x;g ¼ v g R x;g and @G di ðz; t; RjsÞ @s ¼ k iGdi ðz; t; RjsÞ À k iGp;i ðz; t; RjsÞ ð 174Þ and @G S ðz; t; RjsÞ @s ¼ X G g¼1 S g ðsÞG g ðz; t; RjsÞ G S ðz; t; RjsÞ ð 175Þ The final conditions at t ¼ s areG g ðz; t; RjtÞ ¼ ð1 À zÞ DðE g ; U v Þ;G di ðz; t; RjtÞ ¼ 0 and G S ðz; t; RjtÞ ¼ 1 and since all values of E g lie within the range U v , DðE g ; U v Þ ¼ 1. If, however, we wish to know the number of neutrons in the energy range DE g y , then the final condition would becomeG g ðz; t; RjtÞ ¼ ð1 À zÞd g;g y . As an example, let us consider a two group approximation such that À @G 1 ðz; t; RjsÞ @s ¼ Àðk a;1 ðsÞ þ k s;1!1 ðsÞ þ k s;1!2 ðsÞÞG 1 ðz; t; RjsÞ þ k f ;1 ðsÞ þ k s;1!1 ðsÞG 1 ðz; t; RjsÞ þ k s;1!2 ðsÞG 2 ðz; t; RjsÞ À k f ;1 ðsÞH 1 ðG p ;G di Þ ð176Þ and À @G 2 ðz; t; RjsÞ @s ¼ Àðk a;2 ðsÞ þ k s;2!1 ðsÞ þ k s;2!2 ðsÞÞG 2 ðz; t; RjsÞ þ k f ;2 ðsÞ þ k s;2!1 ðsÞG 1 ðz; t; RjsÞ þ k s;2!2 ðsÞG 2 ðz; t; RjsÞ À k f ;2 ðsÞH 2 ðG p ;G di Þ ð177Þ Detailed numerical results on multigroup problems will be reported in a forthcoming publication, but an initial assessment based on the two group model suggests that the value of Q, i.e. the goal of the exercise is not very sensitive to energy, especially over the range ðQ : 10 À8 ; 10 À5 Þ but this result should not be taken as definitive until more detailed multigroup calculations have been carried out. Some work on the effect of energy dependence on the extinction probability has been carried out by Saxby et al. (2016).

A practical problem
The practical aspects of this work involve start up problems and 'just safe' deterministic values of source strength S m and reactivity rate R m . We will now pose a simple start up problem in which the power is allowed to rise to a prescribed value (say 1 MW) and then control rods are tripped and shut down occurs rapidly. We calculate the amount of energy deposited in the reactor during this tran-sient and establish that no damage is done. But any further energy deposition would be unacceptable. Knowing the source strength S m n/s and the reactivity insertion rate (assumed linear ramp) R m $/s, we can then define the 'just safe' condition in a deterministic sense. Now we ask for the probability that this combination of source and reactivity (S m , R m ) is statistically safe. Namely, what is the probability that the source strength is so low as to lead to a significant fraction of neutrons being below some fiducial level, e.g. below the maturity level which defines the boundary between the stochastic and deterministic regions. As we know, if too large a fraction of neutrons is below this level then there is a non-negligible possibility that the associated neutron density will be dangerously low (due to fluctuations) and hence cause a more severe transient than the deterministic calculation would predict.
The data which define our problem are:  In Fig. 16 we show the relative standard deviation of the neutron population defined by rðtÞ= NðtÞ for one and six groups of delayed neutrons. Note that the precursor and neutron values go initially in different directions but all eventually coalesce at the 'maturity point'. i.e. where the neutron population has such a small relative standard deviation that it may be regarded as deterministic. The maturity time, asymptotic relative standard deviation and power for one and six groups are, respectively, (299.3 s, 0.210, 39.8 W) and (290.2 s, 0.267, 5.06 W). Maturity is seen to develop very rapidly as all the curves coalesce.
We shall regard a source strength of 1000 n/s and ramp rate of 0.02 $/s as just safe values (S, R). The problem now is to see by how much one should increase the source strength, or reduce the ramp rate, to achieve a small probability that the neutron density is too low and is likely to lead to an unacceptable surge. To calculate this, we return to the saddlepoint method described above. For more generality we have calculated the cumulative pdf for the additional values of S = 500 n/s, 2000 n/s and 4000 n/s and to be more realistic we start the system from a subcriticality of 5$. Also, to show the error in using one group of delayed neutrons, we show the six group values. They are shown graphically in Fig. 17 and in Table 4. We note that one group values significantly underestimate the values of Q and hence are non-conservative, i.e. Q ð1 gpÞ << Q ð6 gpÞ. We believe therefore that the results of Hurwitz et al. (1963a,b) should be viewed with some caution and Bell et al. (1963) and MacMillan (1970) have shown the significance of six groups of delayed neutrons and arrive at conclusions that are similar to ours. In the region of importance for low source startup, i.e. Q < 10 À5 , six groups of delayed neutrons are always required. Table 4 also shows in brackets the magnification factor for a ramp rate of 0.007 $/s. It is observed that in order to decrease Q we must either use a stronger source or a smaller ramp rate. In the work of Hurwitz et al., in order to obtain the new source strength a relationship of the form n Ã ðt mat Þ= nðt mat Þ ¼ S m =S is used (see Section 4); in this case it is easier to simply calculate Q directly from the figure for each 'just safe' set (S, R), rather than from the 'Hurwitz curve' which displays nðt mat Þ=n Ã ðt mat Þ vs Àlog 10 ðQ Þ.

Summary and conclusions
It is many years since Hurwitz and co-workers published their classic work on low source startup (Hurwitz et al., 1963a,b). Their work arose from a requirement to establish the minimum strength of the source to use to avoid any stochastic effects leading to 'rogue' transients. Too low a source strength and the possibility would arise of rogue transients caused by first persistent chains. The original work performed by Hurwitz used a point model with one group of delayed neutrons and the forward form of probability balance. Although the numerical results obtained by Hurwitz are useful, they are by no means accurate; not because of any innate error but because they only used one group of delayed neutrons rather than six. To obtain numerical results, it was necessary for Hurwitz to use the properties of the precursors rather than the neutrons to obtain the stochastic results, mainly because in the early 1960's computing power was, by today's standards, primitive. Also the point model did not allow for the possibility of the  source being localised unless some separate importance weighting was carried out. The work was also one speed and no indication is given as to how this may influence the outcome. Hurwitz, however, did use five of the multiplicity coefficients v n , but our calculations show that only terms up to quadratic are necessary which can speed up numerical work significantly.
To remedy the neglect of space and energy effects, it is the view of the authors that it is no longer advisable to use the forward form of the probability balance equation because there is no practical way to include these additional variables. Stacey (1969) has developed a formalism using the forward equation for dealing with space and energy but has only shown its viability for calculating moments of the distribution and not the probability distribution itself. The only comprehensive study of these matters must be based on the Pál-Bell equation which uses a process that is backward in time, i.e. one starts from the final time and then works back to the intial one. Such an approach, developed initially in cosmic ray studies (Janossy, 1950), was extended by Pál in the 1950's to describe the behaviour of neutrons and precursors in a nuclear reactor. In 1965 Bell also published a paper on this matter. The method is powerful because it is based upon tracing the history of a single neutron from its birth in fission or from a source to its ultimate demise by capture or leakage. Having got the single neutron pdf it is relatively straightforward to use that to construct the actual pdf when a source is present; it is analogous to the conventional Green's function technique (but not identical). Also the formalism allows the introduction of space and energy because the initial condition on the single neutron allows it to have a specific energy and a prescribed position. It is this procedure that we have used here and have been able to replicate the Hurwitz work in a much simpler fashion. One very important aspect of the earlier work has in fact been retained, namely the use of the saddlepoint method to obtain an approximate inversion of the generating function. By comparison with an exact inversion method, we have shown that the saddlepoint approach is very fast and very accurate and can be used with confidence for space and energy problems. It is this fact that makes the calculations feasible for complex geometries and energy ranges. Use of the exact inversion techniques would lead to execution times being several orders of magnitude greater. We have given a range of numerical results to illustrate the startup procedure and of the probability that, for a given source strength and rate of reactivity introduction, the probability of a rogue transient between startup time and the time at which the system becomes deterministic is a specified value, e.g. 10 À8 . We have also developed a method which enables one to calculate by how much the just safe source strength must be multiplied to obtain a new source strength which results in a given probability that no rogue transient will arise.

Acknowledgement
This work was supported by Rolls Royce Ltd and the authors thank the Company for giving them permission to publish the paper. The authors would also like to thank Dr Chris Cooling for some stimulating discussions. Dr M.D. Eaton would like to thank EPSRC for their support through the following Grants: ''Adaptive Hierarchical Radiation Transport Methods to Meet Future Challenges in Reactor Physics" (EPSRC Grant No.: EP/J002011/1) and ''Nuclear Reactor Kinetics Modelling and Simulation Tools for Small Modular Reactor (SMR) Start-up Dynamics and Nuclear Criticality Safety Assessment of Nuclear Fuel Processing Facilities" (EPSRC Grant No.: EP/K503733/1). Finally, MMRW would like to record his indebtedness to Professor Imre Pazsit for always being there to answer his somewhat arcane questions on stochastic processes.

Data statement
In accordance with EPSRC funding requirements all supporting data used to create figures and tables in this paper may be accessed at the following URL: https://doi.org/10.5281/zenodo.193040.
Appendix A. Inversion of generating function by the Abate-Whitt method In general, it is unlikely that an analytical solution for the generating function will be available and so we must start from the differential equations for the generating functions, solve them for a range of z and then reconstruct the pdf. For the probability that there are no neutrons present, i.e. Pð0; tÞ we need only set z = 0 in the Pál-Bell equation and solve for Gð0; tjsÞ. For the first few values of Pðn; tÞ; n ¼ 1; 2; 3 . . . we may use the relation Pðn; tjsÞ ¼ 1 n! @ n Gðz; tjsÞ @z n z¼0 ðA1Þ and apply the operation directly to the Pál-Bell equation. However for large values of n, this is impractical and we must use the method described below. Thus returning to Eqs. (53), (54) and (56) of the text, we set z ¼ re i# and considerGðz; tjsÞ ¼Gðr; #; tjsÞ although we will retain the symbol z for convenience. Splitting the solution into real and imaginary parts we find Gðz; tjsÞ ¼ G 1 ðz; tjsÞ þ iG 2 ðz; tjsÞ and G dj ðz; tjsÞ ¼ G j1 ðz; tjsÞ þ iG j2 ðz; tjsÞ ð A2Þ and G S ðz; tjsÞ ¼ G S1 ðz; tjsÞ þ iG S2 ðz; tjsÞ ð A3Þ We insert these expressions into the Eqs. (53), (54) where q n ¼ ðÀ1Þ n v n =n! and  In the subcritical region, by definition, the associated average density is independent of time as is the variance. In addition, the macroscopic cross sections will be independent of time thereby rendering the solution of the backward generating function equations much easier to solve. The generating function Gðz; tjsÞ ¼ Gðz; t À sÞ ¼ Gðz; wÞ. From this we note that the source generating function equation If it is assumed that the system has been sub-critical for a long time then one may set s ¼ À1, and G S ðz; 1Þ G S ðzÞ ¼ exp ÀS 0 Thus G S ðzÞ is independent of time, but the single neutron generating function is not. For the case of no delayed neutrons it is straightforward to solve the equation for Gðz; wÞ in the quadratic approximation and we find with a ¼ k a À mk f where q 0 ¼ 1 À k a = mk f and g 0 ¼ 2S 0 =k f v 2 . The moments are v 2 k f 2a , which leads to the variance We can therefore write G S ðzÞ as For a case where there are no delayed neutrons, G S ðzÞ would provide the initial condition for the transient studies. In all practical cases, there are delayed neutrons and we have shown that to include these, D5 may be written in the same form but with a modified value of g 0 which incorporates the influence of delayed neutrons. The details of this derivation will be published elsewhere. In practice, the system is usually at a high subcriticality before startup and it is a very good approximation to assume that there are no neutrons, thereby enabling us to set G S ðzÞ ¼ 1. We have done this in our calculations in the text. We may also arrive at similar results by starting our calculation at an earlier time with a constant sub-critical reactivity until a steady state has been reached.

Appendix E. The mean and variance with space and energy dependence
In section 8, we derived the mean and variance for a one speed, point model. Here we wish to extend that to the more general case. Many of the manipulations are very similar because the neutron fission spectra for both prompt and delayed neutrons are assumed to be isotropic.