Multiplicity theory beyond the point model

Passive methods of nuclear safeguards determine the important parameters of an unknown sample from the statistics of the detection of the neutrons emitted from the item. These latter are due to spontaneous fissions and (α,n) reactions, enhanced by internal multiplication before leaking out. Based on the original work of Bohnel, the methodology of traditional multiplicity counting is based on the first three factorial moments of the number of neutrons, emitted from the sample due to one source event. These “Bohnel moments” were derived in the so-called “point model”, in which no space-dependence is accounted for, rather a uniform first collision probability is assumed for each neutron, irrespective of the position of its birth and its velocity direction, and, more important, it is assumed to be the same for all generations in the fission chain as for the source neutrons. The purpose of the present work is to derive the same factorial moments in a one-speed space-dependent model, in which the position and direction of the neutrons is accounted for, but (similarly to the original Bohnel model), no energy dependence is assumed. The integral equations for the moments are solved numerically with a collision number expansion. It is shown that compared to the space-dependent calculations, the unfolding method using the point model underestimates the fissile mass and the underestimation increases with increasing both of fissile mass and the value of α.


Introduction
The principles of determining the singles, doubles and triples neutron count rates used in identifying and quantifying unknown items of special nuclear materials are well known (Böhnel, 1985;Ensslin et al., 1998;Pázsit and Pál, 2008;Pázsit et al., 2009). The basic step of the derivation consists of the setting up a backward-type probability balance equation (master equation) for the number distribution of neutrons leaving the sample due to a single source event (spontaneous fission or an ða; nÞ event). It is assumed that the only reaction which the neutrons can undergo is induced fission, and that there exists a uniform first collision probability p which is the same for all neutrons, arising both from the source event and from induced fission, and irrespective of the position of the neutron. Since the effect of spatial or directional dependence of the neutrons on the results is not taken into account, this model is often referred to as the point model.
Because of using the backward equation approach (deriving a forward equation is not possible for this case), one has to go in two steps: first an equation for the number distribution of neutrons leaving the sample due to one single starting neutron is written down; then another master equation is derived which connects the source-event induced distribution with the single starting neutron induced distribution. Turning to the generating functions of these quantities, one can derive equations for the factorial moments of these distributions, out of which the first three of the source-induced event are of interest. In the point model, these are simple algebraic equations, which can be solved explicity in a recursive way (the equations for the higher order moment are linear, but contain powers of all the lower order moments).
In this paper we present a derivation of the equations for the same factorial moments, but instead of using the point model and postulating a universal first collision probability, the interactions of the neutrons will be treated in a space-angle dependent one-speed transport model. More precisely, the position and direction dependence of both the source neutrons and those propagating in the sample will be accounted for; however, the factorial moments of the number of neutrons leaving the item will be calculated irrespective of their exiting direction, i.e accounting for the detection process (the geometry and material composition of the detectors is not modelled). Extension of the model to include energy dependence is straightforward, similarly to the case of stochastic neutron transport in reactors (Pál, 1958;Bell, 1965). Likewise, the detector geometry can also be taken into account, although all these generalisations incur a substantial complication of the calculations.
The general theory, valid for all sample geometries will be given first. Then a quantitative solution will be given for a spherical item, since the spherical geometry simplifies the numerical work significantly. This restriction can be easily abandoned, at the price of requiring a larger effort for obtaining quantitative results. Work on extending the calculations to cylindrical geometry is under way and will be reported in a forthcoming publication. The quantitative solutions in this paper are obtained by a Neumann-series (collision number) iterative expansion scheme, applied directly (and separately) on the three moment equations.
With this technique, quantitative results are given for the first three factorial moments. The factorial moments for a given sample can be calculated both in the point model and in the spacedependent model, and compared to each other. In the safeguards context, the interesting question is how much deviation the application of the point model will incur as compared to the more realistic space-dependent model. This question can also be answered by converting the factorial moments to singles, doubles and triples rates, and performing the inversion procedure to extract the fissile mass from the measured multiplicity rates. Such a comparison was also performed. The results show that the difference in the fissile mass, as determined from the point model instead from the more accurate space-angle dependent model is relatively moderate, but it increases with both the sample mass and with the value of a (the relative contribution of ða; nÞ neutrons to the source neutrons). Although the bias is relatively moderate (up to 25%, depending on the item mass and the value of a), an important point is that the bias of the point model is not conservative, i.e. it underestimates the fissile mass. Based on the results, it is possible to suggest a correction procedure applied to the point model equations to extract the correct fissile mass.
One might object that in reality, there are many reasons other than the space-and angular dependence of the neutron transport inside the sample, which are not accounted for by the point model, such as energy dependence, reactions other than fission, accounting for the detection process etc., and that the effect of all those aspects can be readily simulated by Monte-Carlo calculations. Indeed, work to determine the bias of the point model with Monte-Carlo simulations, and hence to derive empirical correction factors has been performed in the past by introducing the so-called Weighted Point Model (Burward-Hoy et al., 2004). Nevertheless, this fact does not diminish the value of the type of investigations presented in this paper, whose goal is to provide both quantitative results, as well as insight. The need for getting insight into the origin of the bias of the point model was also expressed in the work of Croft et al. (2007) in which one of the contributing factors, namely the dependence of the leakage multiplication on the position of the source neutrons was investigated. Our goal here is to investigate the effect of fully accounting for the spatial and angular aspects of the neutron transport inside the item. This goal can be achieved in a fundamental way by the transport model applied here, not only quantifying the bias, but also lending insight into its origin in terms of neutron transport theory.
Another aspect is that by using the same stochastic formalism as the point model, but extended to spatial and angular effects, the point model and the space dependent model are based on equivalent premises, hence it is easier to ensure a consistent comparison between the two models. Last, but not least, the insight given by the analytical model, and the ease and speed of a deterministic calculation, makes it easier to suggest a correction procedure for the bias which is due to the neglection of space dependence in the point model.

Preliminaries
The general assumptions of the space-angle dependent model are as follows. The source of the neutrons in the item, consisting of fissile material, consists partly of spontaneous fission, leading to the emission of a random number of neutrons with a given number distribution p sf ðnÞ and intensity Q sf F, and partly of single neutrons generated in ða; nÞ reactions, with an intensity Q a . The only reaction the neutrons can undergo is fission (absorption and scattering can be included in a straightforward way). The macroscopic reaction cross section is equal to R f . The angular distribution of the neutrons from induced fission will also be assumed isotropic.
The source events occur randomly in space within the item with a uniform spatial probability density p R ðrÞ inside the volume V of the item, i.e.
The source neutrons will have independent, identical uniform angular distribution with a density To simplify the derivations, it is customary to introduce the number distribution p s ðnÞ and the intensity Q s of the (total) source events, which includes both the fission and the ða; nÞ neutrons in a weighted sum, as (Pázsit and Pál, 2008;Pázsit et al., 2009) For further simplification of the formalism, one introduces the factor a as the ratio of the rate of production of neutrons by ða; nÞ reactions and spontaneous fissions: The generating function q s ðzÞ of the number distribution p s ðnÞ of neutrons from a source event, and the factorial moments of the latter are defined as q s ðzÞ ¼ X 1 n¼0 z n p s ðnÞ; d n q s ðzÞ d z n j z¼1 m s;n These factorial moments are not known, because they include the unknown factor a (cf. (7)).
On the other hand, the number distribution, and hence the factorial moments of the spontaneous and induced fission, respectively, are known, since these are known nuclear physics quantities. These are denoted for the spontaneous fission as It is worth mentioning that in the safeguards literature, usually an alternative notation (subscript i) is used for the induced fission (Ensslin et al., 1998;Pázsit et al., 2009). However, here we follow the notations of Pázsit and Pál (2008) (subscript r), for the same reason as in the latter publication, i.e. to avoid confusion when i is used as a running index or in summation formulae.
In the derivations which will follow, the moments of the source distribution and that of the induced fission will appear as coefficients in the equations of the sought moments of the neutrons emitted from the item. In the end, for the unfolding of the fissile mass, the factorial moments m i ði ¼ 1 . . . 3Þ of the number of particles in the analytical formulae of the point model, these are rewritten in terms of the (known) moments of the spontaneous fission and the (unknown) parameter a.
3. Derivation of the master equations for arbitrary item shape 3.1. Single particle induced distribution The following integral master equation can be written down for the number distribution pðnj r; XÞ of the neutrons, leaving the sample due to one starting neutron inside the item with position r and direction X, by adding the probabilities of not having, respectively having a first collision somewhere on the path towards the surface of the item: pðn 1 j r þ sX; X 1 Þ pðn 2 j r þ sX; X 2 Þ . . . pðn k j r þ sX; X k Þ: Here, the quantity 'ðr; XÞ denotes the distance from the point r along direction X to the surface (boundary) of the item.
In one-speed transport theory, and especially in cases with some degree of symmetry, such as in slab and spherical geometry, it is convenient to express length scales in units of the collision mean free path (here being equal to R À1 f ), i.e. to swith to optical path length units. Since we consider here a homogeneous item, even in our case it simplifies the notations to switch to optical units. This means that we introduce the notations for the dimensionless quantitieŝ r ¼x;ŷ;ẑ f g ð12Þ and then re-denoter as r. With this, (11) will read as pðnj r; XÞ ¼ e À'ðr;XÞ d n;1 þ Z 'ðr;XÞ 0 ds e Às X 1 0 p r ðkÞ X pðn 1 j r þ sX; X 1 Þ pðn 2 j r þ sX; X 2 Þ . . . pðn k j r þ sX; X k Þ: where now 'ðr; XÞ; r and s are dimensionless quantities, the first of these denoting the distance to the boundary in optical path units. Introducing the generating function gðzj r; XÞ of the probability pðnj r; XÞ of Eq. (14) where q r ðzÞ is the generating function of the number distribution of the induced fission neutrons, defined in Eq. (10), with the quantity in the square brackets being its argument.

Source event induced distribution
In a similar way, for the number distribution PðN j SÞ of the neutrons leaving the item due to one source event S, and its generating function GðzÞ, the following equations are obtained: Hence the equation for the generating function GðzÞ reads as Eqs. (16) and (18) are the space-dependent analogues of the point-model based Böhnel equations,see Eqs. (11.38) and (11.39) in Pázsit and Pál (2008) for the generating functions hðzÞ (the equivalent of gðzÞ here) and HðzÞ (the analogue of our GðzÞ here). Of course our present equations are now more complicated due to the spatial and angular dependence of the neutron positions and directions, but the correspondence between the respective non-collided and collided terms in the equations for g and h is quite conspicuous. The conceptual similarity, as well as the difference between the point model and the space-dependent model will be returned to soon.
Eqs. (16) and (18) can be written in a more compact form by introducing the ''scalar" (angularly integrated) generating function gðzj rÞ as This quantity might appear as an analogue to the scalar (angularly integrated) angular flux (or rather its generating function), and we shall refer to its moments as ''scalar moments" in this paper. However, it has to be kept in mind that the direction X is that of the starting neutron and not the exiting one (in other words, nðr; XÞ is the analogue of the importance of the starting neutron), and the expression ''scalar" has to be understood in this sense. With the definition (19), Eqs. (16) and (18) can be written in the more compact form gðzj r; XÞ ¼ z e À'ðr;XÞ þ Z 'ðr;XÞ 0 ds e Às q r gðzj r0ðsÞÞ ½ ð 20Þ and where also the notation r0ðsÞ r þ sX ð22Þ was introduced. Eq. (20) can be integrated w.r.t. X to obtain an equation directly to the scalar generating function gðzj rÞ as is the generating function of the uncollided part. The above equations are valid for any item shape as long as only fission with isotropic distribution of the directions of the fission neutrons, and isotropic scattering (not included here) is accounted for. The only quantity which is dependent on the geometry is the distance to the boundary, 'ðr; XÞ. We will see a concrete expression for this quantity in a subsequent section when considering a spherical item.

First moments
The derivation of the equations for the factorial moments from the equations for the generating functions goes on the same lines as in the point model case. Denote the angular and scalar first moments as hnðr; XÞi ¼ @ @z gðzj r; XÞ j z¼1 nðr; XÞ; hnðrÞi ¼ @ @z gðzj rÞ j z¼1 nðrÞ; ð25Þ Then, for nðr; XÞ, one obtains from (20)  As is seen from Eq. (28), for the calculation of the expectation N of the source event induced emission, it is sufficient to know the scalar expectation nðrÞ of the single-particle induced emission number. This is of course a result of the isotropic angular distribution of the source and induced fission neutrons; accounting for anisotropic scattering would necessitate keeping the angular quantities. Keeping only scalar quantities will simplify the notations in the continuation. Integration of (27) yields an equation directly for the scalar expectation as being the uncollided contribution to the expectation of the singleparticle induced number of neutrons leaving the sample, for a starting neutron with isotropically distributed angular direction. It is a known quantity, and it will be the starting term of the collision number expansion form of the numerical solution.

Second factorial moments
From now on we will only keep the scalar forms of the singleparticle induced moments. Let us introduce the notations hnðrÞðnðrÞ À 1Þi mðrÞ ð 31Þ After a twofold derivation of the equations for g and G, respectively, one obtains ds e Às m r;2 n 2 ðr0ðsÞÞ þ m r;1 mðr0ðsÞÞ which can be re-written as

Equations for a spherical item
The calculations simplify significantly if one considers a spherical item. In such a case, due to symmetry considerations, it is easy to see that all appearing quantities will only depend on r ¼ jrj and l, where l is the cosine of the angle between r and X, i.e.

Imre Pázsit and Lénárd Pál
Annals of Nuclear Energy 154 (2021) 108119 To emphasise in the notations that length scales are expressed in units of the mean free path, we will use the variable x for the radial position in optical path length units, i.e.
instead of using the variable r of the original geometrical distance to the origin. The physical radius R of the spherical item is then given It is seen that in this model, unlike in the point model, there is no need to introduce a first collision probability p, since the collisions are described at the level of the individual neutrons due to the transport model applied. However, the transport model also defines a first collision probability for the isotropically distributed source neutrons, which can be calculated from the non-collided part. A numerical calculation of the first collision probability is readily possible for any geometry; however, for spherical geometry, a simple closed form analytic formula is available. This formula has long been known, and is given as (Bell and Glasstone, 1970;Williams, 1971) When making comparison with the results of the point model, Eq. (45) will be used to ensure that the calculations of the factorial moments refer to the same system, since the first collision probability of the source neutrons must be the same in both models. The difference between the more fundamental transport model and the point model is that in the latter, the first collision probability is assumed to be the same for all generation of neutrons, which is not true in reality. Although the source neutrons are distributed uniformly in the sphere, the spatial distribution of the neutrons emerging from the first, second etc. collisions is not uniform. The transport model used here accounts properly for this difference automatically, whereas the point model does not.

Equations for the generating functions
where q r is the generating function of the number distribution of the induced fission neutrons (see Eq. (8)), is the ''scalar" generating function, is the radial position of the neutron at a distance s away from x along direction l, and is the distance to the boundary of the sphere from the radius x along direction l (see Fig. 1).
Similarly to the general case, Eq. (23), an equation for the scalar generating function can be derived directly as Further, the relationship between the generating function GðzÞ of the probability PðN j SÞ of emitting N neutrons due to one source emission event S, and that of the single particle induced scalar distribution gðzj xÞ is given as In the above, account is taken for the fact that the singe event distributions emit neutrons isotropically, and the source event distribution is spatially homogeneous.
Eq. (53) shows that once the single neutron induced distribution (or its moments) are known, the source induced moments can be easily obtained by pure quadrature, without the need of solving any equations. Hence, only the calculation of the single particle induced moments constitutes a computational challenge.

Equations for the factorial moments
The equations for the angularly dependent moments can be obtained by taking derivatives of Eq. (47) w.r.t. z and those for the scalar moments from Eq. (51). For the first moment, we will write down both; for the rest of the moments only the equations for the scalar moments will be given.

Solution with a collision number expansion
It is seen from the foregoing that only the single neutron induced scalar moments obey an integral equation which has to be solved, the source event induced moments (which are the main interest of the calculations) can be obtained by quadrature from the single-particle induced scalar moments. In the simplified case of spherical geometry of the item, obeying a high order symmetry, the solution of the moments due to a single initiating particle can be obtained relatively easily by various deterministic numerical techniques, similarly to other cases of space-and energy dependent stochastic transport problems, such as the S N or finite difference methods (Endo et al., 2008;Saxby, 2018).
Here we suggest the use of a collision number (Neumann series) expansion method (Case et al., 1953), which was used earlier by one of the present authors for solving the non-linear integral equations of invariant embedding for electron and positron backscattering from surfaces (Glazov and Pázsit, 2004;Glazov and Pázsit, 2007). This is a straightforward and easy-to-apply procedure, as will be shown below. As it is mentioned in Case et al. (1953), this method works well only for subcritical cases with short fission chains, where the collision number expansion converges fast. This is actually the case in nuclear safeguards where the fissile items are far away from criticality.

First moment
For the scalar mean nðxÞ, the procedure of solving Eq. (55) goes as follows. The starting (zeroth) term consists of the expectation of the scalar non-collided part, Eq. (56) The scalar once-collided part n 1 ðxÞ is then given as The full solution for the scalar ''flux" is obtained by summing up to all collision numbers, i.e.
nðxÞ ¼ X 1 k¼0 n k ðxÞ ¼ n 0 ðxÞ þ X 1 k¼1 n k ðxÞ n uncoll ðxÞ þ n coll ðxÞ: Here, the first term is the uncollided part, and the rest is the collided part. An approximate solution with desired accuracy can be obtained by a proper choice of the cut-off of the iteration at k ¼ k max .
At this point it is worth pointing out an analogy with the equations and the corresponding results of the point model. In the point model equations, Eq. (11.41) in Pázsit and Pál (2008), the expectation h 1 of the neutrons leaving the item induced by one starting neutron (the point model analogue of our nðxÞ) is given as where p is the (uniform) first collision probability, and 1 À p is the non-collision probability. Since the solution (72) is only valid when p m r;1 < 1 (i.e. when the system is subcritical), the denominator can be expanded into a Taylor series, yielding h 1 ¼ ð1 À pÞ½1 þ p m r;1 þ ðp m r;1 Þ 2 þ ðp m r;1 Þ 3 þ . . .

ð 73Þ
In this form the similarity with Eqs. (70) and (71) is obvious, since the r.h.s. of (73) also contains the exactly none, once, twice etc. collided neutrons, multiplied with the leakage probability collision and non-collision probabilities are dependent on the position and the direction of the neutron, and, this way, they are different for different generations (e.g. for neutrons from the source and from induced fission).
The similarity actually can be extended even to the solution method. Namely, the equation for h 1 , given as which of course can be solved directly leading to the compact form (72), can also be solved by the same iteration technique. Noticing that p m r;1 < 1, on can write that which will lead to the same solution as (73), which is equivalent to (72).

Solution for the second and third moments
The equations for the higher order moments can be solved with the same techniques sequentially, if the solution of the lower order moments is already available. Hence, the solution of Eq. (60) Again, the full solution is given by where the summation can be cut after a suitable value of k max when the series expansion has converged with the desired accuracy. The expansion for the third moment wðxÞ goes in a completely similar manner, assuming that the inhomogeneous term BðxÞ of Eq. (66) is already known from the solutions of the first two moment equations for nðxÞ and mðxÞ.
This way, the solution of all three moments can be obtained solving them in sequence, starting with the first order moments and progressing upwards. Actually, the method is readily applicable to any order moments, at least in principle, since the analytical form of the inhomogeneous part of the model equations, similarly to AðxÞ and BðxÞ, can be obtained from repeated differentiation of Eq. (51). The homogeneous part of all moment equations is the same, which is a general property of the backward master equation. In pure time-dependent problems, this property has the consequence that the solution of the first moment equation serves as the Green's function of the higher order moments (Pázsit and Pál, 2008). The Green's function technique can also be applied in the present problem, but it is more involved than in the pure timedependent cases, and the calculation and application of the Green's function for the present problem will be described in a separate publication.

Quantitative results and comparison with the point model equations
For the numerical work, one needs the material constants (the factorial moments of the number of neutrons emitted in spontaneous and induced fission) and the size of the item. The material constants were taken from Enqvist et al. (2006), corresponding to spontaneous and induced fission in a sample of 20 wt% 240 Pu and 80 wt% 239 Pu. In Enqvist et al. (2006) the full distributions p sf and p r were determined from Monte-Carlo simulations of the above item compositions, from which then the factorial moments could be determined. The first three factorial moments of the spontaneous and induced fission neutron distributions for such an item are given in Table 1 below. With the given material mix, the distribution of the spontaneous fission neutrons should be essentially determined by that of 240 Pu, since the spontaneous branching ratio of 239 Pu is rather small. Indeed, a comparison with the consensus values of both the number distribution and the factorial moments of spontaneous fission of 240 Pu given by Santi and Miller (2008) shows an excellent agreement. Since in the numerical work also calculations will be made for a case with a ¼ 0:5, the factorial moments m s;i ; i ¼ 1 . . . 3 of the number distribution p s ðnÞ of the source event, from Eqs. (7) and (8), are also given in the Table. Here it was assumed for simplicity that the presence of oxides or other light atomic number materials, leading to the emission of ða; nÞ neutrons, will not alter the number distribution of the spontaneous and induced fission neutrons.
Regarding the size of the sphere, calculations were made with several radii X, starting with very small item sizes where the first collision probability is close to zero (and hence the results from the point model and space dependent model are identical) to large items, which approach criticality; where one can expect larger differences between the two models. These material and geometrical data are sufficient for the calculation of the factorial moments in the space dependent case.
Calculation of the factorial moments in the point model requires the specification of the first collision probability p. As mentioned earlier, in order to have a correct comparison between the space dependent and point model values, the first collision probability in the point model should correspond to that given by the transport model of the space-dependent calculations. As was also remarked, for the homogeneous sphere with isotropic source and fission neutron directions, an analytical formula is available (Bell and Glasstone, 1970;Williams, 1971) as given by Eq. (45). From this formula, the first collision probability for a sphere with a given radius can be calculated numerically, and vice versa.
In the calculations of the moments of the space-dependent model with the collision number expansion method, we used the numerical options of Mathematica, version 12.1.1.0 (Wolfram Research, 2020). For the numerical work we used the high performance computing cluster Hebbe of Chalmers University of Technology. Using Mathematica has the advantage that extensions to other geometries, even irregular item shapes, are straightforward by using the analytical geometric functionality of Mathematica.

Numerical results
In order to be able to compare the influence of the ða; nÞ neutrons on the results, calculations were made for two scenarios for each sphere size: one with a ¼ 0, which corresponds to a metallic Pu item, and another one with a ¼ 0:5, corresponding to an oxide Pu item (or one with a matrix material of light nuclei).
The dependence of the mean values (first factorial moments) of the emitted particles, as functions of the radius of the sphere in optical units, is shown in Fig. 2  It is also seen that with the increase of the sphere size, the first moment of the space dependent model increases faster than that of the point model, for both values of a. This indicates that the space dependent model accounts for a larger internal multiplication than the point model, which must be due to the fact that the second and further collision probabilities must be higher than the first collision probability. However, the differences between the two models become significant only for very large items, close to criticality, which is not relevant for safeguards applications. The maximum item diameter in the figure, X ¼ 0:45, corresponds to a first collision probability p ¼ 0:265, which, as will be shown below, is close to the critical value. This value corresponds to a value of the leakage multiplication M over 4 in the point model, which is beyond the practically interesting cases in safeguards applications.
The expectation of the number of emitted neutrons is larger for the case of a ¼ 0 than for a ¼ 0:5 for any sphere size, which is understandable, since the source multiplicity is lower for a ¼ 0:5 (see also Table 1). Quantitative values for the doubles (second factorial moments) are shown in Fig. 3. Again, for small item sizes, they are both equal to the second factorial moment of the source for both a values, and with increasing item size they both increase. Also, the results obtained from the space dependent model are higher than those from the point model. It is also seen that the differences between the point model and the space-dependent model become larger and start already at smaller sample sizes than for the first moment (note the different range on the X-axis compared to Fig. 2). Similarly for the case of the first moment, the second moment values are larger for the case of a ¼ 0 at any sphere size.
The trend for the third factorial moments continues, as is shown in Fig. 4. Here the dependence of the factorial moments is shown only up to X ¼ 0:3 (corresponding to a first collision probability p about 0.2). It is seen that the deviation between the point model and the space dependent model starts at even smaller item sizes (first collision probabilities) than for the first two moments. At X ¼ 0:3 the result of the space-dependent model is nearly 50% larger than that of the point model.
The fact that the space dependent model predicts larger moments than the point model can be understood by the difference in the critical sizes for the given material properties that the two different models yield. In the point model, from the leakage multiplication M (not to be mixed up with the second factorial moment M of the number of emitted neutrons due to one source event in the space dependent model), given as (Pázsit and Pál, 2008): it is seen that the first collision probability corresponding to a critical system is given as With the numerical values used, this gives a critical p as 0.318, and a corresponding critical size On the other hand, an exact one-speed transport theory calculation yields the critical size of the sphere as (O'Rourke, 2019) This indicates that the space-dependent model accounts for a higher internal multiplication than the point model, since a smaller sphere radius is needed to make the system critical. To put it another way, the collision probability is higher for the later generations than the first collision probability p.

The bias of the point model
So far we have only compared the factorial moments of the number of neutrons leaving the item due to one source event. The interesting question is what effect the difference in these factorial moments has on the accuracy of the determination of the parameters of the item, primarily the fissile mass. The fact that the moments of the space dependent model are systematically higher than those of the point model does not immediately tell what error, or bias, the application of the point model yields when it does not predict the moments accurately.
To find out the answer to this question, one has to convert the above factorial moments into singles, doubles and triples count rates (S; D and T rates) (Ensslin et al., 1998), since only these multiplicity rates are measurable and, unlike the factorial moments, they carry information on the fissile mass through the sample fission rate F. The method how the fissile mass is recovered from the point model is also based on the multiplicity rates, as will be shown shortly.
The relationship between the factorial moments and the multiplicity rates is rather simple, and it is the same for both models. For the space dependent model, they are given as where N; M and W are calculated by the space-dependent model. In the above, F is the fission rate which is directly proportional to the fissile mass, m sf ;1 is the first factorial moment (expectation) of the spontaneous fission neutrons, and F ð1 þ a m sf ;1 Þis the total source intensity (see also Eq. (6)). Further, e is the detector efficiency, f d and f t are the so-called doubles and triples gate factors (Ensslin et al., 1998). Since these latter affect the moments of both the point model and the space-dependent model exactly the same way, their actual values will not affect the bias of the point model, and hence they can be assumed to be unity in the quantitative work. Since the bias of the point model will be defined as the ratio of the fission rate obtained by the application of the point model on the ''true" measured data (represented here by the multiplicity rates of the space-dependent model) to the true fission rate (input parameter to the calculations), the fission rate can also be taken as unity, and the bias of the point model will be simply given by the result of the point model by using the multiplicities of the space dependent model.
To determine the bias of the point model we hence recite briefly the procedure of how it is used to determine of the fission rate. As is known (Böhnel, 1985;Pázsit et al., 2009) where the quantity M is the leakage multiplication defined in Eq. (82).
Based on the point model formulae above, the unfolding of the fission rate from the measured S; D and T rates goes in two steps (Ensslin et al., 1998). First, from (89)   Here the multiplicity rates S; D and T are now denoted by boldface letters, because they have a key role in the determination of the bias of the point model. The leakage multiplication M is then obtained as the real root of Eq. (92). In possession of M, the fission rate is obtained from (89) and (90) as (Ensslin et al., 1998) According to the foregoing, the bias of the point model can be obtained as follows. The multiplicities S sp ; D sp and T sp of the space dependent model are calculated from (86)-(88) with F F space ¼ 1. These are then substituted for S; D and T in (93)-(95), and Eq. (92) is solved for M M space . This value is then used, together with S sp and D sp in the right hand side of (96), and the F so obtained is equal to the bias C.
The bias C of the point model was calculated according to this procedure for both a ¼ 0 and a ¼ 0:5. The dependence of the bias on the sphere size X is shown in Fig. 5. The figure shows that the bias of the point model increases with both the item size and the value of a. With a ¼ 0, the bias remains quite moderate even for X ¼ 0:3 (which corresponds approximately to p ¼ 0:2). However, for a ¼ 0:5, the bias is nearly 25% at X ¼ 0:3 (p ¼ 0:2). Although the bias is relatively moderate, it is non-conservative, i.e. application of the point model will underestimate the fissile mass. These results are consistent with the earlier observed bias of the point model, and in particular with the results of Burward-Hoy et al. (2004), who also note the non-conservative bias of the point model, i.e. the underestimation of the fissile mass. In this latter work the empirically calculated bias values are given as functions of the leakage multiplication. Since the leakage multiplication is different for the point model and for the space-dependent model, we have instead chosen the sample size (uniquely related to the sample mass), which makes it more straightforward to define a correction factor for the determination of the sample mass via the (biased) point model results. If the value of a is known, then for the fission rate predicted by the point model, one takes the corresponding X value and the curve corresponding to the a value in Fig. 5, to find the correction factor to be applied to get the correct fission rate.
In reality the value of a is not known, rather it has also to be determined by the same (biased) procedure by the point model, hence the determination of the correction factor requires a multivariate analysis. An alternative approach could be that, instead of determining correction formulae to the point model, to develop machine learning methods to extract the fission rate and the a factor directly from the measured multiplicities as input. The machine learning based unfolding method, such as an artificial neural network, can be trained by a large data set produced by the space dependent model for various item sizes, compositions and a values, similar to methods developed to unfold item data from combined neutron-gamma measurements (Enqvist et al., 2010). The investigation of the bias of the point model in the determination of a, and the feasibility of determining bias factors by multivariate analysis or developing machine learning techniques by training sets provided by the space-dependent calculations will be deferred to later work.
It might be interesting to understand the reason why the point model underestimates the fissile mass. This is the easiest done by considering the case when a ¼ 0. In that case the point model multiplicities (89)-(91) simplify, and the F value given by the point model can be simply expressed from (89) From this it immediately follows that the bias of the point model can simply be expressed as In the above, the ratio between the singles for the space dependent and point model, respectively, in the numerator in the last equality is the same as the ratio of the first moments of the space dependent model and the point model, respectively. This has already been calculated, and the results are shown on the left side of Fig. 2. It shows that the space dependent S values are increasingly larger than the point model ones, hence the ratio S space =S point is larger than unity and increases with the size of the item. However, a simple quantitative analysis shows that the same is valid for the ratio of the leakage multiplicities M space =M point too, except that the ratio between these two latter is much larger. This is shown on the right hand side diagram of Fig. 6. For comparison, the left hand side shows the dependence of the singles rates for the same range of the item sizes X. It is seen that the leakage multiplication of the space dependent model (which is in the denominator of (99)) exceeds that of the point model much more than the singles rates of the space dependent model that of the point model, and this leads to the underestimation of the fissile mass by the point model.
The reason that the deviation between the leakage multiplications is larger than that of the singles (first moments) between the two models is that the leakage multiplication depends on all three multiplicities (factorial moments). This is seen from the way it is determined through (92)-(95). Whereas the first moments (singles) differ only slightly from each other, the second and third moments differ more significantly, as indicated by Figs. 2-4. Hence the larger deviation between the leakage multiplications than between the first moments.
Although a similar simple transparent interpretation cannot be made for the bias of the point model for the case when a > 0, an intuitive interpretation still can be made. It is a general tendency in multiplicity counting that the higher moments are both larger than the lower moments, and also that they increase faster than the lower order moments with the sample mass. The moments of the number of particles leaving the sample reflect both the multiplicities of the source and the multiplicities of the internal multiplication in the sample. The multiplicity of the source neutrons is not related to the fissile mass, it is only the internal multiplication which is related to it. An increasing a means a decreasing source

Conclusions
Equations were derived for the probability distribution of the number of neutrons, leaving a sample of fissile material with spontaneous fission as the neutron source, in a space and angle dependent one-speed transport model. From these the first three factorial moments, which are used to calculate the singles, doubles and triples rates of multiplicity counting of nuclear safeguards, were calculated by a collision number type expansion. The results show that for small item sizes (small first collision probabilities) the point model and the space dependent model give very similar results, as expected. With the increase of the item size, the deviation between the point model and the space dependent model increases, and the deviation becomes larger for the higher order moments. A quantitative evaluation of the error in estimating the sample size with assuming the point model was made based on these results, which showed that compared to the case when one accounts for the spatial aspects of the internal transport of neutrons inside the item before leaking out, the point model underestimates the fissile mass of the item. The tendencies found, as well as the quantitative values, are in agreement with earlier experience in the field Burward-Hoy et al. (2004).
The formalism developed in this paper can readily be extended to more complicated geometries. In particular, calculations for a cylindrical geometry only require the introduction of two more parameters, which will still lead to manageable running times in the numerical work. Work is already continued into this direction. Extension to the calculation of higher order factorial moments is also simple. Further plans include the investigation of the use of the Green's function technique for the calculation of the factorial moments, as well as elaboration of multivariate unfolding techniques, probably based on machine learning methods, to elaborate correction factors to the point model equations, or to extract the correct fissile mass from the measured multiplicities directly by machine learning methods. Taking also energy dependence into account would be substantially more complicated, and for such cases numerical methods other than the collision number expansion might appear more feasible.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.