Spatial Effects in Low Neutron Source Start-up and Associated Stochastic Phenomena

This work concerns the calculation of the neutron source strength necessary to start up a nuclear reactor such that the likelihood of an undesirable stochastic transient is reduced to a specified value (e.g. 8 10 ). We extend our earlier point model work on low source calculations to include the spatial variation of the neutron source. Results for the source multiplier for a given safety factor are obtained for slab, cylindrical and spherical systems. The spatial term in the Pál-Bell equation is dealt with by Chebyshev-Gauss-Lobatto collocation methods and this enables an extrapolation distance to be included, thereby simulating a reflector. Results are given for a range of system sizes, and corresponding source multipliers for safe source determination are obtained. The saddlepoint method is used to invert the generating function. In addition to the low source calculations, we have also tested the collocation method on the survival probability in a sphere which demonstrates excellent convergence. We also comment on the usefulness of the Gamma pdf for spatially dependent problems. For clarity of presentation, some of the detailed mathematical work is relegated to Appendices.


Introduction
This paper extends our earlier work , henceforth referred to as I, on low source effects during start-up to include a detailed treatment of the spatial variation of the source on the calculation of the safety factor Q, where Q is the likelihood of an undesirable stochastic transient, and the resulting safe source magnitude. We will briefly review the main points of the earlier work before presenting the detailed mathematical treatment based on the Pál-Bell equations for dealing with space dependence.
The safe start-up of a nuclear reactor depends upon the presence in the core of a steady neutron source. As well as the fixed source, there also exist intrinsic natural neutron sources from spontaneous fission, cosmic rays, photo neutrons, fission products, etc, and it is necessary to account for these when assessing the overall effect. From an operational point of view, it is important to assess the strength of the natural sources to see if they will be sufficient in magnitude to ensure safe stochastic start-up without the addition of an extraneous source. The most important case for source evaluation is that of a reactor starting up with fresh, un-irradiated fuel because then the natural sources will be at a minimum. The theory developed below will be a powerful tool in that respect.
To understand how the reactor behaves as the control rods are withdrawn or reactivity is increased in some other way, it is necessary to solve the traditional equations of reactor kinetics. A full space-energy-time dependent case should be carried out if possible, but initially a point model approach will be sufficient for guidance. Now, in conventional reactor kinetics, it is implicitly assumed that the neutron source strength and the associated neutron density are high enough to reduce any statistical fluctuations arising from the underlying random processes to negligible proportions. However, for a sufficiently low source strength this may not be the case and the possibility of large fluctuations can arise. This matter is better understood 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, may assume that the rod has 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 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. Starting a reactor with no source, or a source of unknown strength, is known as a 'blind start' and is to be avoided for the reasons explained above (Shaw, 1969).
It is useful to explain physically why the stochastic behaviour differs from the deterministic one. First we note that the chain reaction consists of many independent chains of neutrons which are initiated and, depending on the multiplying properties of the medium, will either die out from leakage or capture or will multiply due to fission. But even in the case of multiplication, there is the possibility that the chain will eventually die out or, in the language of statistics, become extinct. It is only when a persistent chain arises, i.e. a chain that continues to multiply and increase the neutron number indefinitely, that a genuine supercritical transient will arise. With a weak source, such persistent chains can take several seconds to develop, during which time the reactivity can increase to a high value so that when the persistent chain does appear, a very rapid increase in power and energy release arises. This energy release can be several orders of magnitude greater than that predicted by deterministic reactor kinetics (Hansen, 1960, Cooling et al, 2016. If the source is strong enough to make fluctuations unlikely, it means physically that the independent chains are overlapping in time and hence the likelihood of the creation of a single dominant persistent chain is remote. 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 (1963 a, b, c), MacMillan (1960MacMillan ( , 1970, and that of Bell et al (1963) to invert the associated generating function. The accuracy of the saddlepoint method has been assessed by comparison with an exact inversion formula (I). We include delayed neutrons, both one and six groups, and come to the conclusion that six groups are essential for an accurate probability evaluation. Results are presented in graphical form and a computer code will eventually be available to enable the user to assess quickly and accurately the low source behaviour with a spatially movable source. A computer code, CALLISTO, based on the multi-group point model described in Cooling et al (2017) is already available. It has recently been pointed out to the authors (Pazsit, 2017) that the saddlepoint method was first used by Jánossy (1946Jánossy ( , 1950 and Jánossy and Messel (1950) to invert generating functions in connection with electron-photon cascades.
We have also found that the Gamma pdf is a reasonable approximation to the actual pdf and can be used for guidance on the influence of a low source without having to calculate the exact pdf. It is not, however, reliable for cases where the value of the pdf is less than 4 5 10 10 − −  , except for very low reactivity insertion rates and strong sources. Nevertheless, it may well be useful as a guide when complex arrays of sources are present both in position and energy because it reduces the computational time significantly.

General theory
The full details of the theory based on the Pál-Bell equation (Pázsit and Pál, 2008) are given in I, but we summarise here the main points, especially insofar as they concern spatial variation. The probability distribution and related generating function are defined below and the equation obeyed by the generating function or, to be precise, the complement of the one speed single neutron generating function ( , , , , ) 1 ( , , , , ) g x R t s g x R t s = − y r y r  , may be written for a homogeneous medium in the diffusion approximation as (Pázsit and Pál, 2008), where the delayed neutron generating functions ( , , , , ) i h x R t s y r  obey ( , , , , ) ( , , , , ) ( , , , , ); 1,2,..
The term ( ) involves the multiplicity processes due to fission and delayed neutrons, it is a complicated function and can be found in Appendix B. Coupled to eqns (1) and (2) is the equation which incorporates an independent source into the calculations. This is written In eqn (3a) ( , ) d S s r is the number of disintegrations of the source per second, each disintegration leading, with a given probability, to n neutrons. ( ) q f x is the associated probability generating function defining the probability that n neutrons are emitted. For simplicity, in this work we assume that one disintegration leads to one neutron, when eqn (3a) becomes To be more explicit regarding the generating functions, we have in standard statistical notation, The generating functions are , 0 , In other words the chain reaction can be triggered off by one prompt neutron, or by the neutron emitted by the decay of any precursor. The boundary conditions will be discussed in section 4.3. In our previous work we have solved eqn (1) for the point model in which the spatial operator was ignored. Here we will explicitly include the spatial variation which will require a solution of eqns (1) and (2) and subsequent quadrature in eqn (3 a,b).

Collocation
In order to solve eqn (1) numerically, it will be necessary to convert the term 2 g ∇  into finite difference form. To do this we will use Chebyshev-Gauss-Lobatto (CGL) collocation (Kopriva, 2009) which enables the first and second derivatives to be written in a very simple and convenient form and also allows the boundary values to correspond to a collocation point. The details of this process for one dimensional systems (slab, sphere and cylinder) are given in Appendix A, but the outcome is that for a system of the form with two boundaries and the following boundary conditions, we can write after the transformation where the collocation points are (11) as they stand refer to a slab, but if we consider a sphere or cylinder with radial symmetry then the boundary condition at 0 x = becomes the centre of the sphere or cylinder and is, with 0 λ = ∞ , 6 0 ( ) 0 x du x dx = = (13).

Quadrature
Eqn (3b) contains a quadrature over the generating function in the form In the case of slab geometry, , this expression becomes for a time independent where j w  is a weight function (Dehghan andSaadatmandi, 2008, Zhou andLi, 2017). For a uniform source of total strength T S , we may write S v S a n s cm This representation is adequate for a smoothly varying source and generating function, for example if the source is uniform, but for a localised source it would converge very slowly. In the case of a localised source which can be represented by a delta function in the form ( ) For a slab a number of sources can be represented across the slab, each located at a collocation point. However, due to the symmetry restriction, it is only possible to have a point or line source at the sphere or cylinder centre; if we wished to move it to a larger radius it would have to be in the form of a concentric annulus. Any asymmetry would require the angular variables to be included in 2 ∇ and that would increase the number of collocation points to unacceptable levels.
where j S is the source strength at each collocation point. We will show below how the safety factor Q depends on the source position and the geometry.
The modifications to include space dependence into the problem are embodied in a computer code and we will give below some numerical results and discuss their implications for low source startup. In order to use the Pál-Bell equation (1) for low source startup problems we refer the reader to I. There it is shown how to calculate the strength of a startup source that will reduce the probability of any rogue transient arising from the low density of neutrons, to a prescribed value. In short, it is necessary to reconstruct the probability distribution of the neutrons from the generating function and hence calculate the cumulative distribution where ( , ) S P n t is the inverse of the generating function ( ) is the fraction of the neutron population with number density less than n*. Later we will see that n* corresponds to the neutron density at the start of the deterministic regime. For safe startup it is essential that Q be very small. In order to perform the inverse, the saddlepoint method is employed as described in detail in I and briefly in Appendix B.

Numerical examples
It is important that the final results of the calculations employ six groups of delayed neutrons but because the running time of the associated computer code is so long with six groups (  hours), we illustrate some general points using one group of delayed neutrons. Later some comparisons with six groups will be given. Also it is useful to exemplify the collocation method on a relatively simple problem before using it for startup; for this we choose the survival probability.

The Survival Probability
The survival probability can be obtained directly from eqn (1) and written as ( , ) (0, , Physically, ( , ) W r t is the probability that a neutron injected at any point on a concentric surface of radius r at t=0, will survive to time t. To illustrate the collocation method, we will calculate the survival probability of a neutron in a reflected sphere with the reflector replaced by an effective boundary condition of the form Note that this boundary condition follows from the adjoint nature of the Pál-Bell equation.
For a radially symmetric sphere, ( , ) An interesting matter arises concerning symmetry. For example, if the initiating neutron is injected at a particular point at a given radius r, then the same subsequent behaviour will arise no matter where on that radial circle the neutron is injected. As long as the medium is homogeneous then spherical symmetry holds and we are justified in neglecting the angular variables in 2 ∇ . Thus, physically, ( , ) W r t is the probability that a neutron injected at point at r, at time zero, will survive to time t.
We have neglected delayed neutrons in eqn (28); this is not essential but simplifies the discussion (Williams and Pázsit, 2015). We also assume that the cross sections are independent of time and set t t s → − . To further simplify the equation we divide by In the above, ( ) 0,1 r ∈ and the initial condition is ( , 0) 1 W r = , i.e. no matter where the initial neutron is in the sphere, its survival at t=0 is certain. Only subsequently does the survival probability decrease.
If we convert eqn (28) to Gauss-Chevyshev-Lobatto form we find, using the transformation (1 ) / 2 where, for simplicity, we have kept only the quadratic term. The form of ( , ) B k j , which incorporates the boundary condition (26), is derived in Appendix A. This equation is solved for N=20 for a range of times and values of extrapolation distances s λ . Fig 1 shows the survival probability at times from 5 to 60 prompt neutron lifetimes and for a range of extrapolation distances from zero up to 50. An important observation is that, as time increases, so the initial flat distribution tends to a convex shape. 60 τ is essentially the asymptotic shape for t → ∞ . The effect of reflection is also interesting and shows marked changes near the boundary as would be expected, but the inside of the sphere (some 20 diffusion lengths from the surface) is virtually unaffected by the boundary condition, as illustrated by the example at 10τ . The data used are

The Startup Source
We now move on to the calculation of the safety factor Q and the magnitude of the associated startup source. The principle involved here depends, as noted above, on the value of Before entering into these calculations it will be helpful to look at the physical situation. We wish to ensure that the probability of a first persistent chain, i.e. the random generation of a very severe transient, will be reduced to a certain value Q. Alternatively, we must ensure that the fraction of neutrons less than n* is very small because it is a low neutron density that leads to fluctuations. First it is necessary to decide what is meant by a deterministic situation. Clearly this means the time at which the fluctuations have reduced to small relative values and we approach this situation as the neutron density or power increases. In practice, it may be shown by simulation that for a low neutron density the fluctuations are very large but, as the density increases, for example due to the insertion of reactivity, so the relative fluctuations which we define by ( ) 10 10 to − − would also suffice and we will have more to say about the practical aspects of this choice below. The important point is that the fluctuations are no longer changing statistically and we are in the so-called deterministic regime; actually it is not quite deterministic as fluctuations will always be present, but in this case they are very small. The value of the mean density at the start of the deterministic regime is important and by definition is is ( ) To be more precise we actually fix the value of Q, say These arguments were first developed by , although it is not easy to extract them from their discussion. As a clear example of how the increase in source strength dampens the fluctuations, we show below a simulation based on a Gamma distribution for the neutron density as a function of time as the reactivity increases in a ramp-like manner. The data used are as follows:  , has been discussed in detail in I which assumed a point model; now we extend that to spatial variation of the source.

The Saddlepoint Method
It will have been noted that eqn (1) is for the generating function of the probability distribution but that in practice we need eqn (30) which is the cumulative distribution.  employed the saddlepoint method to effect the inversion of the generating function to obtain Q directly. Their result is given in the form . These derivatives can be obtained from eqns (1), (2) and (3a) in a straightforward manner, full details are given in I but a summary is given below, viz: In addition, we will need equations for , , The final conditions on these quantities are These conditions apply for the neutron probability distribution. If the pdf of the precursors is required then the final conditions change to The boundary conditions corresponding to zero neutron and precursor density on the surface s r of the system are  which are given in Appendix B. We now have the complete set of equations that enable us to obtain the quantities used in the saddlepoint method, as described in eqns (B7)-(B9). The spatial aspects in the saddlepoint method use the collocation approach as described in Appendix A, see eqns A20-A23.
A further useful observation is that if one has the set of equations for the saddle-point method, i.e.,

Criticality
In order to introduce reactivity into the system we shall use a ramp model. First, however, consider the criticality problem of our system. Let us assume that we have a homogeneous sphere of radius R in which the diffusion equation may be written This may be rewritten as The criticality equation can be written For an arbitrary variation of reactivity with time, ( ) t ρ we may, from eqn (46), write the variation of the capture cross section as The initial value of the reactivity is defined by ( B is calculated from the root of eqn (45a) or (45b) according to geometry.

slab geometry
We illustrate numerically the above equations for one group of delayed neutrons and one speed neutrons for slab and spherical reactors with a series of localised sources and a uniform source. The neutron flux is zero on the boundary. The following data are used in Figs 4-6 for a ramp insertion of reactivity: Although one group of delayed neutrons is sufficient to show the general behaviour of the generating function with regard to system size, extrapolation distance, etc, it is not sufficient to get working values of Q. Thus we will repeat some of the calculations with six groups to show the error involved. The problem with using six groups is that it takes about 10 times longer to run the code which amounts to around 5 hours per run. We also note that due to the slow rate of change of reactivity and the maximum reactivity achieved in such calculations (about 1$), it is only necessary to use the quadratic approximation in the multiplicity terms.
We have checked that using seven terms gives no significant change in Q.
In order to compare the accuracy of the collocation method with that of our earlier method (I), namely the eigenfunction expansion, we have used both methods for a few examples. ). Three configurations are used; the source is at the core centre (a/2) and the source midway between centre and edge (a/4) and finally the source is uniformly distributed across the core. It is noted that the eigenfunction and collocation methods are close to each other especially in the region of interest ( ) We have also checked the result of placing two sources of half-strength in mirror image positions in the core. It appears that two half strength sources at a/4 and 3a/4 give the same value of Q as one full strength source at either a/4 or 3a/4. The execution times of the two methods are broadly comparable for the same accuracy. However, when six groups of delayed neutrons are included and higher multiplicities used, the collocation method proves more effective, especially when the geometry is of a realistic form. The symbol 2 nd N = − in Fig 4, where N is the number of collocation points. Fig 5 shows the variation of Q with reactor width, when the source is uniformly distributed across the core. The behaviour in fig 5 is interesting because it shows that Q tends to a limiting value as the width decreases. This suggests that a point model is developing and might explain why the point model is such a good approximation to assemblies such as Godiva. Conversely, as the width increases, Q again tends to a limiting value which, in the limit of infinite width, does tend to the classic point model as developed in I [see eqn (16)].
Thus the small width limit shows that there is a significant error in the point model for small systems due to the increased leakage. We also note that the results converge well with only a few collocation points; for example nd=5, 7 and 9 values of Q are within 0.01% of each other. When we move on to spherical geometry we will see some rather different behaviour nd=7 nd=5 x Figure 6: spatial variation of single particle density across core

Effective point model
It would be very convenient if we could employ the point model approach to describe the space dependence by simply modifying the strength of the source, i.e. by multiplying the source by a factor which is related to the importance function associated with the source position in the core. As an example of this, we show in are also close at the maturity time. We have not been able to deduce values of these weighting functions from first principles, but it seems likely that they can be obtained from an expression similar to that in the weighted source found from perturbation theory (Akcasu et al, 1971). This would seem to be a profitable line of investigation as it could simplify the solution of the Pál-Bell equation significantly.

Spherical and cylindrical geometry
We pass on now to more realistic geometries, namely the sphere and cylinder. Our first example uses the following data  j mat x mat j n r t g x t R v = = = y and also the corresponding relative standard deviation. It is important to note that the value of (0, ) mat n t decreases as s λ increases thereby reducing the effective source strength, and explaining mathematically why the source multiplication factor increases. The physical reason for this behaviour is more difficult to understand.
The shape of the density in Fig 11 is as we expect, and that of the relative standard deviation increases as we move towards the sphere surface. This is also to be expected because at the sphere surface the neutron density is lower than at the centre and hence involves greater fluctuations. There is also a general flattening of the distribution as s λ increases.  It is satisfying to note that despite the radius being 50 cm, five collocation points give a reasonably accurate representation.
In order to illustrate the relative differences between slab, cylindrical and spherical systems we show Fig 12. The source is in the form of a plane sheet at the centre of the slab, a line along the centre of the infinite cylinder and a central point for the sphere. A range of slab thicknesses and cylinder and sphere radii are used and it is clear that, for such localised sources, the size of the system is of minor importance as far as the value of Q is concerned in the range of radii 25-100 cm. The curves in Fig 12 are all for one group of delayed neutrons except for one case for the sphere which is for 6 groups. The one group case has a maturity time of 87.2 s. We give this example in order to show the significant, indeed vital, importance of using 6 groups. However, as we intimated earlier, the one group results are still useful for demonstrating general behaviour. We also show the insensitivity to the number of collocation points greater than five. A further result shown in Fig 12 is for the classic point model (bold line) which is significantly different from the spatial models and is strongly conservative. For example, if we require a safety factor of   ). We note that, in contrast to the central source location, the Q value for a uniformly distributed source is very sensitive to size. Physically, this is easy to understand because central sources are always in the position of maximum importance whereas for uniform sources each point has a gradually decreasing importance as the boundary is approached. A combination of localised and and uniform sources, as will be met in practice, indicates that system size will have an effect on Q. The traditional point model results lie in between the uniform and central values and there are significant differences between these values. The new source multiplication factor is around 2, which is physically acceptable because there is always some stochastic residual behaviour, but in this case it is so small as to indicate essentially deterministic behaviour. The point model result in Fig 14 shows, that for a centrally located source, the point model is highly conservative, whereas in the case of uniformly distributed source the point model is highly non-conservative. We also show (indicated by "space λ = ∞ ") the case of the spatial model with an infinite extrapolation distance. As can be seen, and as expected, this closely resembles the point model. It is tempting to approximate leakage by introducing a 2 DB term but this has little practical effect as its contribution to neutron losses in a typical reactor system is much less that that due to absorption. The main effect of spatial dependence is via the source position and distribution, as indicated by the graphical results above. , the slab has the smallest source multiplier of 32 followed by the cylinder at 37 and then slab at 71.

Summary and conclusions
The problem of dealing with the calculation of the magnitude of the neutron source in the startup of a nuclear reactor is of great practical importance. If the source is too weak then severe statistical fluctuations can arise, and if too strong unnecessary expense is incurred. The theory described above is based on the Pál-Bell equations which comprise methods of representing the neutron and precursor probability distribution in a compact form using generating functions. The main goal of the calculations is to obtain the source strength necessary to restrict the probability of any adverse stochastic effect to a prescribed value, say 8 10 − . We have described the way in which the Pál-Bell equation can be used and the methods of inverting the generating function to obtain the probability distribution function. From this one can ensure that the fraction of neutrons less than a prescribed value is very small thereby reducing the likelihood of unwanted fluctuations. To put it otherwise the probability of there being very few neutrons below some fiducial level is small. In order to carry out the calculations, we have relied on the general philosophy of earlier workers in the field, mainly Hurwitz and co-workers (1963), but with a significant variation in method. These early workers used the forward form of probability balance which, although applicable to a point model, is not suitable for the case of energy dependence or spatial variation of the source. The main advances in this paper are extensions to space dependence and to slab, cylindrical and spherical geometries. One group of delayed neutrons and six groups are used in the calculations. It is found that six groups are essential for practical calculations but that using one group is suitable for general scoping studies to investigate changes in geometry and system size. Using six groups does involve long computation times of several hours but this is possible, and results for practical situations are readily available.
Finally we comment on energy dependence. Calculations with the code CALLISTO, and using our earlier two group work (I), indicate that the effect of cross section energy dependence and the associated slowing down process, do not change the value of Q significantly provided the appropriate average thermal cross sections are employed. More will be written about this matter in future.

Appendix A Collocation methods in slabs, spheres and cylinders.
Let us assume that we have a one dimensional system such that ( ) u x , is the solution of the following second order differential equation.
( ) . The above equation then becomes ( ) 2 2 2 4 ( ) 2 ( ) ( ( 1)) ( ( 1)) ( ( 1)) ( ) ( ( 1)) 0; 1,1 2 2 2 2 We may now choose collocation points Eqn (A2) now becomes in an obvious notation The symbols 1 2 ( , ) and ( , ) d k j d k j are defined in Kopriva (2009). The advantage of the CGL method is that the boundary lies on collocation points and so may be separated off in order to apply boundary conditions. Thus we may re-write eqn (A4) as In the simplest case of zero value on the boundary 0 0 N u u = = and the additional terms are zero. Suppose now that we have more general boundary conditions of the form Transforming these conditions to the variable v , as described above, we find ( 1) and Using the CGL approximation for the first derivative leads to where / The terms 0 and N d d are defined as This operator does, of course, include the boundary conditions (A6). The solution of these linear equations gives the solution to eqn (A1) with boundary conditions (A6) to any degree of accuracy as N is increased. In general, for smoothly varying functions, relatively few collocation points ( ) 10 N < give adequate accuracy.
To extend the collocation method to two and three dimensions we refer the reader to Kopriva (2009) who shows that, for two dimensional x-y geometry, Poisson's equation 2 ( , ) ( , ) 0 x y S x y     , may be written in the form where the terms on the right hand side are the source and the known boundary conditions. An analogous expression arises for three dimensions.

Spherical and cylindrical geometry
The results above apply to slab geometry and to symmetric spherical and cylindrical geometry, i.e. no angular terms. For example, if we have a sphere or an infinite cylinder then the centre point is equivalent to the boundary condition (0) S d d N N d N d N d N N V d N U d N N T d d N N d N d N d For a cylinder we have the same value of ( , ) B k j except the numeral 8 in eqns (A15) The collocation form of the generating function is therefore defined by eqn (32) as: Note that 0 and N g g   are given by eqns (A18)  It will be noted that, in the above calculations, the reactivity via the change in capture cross section is inserted in a spatially uniform manner. In practice, the reactivity change would result from a locally placed absorber, i.e. the control rod or blade. This spatial location could be taken into account by regarding the capture cross section to be composed of two parts; a homogeneous part due to fuel and moderator and a highly localised part due to the rod. Thus in eqn A20 above we can write The localisation of the absorber has been discussed by Stacey (1969) and in some cases, for example in large reactors where spatial instability may occur, could have some effect. That is to say the neutron density on one side of the core could be lower than that on the other. This matter should be studied, although preliminary calculations indicate that it will not be a severe threat to stochastic stability.
are avoided unless one wishes to actually write down the equations. Even though this would be a useful exercise it would not be necessary to code them separately as their solution is readily incorporated into the code. Moreover, the only place where x and i y appear explicitly in the backward equations is in the final condition, e.g. ( , , ) 1 = − y  g x t t x , and so we may set x=1 and (1,1, ) 0 =  g t t . Thus we simply run the saddle-point code with new final conditions.

Appendix C: The Gamma distribution
In our previous work (I), we noted the general usefulness of the Gamma pdf which requires only the mean and variance of the distribution to completely specify it. As an example of this we have seen that the Gamma pdf leads to good agreement with the CALIBAN experiments (Williams, 2016). We also note that for a sufficiently small start-up rate it agrees well with the saddlepoint method for predicting values of Q. In addition, the Gamma pdf is a powerful tool for simulating stochastic behaviour by obtaining realisations of the neutron density to explain the physical processes behind the low source problem and the calculation of a safe source.
The conventional method devised by Hurwitz and co-workers to calculate the safe source in reactor startup is to find the maturity time tmat and hence to find the source multiplier and ( , ) P n t S is the neutron density pdf which we will take as the Gamma pdf, and which bypasses the need to solve the Pál-Bell equations. We wish to give some physical explanation to the argument which defines the maturity time as a suitable criterion for defining the onset of the deterministic region and for the calculation of the source multiplier M. In practice, we find a relationship between Q and ( ) / * mat M n t n  so that if S is the original source then M is the factor by which it must be multiplied to give the specified value of Q. For the case of S=426 n/s and a ramp reactivity of 0.017 $/s , in a reactor initially at a sub-criticality of -0.5$, we obtain the Hurwitz curve shown in Fig C1 which The Gamma pdf is known to be a physically acceptable equivalent of the actual neutron pdf and we use it in eqn C1 in the form, x y γ is the incomplete Gamma function. Although the Gamma pdf is physically reasonable, as we shall see it is not always so quantitatively.
In order to understand the actual physical behaviour of the pdf we may regard the Gamma pdf as a close analogy and use the data below to illustrate a number of points.  C2, we see the pdfs at various times after startup as the reactivity increases from its subcritical value at -0.5$ to a supercritical value, passing through criticality at 29.4s. We use five different times 10, 20, 30 , 40 and 50 s. The figure shows that, for the weak source, all pdfs diverge at n=0 and thereby emphasise the low density region. This behaviour is reflected in the large fluctuations in the density observed in Fig 3 of the text. The source S*, which is increased by the factor 43.8, yields pdfs that are zero at n=0 and have a maximum value which moves towards the deterministic mean as time increases. The behaviour of both curves can be explained if we examine the behaviour of ( ) t  as it appears in eqn (C2). Fig   C3 shows for both sources as a function of time, along with the variation of reactivity ( ) t  in $. It is clear from the form of the Gamma pdf that the value of ( ) t  is crucial. If it is less that unity, then the curve will diverge as 0 n  , on the other hand for ( ) t  greater than unity the pdf will go to zero at 0 n  . Numerically, we find that for the weaker source ( ) 1 t   at a time of 58 s. Thus the bulk of the pdfs will be divergent at n=0 until after 58 s and the maturity time of 76 s is reached. For the modified source S* we find that ( ) 1 t   when t= 0.48 s which means that only during a small period of time, very near the beginning of startup, is it liable to undergo divergence at n=0. The value of ( ) t  for S* goes up to nearly 200 and the associated pdf will have a maximum over most of its time period. Note that the relative standard deviation is / 1/ n σ η = . The maturity time mat t is unaffected by source strength.
We conclude by noting that the Gamma pdf is in principle useful for space dependent studies because all that is needed are the mean S N , i.e. a solution of the first derivative with respect . It is clear that the Gamma pdf underestimates the multiplication factor and is therefore nonconservative, also the error increases as the source weakens. In spite of that, the Gamma pdf remains a useful tool for guidance. Increasing the number of collocation points to 11, changes the results by an insignificant amount.

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.888962.