# Two-Group Theory of the Feynman-Alpha Method for Reactivity Measurement in ADS

Authors: Lénárd Pál; Imre Pázsit
Journal: Science and Technology of Nuclear Installations
 (2012)

---

## Abstract

The theory of the Feynman-alpha method, which is used to determine the subcritical reactivity of systems driven by an external source such as an ADS, is extended to two energy groups with the inclusion of delayed neutrons. This paper presents a full derivation of the variance to mean formula with the inclusion of two energy groups and delayed neutrons. The results are illustrated quantitatively and discussed in physical terms.


---

## Body

## 1. Introduction

Methods of online measurement of subcritical reactivity, in connection with ADS, have been studied over a decade by now. Both deterministic methods, such as the area ratio or Sjöstrand method [1] (pulsed measurements), and stochastic or fluctuation-based methods (Feynman- and Rossi-alpha methods) have been investigated. What regards the latter class of methods, the theory of classical systems, based on a stationary source with Poisson statistics, had to be extended to the case of nonstationary (pulsed) source with compound Poisson statistics (spallation source, generating several neutrons simultaneously in one source emission event). Regarding the pulsed sources, both narrow (instantaneous) as well as finite width pulses with various pulse shapes were considered, both with “deterministic’’ (synchronised between source emission and counting interval) and stochastic (nonsynchronised) pulse injection. An overview of the field can be found in [2]. The most general treatment of all the above mentioned cases is found in Degweker and Rana [3] and Rana and Degweker [4, 5].In this paper, we will discuss another aspect of stochastic reactivity measurement methods, which is related more to the system properties than those of the source. The new aspect is to take into account the energy dependence of the neutrons by the use of a two-group approach. All work so far in which compact analytical results could be obtained in this area (calculation of the Feynman- and Rossi-alpha formulae) was made by the use of one-group theory. This is justified by the fact that the methods were used in thermal systems, where the neutron population and hence its dynamics is dominated by thermal neutrons. However, many of the planned ADS concepts will use a core with a fast spectrum, in which the dominance of thermal neutrons will be significantly reduced. In terms of a two-group approach, unlike in a thermal system where there is one time or decay constant in a pulsed experiment, there will be two components in the temporal response with two different time constants and with comparable amplitudes. One indication for this possibility comes from the area of nuclear safeguards, where such effects have already been investigated, as it will be described below.Actually there have been experiments, such as in the EU-supported project MUSE [6] and the Yalina experiment [7–9], where the fitting of Feynman- and Rossi-alpha measurements required more than one exponentials. Although the appearance of more than one decay constant may have also other reasons (e.g., the presence of a reflector, i.e., a multiregion system), the energy aspect is clearly one possibility to lead to the occurrence of two different decay constants. In the theoretical work so far on the explanation of the multiple exponentials or multiple alpha modes, the emphasis was on the spatial effects as being the reason for the appearance of the multiple alpha modes [10–12]. In [10, 11] a general energy dependent framework was used which, in order to arrive to explicit expressions, requires the possession of the fundamental and higher order alpha-eigenfunctions of the space-energy dependent transport equation. However, due to the continuous energy treatment, no discrete modes can be attributed to the energy dependence.The purpose of the present paper is the elaboration of the two-group theory of the Feynman-alpha method, based on the backward master equation approach. Such an energy-dependent extension of the existing theory might have a relevance also to on-going and future ADS experiments, such as the European FP7 project FREYA. This work has actually been started already by the present authors, although in a different setting. In nuclear safeguards, identification of fissile material can be achieved by detecting the temporal decay of fast neutrons from an unknown sample, following irradiation by a pulse of fast neutrons. Appearance of a second decay constant is an indication of the presence of fissile material. This is the so-called differential die-away analysis (DDAA) method. A stochastic generalisation of this method, to the use of a stationary (intrinsic) source and the measurement of the time correlations, the so-called DDSI (differential die-away self interrogation) method, was recently suggested by Menlove et al. [13]. The theory of the DDSI method, by way of the extension of the Rossi-alpha method to two energy groups, without delayed neutrons, was recently given by the present authors based on both the backward [14] and the forward master equation approach [15].In the present paper, we extend the treatment by the inclusion of delayed neutrons in the formalism. Understandably, the calculations get quite involved. Although a full analytical treatment is possible, many expressions in the final results become too extensive to be quoted explicitly. Hence these will not be given and analysed in this paper. Likewise, the question of how to extract the subcritical reactivity from the measurement of the two time constants will not be discussed here, rather it will be given in a subsequent communication. Instead, here we focus on the formulation of the problem and the full derivation sequence until arriving to the final result. The intermediate and final results will be discussed in physical terms.


## 2. General Principles

As usual in the context, we shall assume a homogeneous infinite medium with properties constant in space. We shall use a two-group theory model, that is, describe the neutron population with two type of neutrons: fast and thermal. One group of delayed neutron precursors will be assumed. The possible neutron reactions are absorption of both the fast and thermal neutrons, downscattering (“removal’’) of neutrons from the fast group to the thermal, and thermal fission, which produces a random number of fast neutrons according to a probability distribution function. At such a thermal fission also at most one delayed neutron precursor can be generated with a certain probability. The decay of this precursor will lead to the appearance of one fast neutron. For simplicity, fast fission will be neglected. It can be easily incorporated into the model, at the expense of some further complication of the calculations, but without any essential problem. Figure1 illustrates the possible reactions induced by the fast and thermal neutrons and the delayed neutron precursors, respectively.Figure 1
The possible reactions which the three different particles can undergo. Fast neutrons (red circle) can be detected (yellow), absorbed (white), and thermalized (green). Thermal neutrons (green) can be absorbed, and through fission lead to a random number of delayed neutron precursors (grey) and fast neutrons (red). Delayed neutron precursors decay into a fast neutron.A word on the notations used is in order here. In a two-group model of reactor physics, the indices 1 and 2 are used to denote the fast and the thermal group, respectively. This notation will be used also in this paper, whenever it will not lead to confusion. For practical reasons, the delayed neutron precursors will be taken as group 3. However, often it will be simpler and more practical to refer to the fast and thermal neutrons and the precursors with the notationsF, T, and C, as it is seen also in Figure 1. In traditional reactor physics texts, short-hand notations are used for denoting the first factorial moments of the neutron population and that of the fission neutron distribution, such as ν and the Diven factor Dν for the latter. In the literature on neutron fluctuations, the transition probabilities (more correctly, transition intensities) are usually denoted by λi, where the subscript i stands for the type of reaction (capture, absorption, removal, and fission) (see, e.g., [2, Part II]). These transition intensities are related to the corresponding macroscopic cross sections of the particular reactions and the neutron speed.However, in this paper we will keep a more general mathematical system of notations for the factorial moments and transition intensities. This is partly because the purpose of this paper is to describe the formalism used in a general setting, and partly in order to relate the work reported here to our preceding paper, where the two-group generalisation of the Feynman- and Rossi-alpha methods was introduced [14]. Accordingly, the transition probabilities will be denoted in a way similar to that in [14], that is, for the fast neutrons one has(1)QF=QFA+QFT,
where QFA and QFT are the intensities of the absorption and thermalization of fast neutrons, respectively. Similar, self-obvious notations are used for the thermal neutrons. In the numerical work, there will be no attempt to relate these reaction intensities to cross-sections of a real reactor in this paper; such a coupling to realistic systems will be made in a subsequent work.The final quantity we need in order to formulate a probability balance equation is the number distribution of fast neutrons and delayed neutron precursors in a thermal fission event. Denote byf(k,ℓ) the probability that a thermal neutron produces k≥0 fast neutrons and ℓ≥0 delayed neutron precursors, that is, particles of type C. Suppose that(2)f(k,l)=fp(k)fd(l),
that is, the numbers fast neutrons and delayed neutron precursors created in one reaction are independent. Further, let fd(ℓ) be given by
(3)fd(l)={1-q1(d),ifl=0,q1(d),ifl=1,0,ifl>1,
where q1(d)≤1. The probability that a delayed neutron precursor produced at time t=0 decays to a fast neutron during the time interval not larger than t≥0 is given by(4)Td=1-e-λt,λ≥0.
With this all quantities that are needed to formulate the problem are defined.


## 3. Description of the Basic Process with One Starting Particle

Since we are going to use the backward master equation approach, first we will need the neutron and precursor distributions, generated by one single starting particle (a fast or thermal neutron or a delayed neutron precursor). First the number distribution of generated particles will be studied and explicit results derived for the first two moments. In the next section the detection process will also be accounted for, and in the last section the variance and the mean of the number of neutrons detected in a time interval, induced by a stationary external source of fast neutrons, will be calculated from the results of the preceding two sections.Let us introduce the random functionsnF(t),nT(t), and nC(t), giving the numbers of particles of types F, T, and C, respectively, at the time moment t≥0. Define the probabilities:(5)P{nF(t)=n1,nT(t)=n2,nC(t)=n3∣Sj}=p(n1,n2,n3,t∣Sj),
where the conditions Sj, j={1,2,3},{Sj}={F,T,C}, indicating the type of particle starting the process, are defined as
(6)S1≡F={nF(0)=1,nT(0)=0,nC(0)=0},S2≡T={nF(0)=0,nT(0)=1,nC(0)=0},S3≡C={nF(0)=0,nT(0)=0,nC(0)=1},
respectively. Introduce also the corresponding generating functions:
(7)g(z1,z2,z3,t∣Sj)=∑n1=0∞∑n2=0∞∑n3=0∞p(n1,n2,n3,t∣Sj)z1n1z2n2z3n3.
With these notations, based on the probabilities of the mutually exclusive events of the starting particle not having or having a reaction within the time interval (0,t), and the summing up of the probabilities of the mutually exclusive events generated by the respective reactions, one can write that
(8)g(z1,z2,z3,t∣F)=e-QFtz1+QFA∫0te-QF(t-t′)dt′+QFT∫0te-QF(t-t′)g(z1,z2,z3,t′∣T)dt′,(9)g(z1,z2,z3,t∣T)=e-QTtz2+QTA∫0te-QT(t-t′)dt′+QTF∫0te-QT(t-t′)qp[g(z1,z2,z3,t′∣F)]×qd[g(z1,z2,z3,t′∣C)]dt′,
as well as(10)g(z1,z2,z3,t∣C)=e-λtz3+λ∫0te-λ(t-t′)g(z1,z2,z3,t′∣F)dt′,
where(11)qp(z)=∑k=0∞fp(k)zk,qd(z)=∑l=0∞fd(l)=1-q1(d)+q1(d)z.
For later use, we note that(12)[djqp(z)dzj]z=1=qj(p),[djqd(z)dzj]z=1=qj(d)=δj1q1(d)
are the factorial moments of the number of prompt and delayed neutrons in a fission event. The relationship with the traditional notations is given by(13)q1(p)=νp=(1-β)ν,q1(d)=νd=βν,
with ν=νp+νd, and β=νd/(νp+νd).


### 3.1. Expectations of the Numbers of Particles of Different Types

By using (8)–(10), one can derive equations for the expectations of the numbers of fast and thermal neutrons and the delayed neutron precursors, that is, particles of types F, T, and C, respectively. With obvious notations, for the expectation (first moment) m1(F)(t∣Sj) of the number of fast neutrons, induced by one starting fast neutron, thermal neutron and delayed neutron precursor, respectively, one obtains the equations:
(14)m1(F)(t∣F)=[∂g(z1,z2,z3,t∣F)∂z1]z1=z2=z3=1=e-QFt+QFT∫0te-QF(t-t′)m1(F)(t′∣T)dt′,m1(F)(t∣T)=[∂g(z1,z2,z3,t∣T)∂z1]z1=z2=z3=1=QTF∫0te-QT(t-t′)×[q1(p)m1(F)(t′∣F)+q1(d)m1(F)(t′∣C)]dt′,m1(F)(t∣C)=[∂g(z1,z2,z3,t∣C)∂z1]z1=z2=z3=1=λ∫0te-λ(t-t′)m1(F)(t′∣F)dt′.
Equations for the quantities m1(T)(t∣Sj) and m1(C)(t∣Sj) can be derived in a completely similar manner; these will not be given here, for brevity.The arising integral equation system can be readily solved by Laplace transform methods. Introduce the Laplace transforms:(15)m̃1(Sj)(s∣Si)=∫0∞e-stm1(Sj)(t∣Si)dt,
then, from (14) one obtains
(16)m̃1(F)(s∣F)=1s+QF+QFTs+QFm̃1(F)(s∣T),m̃1(F)(s∣T)=QTFs+QT[q1(p)m̃1(F)(s∣F)+q1(d)m̃1(F)(s∣C)],m̃1(F)(s∣C)=λs+λm̃1(F)(s∣F).
After elementary algebra one obtains(17)m̃1(F)(s∣F)=(s+QT)(s+λ)N(s),
where(18)N(s)=(s+QF)(s+QT)(s+λ)-QFTQTF[q1(p)(s+λ)+q1(d)λ]
is a third-order polynomial of s. It can be proven that the roots of the equation:
(19)N(s)=s3+(QF+QT+λ)s2+[QFQT+λ(QF+QT)-q1(p)QFTQTF]s+λ[QFQT-QFTQTF(q1(p)+q1(d))]=0
are all real. Hence, introducing the notation si=-ωi, i=1,2,3, one can write that(20)N(s)=(s+ω1)(s+ω2)(s+ω3).Clearly, the system is subcritical, if none of theω1,ω2, and ω3 is zero or negative. For the sake of illustration in Figure 2 the dependence of the characteristic function 𝒩(s) on s is shown (for the calculations, the following parameter values were used: QFA=1/3, QFT=2/3, QF=1, QTA=14/10, QTF=3/5, QT=2, q1(p)=3, q1(d)=0.02, and λ=0.1). The right hand side figure is an enlargement of the 𝒩(s) in the interval (-0.4,0).Characteristic function𝒩(s). The right-hand figure shows the two largest roots of the equation.
(a)(b)The algebraic solutions for the Laplace transformsm̃1(T)(s∣Sj) and m̃1(C)(s∣Sj) can be obtained in a similar manner. Here we only list these solutions, which for m̃1(T)(s∣Sj) are given by(21)m̃1(T)(s∣F)=QFT(s+λ)N(s),m̃1(T)(s∣T)=(s+QF)(s+λ)N(s),m̃1(T)(s∣C)=QFTλN(s).
and for m̃1(C)(s∣Sj) as(22)m̃1(C)(s∣F)=q1(d)QFTQTFN(s),m̃1(C)(s∣T)=q1(d)QTF(s+QF)N(s),m̃1(C)(s∣C)=q1(d)(s+QF)(s+QT)-q1(p)QFTQTFN(s).The expectations can be obtained in a rather simple way by inversion of these Laplace transforms. All solutions consist of the sum of three exponential functions, namely, ofe-ω1t,e-ω2t and e-ω3t. As an illustration, we give the expectation m1(F)(t∣F) of the number of fast neutrons, generated by one initial fast neutron injected into the system. One obtains
(23)m1(F)(t∣F)=ω12-(QT+λ)ω1+QTλ(ω1-ω2)(ω1-ω3)e-ω1t-ω22-(QT+λ)ω2+QTλ(ω1-ω2)(ω2-ω3)e-ω2t+ω32-(QT+λ)ω3+QTλ(ω1-ω3)(ω2-ω3)e-ω3t.
The other expectations are obtained in a similar form, that is, as a sum of three exponentials, and they will not be given here. Some quantitative examples of the expectations are shown in Figure 3. The figure shows the expectations of the numbers of fast and thermal neutrons and the delayed neutron precursors versus time, assuming that the starting particle was either a fast or a thermal neutron (for the calculations the following parameter values were used: QFA=1/3, QFT=2/3, QF=1, QTA=14/10, QTF=3/5, QT=2, q1(p)=3, q1(d)=0.02, and λ=0.1).Dependence of the expectations of the numbers of particles of typesF, T,andC on the time, assuming that the starting particle was either F or T.
(a)(b)(c)It is to be mentioned that the above results could also be obtained directly from deterministic equations, namely, from the two-group point kinetic equations with one group of delayed neutrons.


### 3.2. Variances of the Numbers of Particles of Different Types

As it follows from the definitions and formulae in the previous section, the variance of the numbers of, say, fast neutrons at the timet≥0, induced by an initial fast neutron, is given by the formula:
(24)D2{n1(t)∣F}=m2(F)(t∣F)+m1(F)(t∣F)[1-m1(F)(t∣F)].
Similar expressions can be derived for the other 8 variances.


#### 3.2.1. Second Factorial Moments

It is seen that for the determination of variances, one needs the second factorial moments which can be obtained from the generating function (8)–(10). Introducing the notations:
(25)m2(F)(t∣F)=[∂2g(z1,z2,z3,t∣F)∂z12]z1=z2=z3=1,m2(F)(t∣T)=[∂2g(z1,z2,z3,t∣T)∂z12]z1=z2=z3=1,m2(F)(t∣C)=[∂2g(z1,z2,z3,t∣C)∂z12]z1=z2=z3=1,
one can derive the following equations:(26)m2(F)(t∣F)=QFT∫0te-QF(t-t′)m2(F)(t′∣T)dt′,(27)m2(F)(t∣T)=QTF∫0te-QT(t-t′)×{q1(p)m2(F)(t′∣F)+q2[m1(F)(t′∣F)]2+2q1(p)q1(d)m1(F)(t′∣F)×m1(F)(t′∣C)+q1(d)m2(F)(t′∣C)}dt′,(28)m2(F)(t∣C)=λ∫0te-λ(t-t′)m2(F)(t∣F)dt′.Analogous definitions can be introduced form2(T)(t∣Sj), and m2(C)(t∣Sj) and similar equations can be readily derived for these.For the solution of the arising system of integral equations, again the method of the Laplace transform is used. Before performing the transforms, it is practical to introduce short-hand notations for the functions appearing in (27) and in the corresponding equations for m2(T)(t∣Sj) and m2(C)(t∣Sj) as(29)AF(t′)=q2[m1(F)(t′∣F)]2+2q1(p)q1(d)m1(F)(t′∣F)m1(F)(t′∣C),AT(t′)=q2[m1(T)(t′∣F)]2+2q1(p)q1(d)m1(T)(t′∣F)m1(T)(t′∣C),AC(t′)=q2[m1(C)(t′∣F)]2+2q1(p)q1(d)m1(C)(t′∣F)m1(C)(t′∣C).
These functions are known at this stage, since they contain only the expectations of the numbers of the particles. Applying the notation for the Laplace transforms defined already, one obtains(30)m̃2(F)(s∣F)=QFTs+QFm̃2(F)(s∣T),m̃2(F)(s∣T)=QTFs+QT,×[q1(p)m̃2(F)(s∣F)+q1(d)m̃2(F)(s∣C)+ÃF(s)],m̃2(F)(s∣C)=λs+λm̃2(F)(s∣F),
and similar equations for m̃2(T)(s∣Sj) and m̃2(C)(s∣Sj). After simple algebra the solutions can be written in the following form:
(31)m̃2(F)(s∣F)=QTFQFT(s+λ)N(s)ÃF(s)=QTFm̃1(T)(s∣F)ÃF(s),m̃2(F)(s∣T)=QTF(s+QT)(s+λ)N(s)ÃF(s)=QTFm̃1(T)(s∣T)ÃF(s),m̃2(F)(s∣C)=QTFQFTλN(s)ÃF(s)=QTFm̃1(T)(s∣C)ÃF(s).
Similar solutions are found for the other two groups of factorial moments.Based on the form of the solutions in the Laplace domain, one notices that the solutions in the time domain can be written in the form a convolution as follows:(32)m2(F)(t∣F)=QTF∫0tm1(T)(t-t′∣F)AF(t′)dt′,m2(F)(t∣T)=QTF∫0tm1(T)(t-t′∣T)AF(t′)dt′,m2(F)(t∣C)=QTF∫0tm1(T)(t-t′∣C)AF(t′)dt′,
and similarly for the second factorial moments of the thermal neutrons and the delayed neutron precursors. Equations (32) express the fact, known from the backward theory of branching processes [2], that the first moments of the single-particle generated distributions play the role of the Green’s function for the higher-order moments (and also for all order moments of the distributions of the detected neutrons). This is because the higher-order moment equations have the same form as those for the first moment, except that the Dirac delta function in the first moment equations, representing the starting particle, is replaced by some products of known first moments quantities, which play the role of the inhomogeneous r.h.s. of the second- and higher-order moments. These inhomogeneous right hand sides, or “source functions’’ depend on the problem at hand and the type of the moment to be calculated. In the present case they are given by the functions AF(t), and so forth of (29).The convolution integrals over these known functions can be performed analytically and closed form analytical solutions can be obtained for the second factorial moments. However, in the present case these explicit forms are extremely long and complicated expressions containing the exponential functionse-ω1t,e-ω2t, and e-ω3t in various combinations. These will not be given here since there is very little insight one could gain from the analytical form of the coefficients multiplying the exponentials.


### 3.3. Covariances of the Numbers of Particles of Different Types

Although not needed explicitly for the calculation of the two-group version of the Feynman-alpha formula, it might give some insight to calculate the covariances of the numbers of particles of different types at a given time moment, provided that at the time instantt=0, only one particle was in the system. This will be done in this subsection. If the starting particle was a fast neutron, then the following covariances have to be calculated:
(33)Cov{n1(t),n2(t)∣SF}=m2(FT)(t∣F)-m1(F)(t∣F)m1(T)(t∣F),Cov{n1(t),n3(t)∣SF}=m2(FC)(t∣F)-m1(F)(t∣F)m1(C)(t∣F),Cov{n2(t),n3(t)∣SF}=m2(TC)(t∣F)-m1(T)(t∣F)m1(C)(t∣F).
Along the same lines, the other 6 covariances can also readily be written down. However, for the sake of the simplicity only the above covariances will be calculated.


#### 3.3.1. Mixed Second Moments

In order to determine the covariance between two different random functions at a given time moment, one should calculate first the mixed second moments. In the present case, one needs the following moments:(34)m2(FT)(t∣F)=[∂2g(z1,z2,z3,t∣F)∂z1∂z2]z1=z2=z3=1,m2(FC)(t∣F)=[∂2g(z1,z2,z3,t∣F)∂z1z3]z1=z2=z3=1,m2(TC)(t∣F)=[∂2g(z1,z2,z3,t∣F)∂z2z3]z1=z2=z3=1.The calculations are straightforward and similar to the calculation of the second factorial moments of the previous section, but rather involved and lengthy. Hence the details of the calculations will not be given here. We only note that similar to the case of the second factorial moments, it is practical to introduce a shorthand notation for the functions(35)AFT(t′)=q2m1(F)(t′∣F)m1(T)(t′∣F)+q1(p)q1(d)[m1(F)(t′∣F)m1(T)(t′∣C)+m1(T)(t′∣F)m1(F)(t′∣C)],
and similarly for AFC(t′) and ATC(t′) which play the role of the inhomogeneous part of the mixed moment equations and hence appear in the convolution expressions for the solutions. Without going into details, by using the convolution theorem, the formal solutions for the three mixed moments of the distributions of particles induced by one starting fast neutron are quoted as follows:
(36)m2(FT)(t∣F)=QTF∫0tm1(T)(t-t′∣F)AFT(t′)dt′,m2(FC)(t∣F)=QTF∫0tm1(T)(t-t′∣F)AFC(t′)dt′,m2(TC)(t∣F)=QTF∫0tm1(T)(t-t′∣F)ATC(t′)dt′.By using the formulae (36), one can immediately calculate the covariances (33).Figure4 shows the time dependence of the covariances between the numbers of fast and thermal neutrons, as well as between fast neutrons and delayed neutron precursors and also between the thermal neutrons and the precursors, assuming that the starting particle was a fast neutron. It is remarkable that each of them changes sign, but the absolute values of the covariances are very small. One can also show that the variation of the mean decay time λ-1 influences only slightly the values of covariances.Dependence of the covariances between the numbers of particles of typesF and T, as well as F and C, and also between the numbers of particles of types T and C on the time, assuming that the starting particle was F type.
(a)(b)(c)The reason for the initially negative value of the covariances can be explained in the same way as for the covariance between the fast and thermal neutrons without the presence of delayed neutrons, as was discussed in [14]. Namely, the process is started by one single fast neutron, which is the only particle in the system at t=0, hence the joint expectation of having any two particles is zero at the beginning. The expectation of having a thermal neutron or a delayed neutron precursor is also zero at the beginning, but with the thermalization of the initial fast neutron the expectation of having a thermal neutron starts to deviate from zero, whereas the joint expectation of finding both a fast and a thermal neutron is negligible until the thermal neutron induces fission. Before the branching starts with first thermalisation and then a thermal fission, the covariance is negative. When the branching starts, the joint expectation of any two particles starts to increase, so the covariance starts to increase after having reached a local minimum and thereafter becomes positive.


## 3.1. Expectations of the Numbers of Particles of Different Types

By using (8)–(10), one can derive equations for the expectations of the numbers of fast and thermal neutrons and the delayed neutron precursors, that is, particles of types F, T, and C, respectively. With obvious notations, for the expectation (first moment) m1(F)(t∣Sj) of the number of fast neutrons, induced by one starting fast neutron, thermal neutron and delayed neutron precursor, respectively, one obtains the equations:
(14)m1(F)(t∣F)=[∂g(z1,z2,z3,t∣F)∂z1]z1=z2=z3=1=e-QFt+QFT∫0te-QF(t-t′)m1(F)(t′∣T)dt′,m1(F)(t∣T)=[∂g(z1,z2,z3,t∣T)∂z1]z1=z2=z3=1=QTF∫0te-QT(t-t′)×[q1(p)m1(F)(t′∣F)+q1(d)m1(F)(t′∣C)]dt′,m1(F)(t∣C)=[∂g(z1,z2,z3,t∣C)∂z1]z1=z2=z3=1=λ∫0te-λ(t-t′)m1(F)(t′∣F)dt′.
Equations for the quantities m1(T)(t∣Sj) and m1(C)(t∣Sj) can be derived in a completely similar manner; these will not be given here, for brevity.The arising integral equation system can be readily solved by Laplace transform methods. Introduce the Laplace transforms:(15)m̃1(Sj)(s∣Si)=∫0∞e-stm1(Sj)(t∣Si)dt,
then, from (14) one obtains
(16)m̃1(F)(s∣F)=1s+QF+QFTs+QFm̃1(F)(s∣T),m̃1(F)(s∣T)=QTFs+QT[q1(p)m̃1(F)(s∣F)+q1(d)m̃1(F)(s∣C)],m̃1(F)(s∣C)=λs+λm̃1(F)(s∣F).
After elementary algebra one obtains(17)m̃1(F)(s∣F)=(s+QT)(s+λ)N(s),
where(18)N(s)=(s+QF)(s+QT)(s+λ)-QFTQTF[q1(p)(s+λ)+q1(d)λ]
is a third-order polynomial of s. It can be proven that the roots of the equation:
(19)N(s)=s3+(QF+QT+λ)s2+[QFQT+λ(QF+QT)-q1(p)QFTQTF]s+λ[QFQT-QFTQTF(q1(p)+q1(d))]=0
are all real. Hence, introducing the notation si=-ωi, i=1,2,3, one can write that(20)N(s)=(s+ω1)(s+ω2)(s+ω3).Clearly, the system is subcritical, if none of theω1,ω2, and ω3 is zero or negative. For the sake of illustration in Figure 2 the dependence of the characteristic function 𝒩(s) on s is shown (for the calculations, the following parameter values were used: QFA=1/3, QFT=2/3, QF=1, QTA=14/10, QTF=3/5, QT=2, q1(p)=3, q1(d)=0.02, and λ=0.1). The right hand side figure is an enlargement of the 𝒩(s) in the interval (-0.4,0).Characteristic function𝒩(s). The right-hand figure shows the two largest roots of the equation.
(a)(b)The algebraic solutions for the Laplace transformsm̃1(T)(s∣Sj) and m̃1(C)(s∣Sj) can be obtained in a similar manner. Here we only list these solutions, which for m̃1(T)(s∣Sj) are given by(21)m̃1(T)(s∣F)=QFT(s+λ)N(s),m̃1(T)(s∣T)=(s+QF)(s+λ)N(s),m̃1(T)(s∣C)=QFTλN(s).
and for m̃1(C)(s∣Sj) as(22)m̃1(C)(s∣F)=q1(d)QFTQTFN(s),m̃1(C)(s∣T)=q1(d)QTF(s+QF)N(s),m̃1(C)(s∣C)=q1(d)(s+QF)(s+QT)-q1(p)QFTQTFN(s).The expectations can be obtained in a rather simple way by inversion of these Laplace transforms. All solutions consist of the sum of three exponential functions, namely, ofe-ω1t,e-ω2t and e-ω3t. As an illustration, we give the expectation m1(F)(t∣F) of the number of fast neutrons, generated by one initial fast neutron injected into the system. One obtains
(23)m1(F)(t∣F)=ω12-(QT+λ)ω1+QTλ(ω1-ω2)(ω1-ω3)e-ω1t-ω22-(QT+λ)ω2+QTλ(ω1-ω2)(ω2-ω3)e-ω2t+ω32-(QT+λ)ω3+QTλ(ω1-ω3)(ω2-ω3)e-ω3t.
The other expectations are obtained in a similar form, that is, as a sum of three exponentials, and they will not be given here. Some quantitative examples of the expectations are shown in Figure 3. The figure shows the expectations of the numbers of fast and thermal neutrons and the delayed neutron precursors versus time, assuming that the starting particle was either a fast or a thermal neutron (for the calculations the following parameter values were used: QFA=1/3, QFT=2/3, QF=1, QTA=14/10, QTF=3/5, QT=2, q1(p)=3, q1(d)=0.02, and λ=0.1).Dependence of the expectations of the numbers of particles of typesF, T,andC on the time, assuming that the starting particle was either F or T.
(a)(b)(c)It is to be mentioned that the above results could also be obtained directly from deterministic equations, namely, from the two-group point kinetic equations with one group of delayed neutrons.


## 3.2. Variances of the Numbers of Particles of Different Types

As it follows from the definitions and formulae in the previous section, the variance of the numbers of, say, fast neutrons at the timet≥0, induced by an initial fast neutron, is given by the formula:
(24)D2{n1(t)∣F}=m2(F)(t∣F)+m1(F)(t∣F)[1-m1(F)(t∣F)].
Similar expressions can be derived for the other 8 variances.


### 3.2.1. Second Factorial Moments

It is seen that for the determination of variances, one needs the second factorial moments which can be obtained from the generating function (8)–(10). Introducing the notations:
(25)m2(F)(t∣F)=[∂2g(z1,z2,z3,t∣F)∂z12]z1=z2=z3=1,m2(F)(t∣T)=[∂2g(z1,z2,z3,t∣T)∂z12]z1=z2=z3=1,m2(F)(t∣C)=[∂2g(z1,z2,z3,t∣C)∂z12]z1=z2=z3=1,
one can derive the following equations:(26)m2(F)(t∣F)=QFT∫0te-QF(t-t′)m2(F)(t′∣T)dt′,(27)m2(F)(t∣T)=QTF∫0te-QT(t-t′)×{q1(p)m2(F)(t′∣F)+q2[m1(F)(t′∣F)]2+2q1(p)q1(d)m1(F)(t′∣F)×m1(F)(t′∣C)+q1(d)m2(F)(t′∣C)}dt′,(28)m2(F)(t∣C)=λ∫0te-λ(t-t′)m2(F)(t∣F)dt′.Analogous definitions can be introduced form2(T)(t∣Sj), and m2(C)(t∣Sj) and similar equations can be readily derived for these.For the solution of the arising system of integral equations, again the method of the Laplace transform is used. Before performing the transforms, it is practical to introduce short-hand notations for the functions appearing in (27) and in the corresponding equations for m2(T)(t∣Sj) and m2(C)(t∣Sj) as(29)AF(t′)=q2[m1(F)(t′∣F)]2+2q1(p)q1(d)m1(F)(t′∣F)m1(F)(t′∣C),AT(t′)=q2[m1(T)(t′∣F)]2+2q1(p)q1(d)m1(T)(t′∣F)m1(T)(t′∣C),AC(t′)=q2[m1(C)(t′∣F)]2+2q1(p)q1(d)m1(C)(t′∣F)m1(C)(t′∣C).
These functions are known at this stage, since they contain only the expectations of the numbers of the particles. Applying the notation for the Laplace transforms defined already, one obtains(30)m̃2(F)(s∣F)=QFTs+QFm̃2(F)(s∣T),m̃2(F)(s∣T)=QTFs+QT,×[q1(p)m̃2(F)(s∣F)+q1(d)m̃2(F)(s∣C)+ÃF(s)],m̃2(F)(s∣C)=λs+λm̃2(F)(s∣F),
and similar equations for m̃2(T)(s∣Sj) and m̃2(C)(s∣Sj). After simple algebra the solutions can be written in the following form:
(31)m̃2(F)(s∣F)=QTFQFT(s+λ)N(s)ÃF(s)=QTFm̃1(T)(s∣F)ÃF(s),m̃2(F)(s∣T)=QTF(s+QT)(s+λ)N(s)ÃF(s)=QTFm̃1(T)(s∣T)ÃF(s),m̃2(F)(s∣C)=QTFQFTλN(s)ÃF(s)=QTFm̃1(T)(s∣C)ÃF(s).
Similar solutions are found for the other two groups of factorial moments.Based on the form of the solutions in the Laplace domain, one notices that the solutions in the time domain can be written in the form a convolution as follows:(32)m2(F)(t∣F)=QTF∫0tm1(T)(t-t′∣F)AF(t′)dt′,m2(F)(t∣T)=QTF∫0tm1(T)(t-t′∣T)AF(t′)dt′,m2(F)(t∣C)=QTF∫0tm1(T)(t-t′∣C)AF(t′)dt′,
and similarly for the second factorial moments of the thermal neutrons and the delayed neutron precursors. Equations (32) express the fact, known from the backward theory of branching processes [2], that the first moments of the single-particle generated distributions play the role of the Green’s function for the higher-order moments (and also for all order moments of the distributions of the detected neutrons). This is because the higher-order moment equations have the same form as those for the first moment, except that the Dirac delta function in the first moment equations, representing the starting particle, is replaced by some products of known first moments quantities, which play the role of the inhomogeneous r.h.s. of the second- and higher-order moments. These inhomogeneous right hand sides, or “source functions’’ depend on the problem at hand and the type of the moment to be calculated. In the present case they are given by the functions AF(t), and so forth of (29).The convolution integrals over these known functions can be performed analytically and closed form analytical solutions can be obtained for the second factorial moments. However, in the present case these explicit forms are extremely long and complicated expressions containing the exponential functionse-ω1t,e-ω2t, and e-ω3t in various combinations. These will not be given here since there is very little insight one could gain from the analytical form of the coefficients multiplying the exponentials.


## 3.2.1. Second Factorial Moments

It is seen that for the determination of variances, one needs the second factorial moments which can be obtained from the generating function (8)–(10). Introducing the notations:
(25)m2(F)(t∣F)=[∂2g(z1,z2,z3,t∣F)∂z12]z1=z2=z3=1,m2(F)(t∣T)=[∂2g(z1,z2,z3,t∣T)∂z12]z1=z2=z3=1,m2(F)(t∣C)=[∂2g(z1,z2,z3,t∣C)∂z12]z1=z2=z3=1,
one can derive the following equations:(26)m2(F)(t∣F)=QFT∫0te-QF(t-t′)m2(F)(t′∣T)dt′,(27)m2(F)(t∣T)=QTF∫0te-QT(t-t′)×{q1(p)m2(F)(t′∣F)+q2[m1(F)(t′∣F)]2+2q1(p)q1(d)m1(F)(t′∣F)×m1(F)(t′∣C)+q1(d)m2(F)(t′∣C)}dt′,(28)m2(F)(t∣C)=λ∫0te-λ(t-t′)m2(F)(t∣F)dt′.Analogous definitions can be introduced form2(T)(t∣Sj), and m2(C)(t∣Sj) and similar equations can be readily derived for these.For the solution of the arising system of integral equations, again the method of the Laplace transform is used. Before performing the transforms, it is practical to introduce short-hand notations for the functions appearing in (27) and in the corresponding equations for m2(T)(t∣Sj) and m2(C)(t∣Sj) as(29)AF(t′)=q2[m1(F)(t′∣F)]2+2q1(p)q1(d)m1(F)(t′∣F)m1(F)(t′∣C),AT(t′)=q2[m1(T)(t′∣F)]2+2q1(p)q1(d)m1(T)(t′∣F)m1(T)(t′∣C),AC(t′)=q2[m1(C)(t′∣F)]2+2q1(p)q1(d)m1(C)(t′∣F)m1(C)(t′∣C).
These functions are known at this stage, since they contain only the expectations of the numbers of the particles. Applying the notation for the Laplace transforms defined already, one obtains(30)m̃2(F)(s∣F)=QFTs+QFm̃2(F)(s∣T),m̃2(F)(s∣T)=QTFs+QT,×[q1(p)m̃2(F)(s∣F)+q1(d)m̃2(F)(s∣C)+ÃF(s)],m̃2(F)(s∣C)=λs+λm̃2(F)(s∣F),
and similar equations for m̃2(T)(s∣Sj) and m̃2(C)(s∣Sj). After simple algebra the solutions can be written in the following form:
(31)m̃2(F)(s∣F)=QTFQFT(s+λ)N(s)ÃF(s)=QTFm̃1(T)(s∣F)ÃF(s),m̃2(F)(s∣T)=QTF(s+QT)(s+λ)N(s)ÃF(s)=QTFm̃1(T)(s∣T)ÃF(s),m̃2(F)(s∣C)=QTFQFTλN(s)ÃF(s)=QTFm̃1(T)(s∣C)ÃF(s).
Similar solutions are found for the other two groups of factorial moments.Based on the form of the solutions in the Laplace domain, one notices that the solutions in the time domain can be written in the form a convolution as follows:(32)m2(F)(t∣F)=QTF∫0tm1(T)(t-t′∣F)AF(t′)dt′,m2(F)(t∣T)=QTF∫0tm1(T)(t-t′∣T)AF(t′)dt′,m2(F)(t∣C)=QTF∫0tm1(T)(t-t′∣C)AF(t′)dt′,
and similarly for the second factorial moments of the thermal neutrons and the delayed neutron precursors. Equations (32) express the fact, known from the backward theory of branching processes [2], that the first moments of the single-particle generated distributions play the role of the Green’s function for the higher-order moments (and also for all order moments of the distributions of the detected neutrons). This is because the higher-order moment equations have the same form as those for the first moment, except that the Dirac delta function in the first moment equations, representing the starting particle, is replaced by some products of known first moments quantities, which play the role of the inhomogeneous r.h.s. of the second- and higher-order moments. These inhomogeneous right hand sides, or “source functions’’ depend on the problem at hand and the type of the moment to be calculated. In the present case they are given by the functions AF(t), and so forth of (29).The convolution integrals over these known functions can be performed analytically and closed form analytical solutions can be obtained for the second factorial moments. However, in the present case these explicit forms are extremely long and complicated expressions containing the exponential functionse-ω1t,e-ω2t, and e-ω3t in various combinations. These will not be given here since there is very little insight one could gain from the analytical form of the coefficients multiplying the exponentials.


## 3.3. Covariances of the Numbers of Particles of Different Types

Although not needed explicitly for the calculation of the two-group version of the Feynman-alpha formula, it might give some insight to calculate the covariances of the numbers of particles of different types at a given time moment, provided that at the time instantt=0, only one particle was in the system. This will be done in this subsection. If the starting particle was a fast neutron, then the following covariances have to be calculated:
(33)Cov{n1(t),n2(t)∣SF}=m2(FT)(t∣F)-m1(F)(t∣F)m1(T)(t∣F),Cov{n1(t),n3(t)∣SF}=m2(FC)(t∣F)-m1(F)(t∣F)m1(C)(t∣F),Cov{n2(t),n3(t)∣SF}=m2(TC)(t∣F)-m1(T)(t∣F)m1(C)(t∣F).
Along the same lines, the other 6 covariances can also readily be written down. However, for the sake of the simplicity only the above covariances will be calculated.


### 3.3.1. Mixed Second Moments

In order to determine the covariance between two different random functions at a given time moment, one should calculate first the mixed second moments. In the present case, one needs the following moments:(34)m2(FT)(t∣F)=[∂2g(z1,z2,z3,t∣F)∂z1∂z2]z1=z2=z3=1,m2(FC)(t∣F)=[∂2g(z1,z2,z3,t∣F)∂z1z3]z1=z2=z3=1,m2(TC)(t∣F)=[∂2g(z1,z2,z3,t∣F)∂z2z3]z1=z2=z3=1.The calculations are straightforward and similar to the calculation of the second factorial moments of the previous section, but rather involved and lengthy. Hence the details of the calculations will not be given here. We only note that similar to the case of the second factorial moments, it is practical to introduce a shorthand notation for the functions(35)AFT(t′)=q2m1(F)(t′∣F)m1(T)(t′∣F)+q1(p)q1(d)[m1(F)(t′∣F)m1(T)(t′∣C)+m1(T)(t′∣F)m1(F)(t′∣C)],
and similarly for AFC(t′) and ATC(t′) which play the role of the inhomogeneous part of the mixed moment equations and hence appear in the convolution expressions for the solutions. Without going into details, by using the convolution theorem, the formal solutions for the three mixed moments of the distributions of particles induced by one starting fast neutron are quoted as follows:
(36)m2(FT)(t∣F)=QTF∫0tm1(T)(t-t′∣F)AFT(t′)dt′,m2(FC)(t∣F)=QTF∫0tm1(T)(t-t′∣F)AFC(t′)dt′,m2(TC)(t∣F)=QTF∫0tm1(T)(t-t′∣F)ATC(t′)dt′.By using the formulae (36), one can immediately calculate the covariances (33).Figure4 shows the time dependence of the covariances between the numbers of fast and thermal neutrons, as well as between fast neutrons and delayed neutron precursors and also between the thermal neutrons and the precursors, assuming that the starting particle was a fast neutron. It is remarkable that each of them changes sign, but the absolute values of the covariances are very small. One can also show that the variation of the mean decay time λ-1 influences only slightly the values of covariances.Dependence of the covariances between the numbers of particles of typesF and T, as well as F and C, and also between the numbers of particles of types T and C on the time, assuming that the starting particle was F type.
(a)(b)(c)The reason for the initially negative value of the covariances can be explained in the same way as for the covariance between the fast and thermal neutrons without the presence of delayed neutrons, as was discussed in [14]. Namely, the process is started by one single fast neutron, which is the only particle in the system at t=0, hence the joint expectation of having any two particles is zero at the beginning. The expectation of having a thermal neutron or a delayed neutron precursor is also zero at the beginning, but with the thermalization of the initial fast neutron the expectation of having a thermal neutron starts to deviate from zero, whereas the joint expectation of finding both a fast and a thermal neutron is negligible until the thermal neutron induces fission. Before the branching starts with first thermalisation and then a thermal fission, the covariance is negative. When the branching starts, the joint expectation of any two particles starts to increase, so the covariance starts to increase after having reached a local minimum and thereafter becomes positive.


## 3.3.1. Mixed Second Moments

In order to determine the covariance between two different random functions at a given time moment, one should calculate first the mixed second moments. In the present case, one needs the following moments:(34)m2(FT)(t∣F)=[∂2g(z1,z2,z3,t∣F)∂z1∂z2]z1=z2=z3=1,m2(FC)(t∣F)=[∂2g(z1,z2,z3,t∣F)∂z1z3]z1=z2=z3=1,m2(TC)(t∣F)=[∂2g(z1,z2,z3,t∣F)∂z2z3]z1=z2=z3=1.The calculations are straightforward and similar to the calculation of the second factorial moments of the previous section, but rather involved and lengthy. Hence the details of the calculations will not be given here. We only note that similar to the case of the second factorial moments, it is practical to introduce a shorthand notation for the functions(35)AFT(t′)=q2m1(F)(t′∣F)m1(T)(t′∣F)+q1(p)q1(d)[m1(F)(t′∣F)m1(T)(t′∣C)+m1(T)(t′∣F)m1(F)(t′∣C)],
and similarly for AFC(t′) and ATC(t′) which play the role of the inhomogeneous part of the mixed moment equations and hence appear in the convolution expressions for the solutions. Without going into details, by using the convolution theorem, the formal solutions for the three mixed moments of the distributions of particles induced by one starting fast neutron are quoted as follows:
(36)m2(FT)(t∣F)=QTF∫0tm1(T)(t-t′∣F)AFT(t′)dt′,m2(FC)(t∣F)=QTF∫0tm1(T)(t-t′∣F)AFC(t′)dt′,m2(TC)(t∣F)=QTF∫0tm1(T)(t-t′∣F)ATC(t′)dt′.By using the formulae (36), one can immediately calculate the covariances (33).Figure4 shows the time dependence of the covariances between the numbers of fast and thermal neutrons, as well as between fast neutrons and delayed neutron precursors and also between the thermal neutrons and the precursors, assuming that the starting particle was a fast neutron. It is remarkable that each of them changes sign, but the absolute values of the covariances are very small. One can also show that the variation of the mean decay time λ-1 influences only slightly the values of covariances.Dependence of the covariances between the numbers of particles of typesF and T, as well as F and C, and also between the numbers of particles of types T and C on the time, assuming that the starting particle was F type.
(a)(b)(c)The reason for the initially negative value of the covariances can be explained in the same way as for the covariance between the fast and thermal neutrons without the presence of delayed neutrons, as was discussed in [14]. Namely, the process is started by one single fast neutron, which is the only particle in the system at t=0, hence the joint expectation of having any two particles is zero at the beginning. The expectation of having a thermal neutron or a delayed neutron precursor is also zero at the beginning, but with the thermalization of the initial fast neutron the expectation of having a thermal neutron starts to deviate from zero, whereas the joint expectation of finding both a fast and a thermal neutron is negligible until the thermal neutron induces fission. Before the branching starts with first thermalisation and then a thermal fission, the covariance is negative. When the branching starts, the joint expectation of any two particles starts to increase, so the covariance starts to increase after having reached a local minimum and thereafter becomes positive.


## 4. Moments of the Detection of Neutrons

In order to calculate the variance to mean of detected particles, induced by a stationary extraneous source, two steps remain. One is the introduction of the detection process, which is treated in this section. The second is the introduction of a stationary source, which will be treated in the next section.For brevity, in the forthcoming only the detection of the fast neutrons will be discussed. For the application of the Feynman-alpha method, even in a fast system, presumably the detection of the thermal neutrons will be more practical. However, the calculation goes exactly along the same lines, hence for illustration it is sufficient to discuss the detection of fast neutrons.Denote byQFDdt+o(dt) the probability that a fast neutron is detected in the time interval (t,t+dt). Obviously, the detected neutron is also absorbed, that is, removed from the branching process. Let N(t,u) denote the number of fast neutrons detected in the time interval [t-u,t]. If t<u, then N(t) stands for the number of fast neutron detections in the time interval [0,t].


### 4.1. Detection in Time Interval[0,t]

Define the probabilities:(37)P{N(t)=n∣F}=pF(n,t),P{N(t)=n∣T}=pT(n,t),P{N(t)=n∣C}=pC(n,t),
and introduce the generating functions:
(38)E{zN(t)∣F}=∑n=0∞pF(n,t)zn=gF(z,t),E{zN(t)∣T}=∑n=0∞pT(n,t)zn=gT(z,t),E{zN(t)∣C}=∑n=0∞pC(n,t)zn=gC(z,t).
The equations determining the generating function can be easily written down. One obtains
(39)gF(z,t)=e-QFt+QFA∫0te-QF(t-t′)dt′+zQFD∫0te-QF(t-t′)dt′+QFT∫0te-QF(t-t′)gT(z,t′)dt′,gT(z,t)=e-QTt+QTA∫0te-QF(t-t′)dt′+QFT∫0te-QF(t-t′)qp[gF(z,t′)]qd[gC(z,t′)]dt′,
(40)gC(z,t)=e-λt+λ∫0te-λ(t-t′)gF(z,t′)dt′.
Taking into account the second formula in (11), one has(41)qd[gC(z,t′)]=1-q1(d)+q1(d)gC(z,t′).
The notation applied does not show that the delayed neutrons also participate in the process.


#### 4.1.1. Expectations

The expectations of detected number ofF type particles can be easily calculated by using the formulae:
(42)E{N(t)∣F}=[∂gF(z,t)∂z]z=1=n1(F)(t),E{N(t)∣T}=[∂gT(z,t)∂z]z=1=n1(T)(t),E{N(t)∣C}=[∂gC(z,t)∂z]z=1=n1(C)(t).After simple considerations one obtains the Laplace transforms of the expectations given by(43)ñ1(F)(s)=QFD(s+QT)(s+λ)sN(s),ñ1(T)(s)=QFDQFT[q1(p)(s+λ)+q1(d)λ]sN(s),ñ1(C)(s)=QFDλ(s+QT)sN(s),
where 𝒩(s) is defined by (20). Clearly, these expressions are the Laplace transforms of the following integrals:
(44)n1(F)(t)=QFD∫0tm1(F)(t′∣F)dt′,n1(T)(t)=QFD∫0tm1(F)(t′∣T)dt′,n1(C)(t)=QFD∫0tm1(F)(t′∣C)dt′.
By using the well-known Tauberian theorem [16], from (43) one obtains immediately the expectations of the number of F type particles detected in time interval [0,∞]. They are given by
(45)n1(F)(∞)=QFDQTλω1ω2ω3,n1(T)(∞)=QFDλq1(p)+q1(d)ω1ω2ω3n1(C)(∞)=QFDQTλω1ω2ω3,,
where(46)ω1ω2ω3=λ[QFQT-QFTQTF(q1(p)+q1(d))].
It is worth to note that the total number or detected fast neutrons is the same whether the starting particle is a fast neutron or a delayed neutron precursor. However, if the starting particle is a thermal neutron, then, as expected, one obtains a different expectation for the number of detected fast neutrons.


#### 4.1.2. Second Factorial Moments and the Variances

For the characterization of the detecting process, one needs the variances of the number of the fast neutrons, counted in the time interval[0,t], for the three different types of starting particles. In order to determine the variances, one has to calculate the second factorial moments. From equations (9) and (40), it follows that
(47)n2(F)(t)=QFT∫0te-QF(t-t′)n2(T)(t′)dt′,n2(T)(t)=QTF∫0te-QF(t-t′)[q1(p)n2(F)(t′)+q1(d)n2(C)(t′)]dt′+QTF∫0te-QF(t-t′)×{q2[n1(F)(t′)]2+2q1(p)q1(d)n1(F)(t′)n1(C)(t′)}dt′,n2(C)(t)=λ∫0te-λ(t-t′)n2(F)(t′)dt′.
Introducing the notation(48)BFC(t′)=q2[n1(F)(t′)]2+2q1(p)q1(d)n1(F)(t′)n1(C)(t′),
and performing a Laplace transformation on (47), after the usual algebraic manipulations including the application of the convolution theorem, one obtains the solutions as
(49)n2(F)(t)=QTF∫0tm1(T)(t-t′∣F)BFC(t′)dt′,(50)n2(T)(t)=QTF∫0tm1(T)(t-t′∣T)BFC(t′)dt′,(51)n2(C)(t)=QFT∫0tm1(T)(t-t′∣C)BFC(t′)dt′,With the help of these expressions, the variances of the number of fast neutron detections for the three starting particle types are given by(52)D2{N(t)∣F}=n2(F)(t)+n1(F)(t,∣F)[1-n1(F)(t,∣F)],D2{N(t)∣T}=n2(T)(t)+n1(F)(t,∣T)[1-n1(F)(t,∣T)],D2{N(t)∣C}=n2(C)(t)+n1(F)(t,∣C)[1-n1(F)(t,∣C)].


### 4.2. Detection in the Time Interval[t-u,t],t>u

Since the detection in the time interval[0,t-u) is excluded and only in the time interval [t-u,t] is permitted, one can write that
(53)P{N(t,u)=n∣SF}=pF(n,t,u)=∑j+k+l=n∑n1=0∞∑n2=0∞∑n3=0∞p(n1,n2,n3,t-u∣F)×Z(1)(j,u∣n1)Z(2)(k,u∣n2)Z(3)(l,u∣n3),
provided that at the time moment t=0 one fast neutron was in the system. It is easy to prove that the probabilities Z(1)(j,un1),Z(2)(k,un2), and Z(3)(ℓ,u∣n3) are given by the following formulae:
(54)Z(1)(j,u∣n1)=∑j1+⋯+jn1=j∏i=1n1pF(ji,u),Z(2)(k,u∣n2)=∑k1+⋯+kn2=k∏i=1n2pT(ki,u),Z(3)(l,u∣n3)=∑l1+⋯+ln3=l∏i=1n3pC(li,u),
where pF(ji,u),pT(ki,u), and pC(ℓi,u) are defined by (37). From (54) one can immediately see that
(55)∑j=0∞Z(1)(j,u∣n1)zj=[gF(z,u)]n1,∑k=0∞Z(2)(k,u∣n2)zk=[gT(z,u)]n2,∑l=0∞Z(3)(l,u∣n3)zl=[gC(z,u)]n3.
Thus, the generating function(56)gF(z,t,u)=∑n=0∞pF(n,t,u)zn,t>u
can be written in the form:(57)gF(z,t,u)=g[gF(z,u),gT(z,u),gC(z,u),t-u∣SF].
Since, as it is seen from(8)(58)g(z1,z2,z3,0∣F)=z1,
from (57), one obtains(59)gF(z,u,u)=g[gF(z,u),gT(z,u),gC(z,u),0∣SF]=gF(z,u),
which corresponds to the condition(60)lim⁡t↓ugF(z,t,u)=gF(z,u).
Expressions similar to (57) can be obtained for the cases when the starting particle is a fast neutron or a delayed neutron precursor.


#### 4.2.1. Expectation and Second Factorial Moment

For later use, let us determine the expectation:(61)n1(F)(t,u)=[∂gF(z,t,u)∂z]z=1
of the number of fast neutron detections in the time interval [t-u,t], where t>u, provided that the starting particle was also a fast neutron. By using (57) one obtains
(62)n1(F)(t,u)=m1(F)(t-u∣F)n1(F)(u)+m1(T)(t-u∣F)n1(T)(u)+m1(C)(t-u∣F)n1(C)(u),
where n1(F)(u),n1(T)(u), and n1(C)(u) are defined by (44).The calculation of the second factorial moment:(63)n2(F)(t,u)=[∂2gF(z,t,u)∂z2]z=1
is straightforward, but rather tedious. One finds
(64)n2(F)(t,u)=m1(F)(t-u∣F)n2(F)(u)+m1(T)(t-u∣F)×n2(T)(u)+m1(C)(t-u∣F)n2(C)(u)+m2(F)×(t-u∣F)[n1(F)(u)]2+m2(T)(t-u∣F)×[n1(T)(u)]2+m2(C)(t-u∣F)[n1(C)(u)]2+2m2(FT)(t-u∣F)n1(F)(u)n1(T)(u)+2m2(FC)×n1(F)(u)n1(C)(u)+2m2(TC)(t-u∣F)n1(T)(u)n1(C)(u).
In this expression, the second factorial moments n2(F)(u),n2(T)(u), and n2(C)(u) have been determined already by (49), (50), and (51) respectively hence they are already known. With the substitution of all the known functions, an explicit expression can be obtained which will contain three exponentials with the known exponents. The coefficients multiplying the exponents would take to much space to display and hence will not be shown here.


## 4.1. Detection in Time Interval[0,t]

Define the probabilities:(37)P{N(t)=n∣F}=pF(n,t),P{N(t)=n∣T}=pT(n,t),P{N(t)=n∣C}=pC(n,t),
and introduce the generating functions:
(38)E{zN(t)∣F}=∑n=0∞pF(n,t)zn=gF(z,t),E{zN(t)∣T}=∑n=0∞pT(n,t)zn=gT(z,t),E{zN(t)∣C}=∑n=0∞pC(n,t)zn=gC(z,t).
The equations determining the generating function can be easily written down. One obtains
(39)gF(z,t)=e-QFt+QFA∫0te-QF(t-t′)dt′+zQFD∫0te-QF(t-t′)dt′+QFT∫0te-QF(t-t′)gT(z,t′)dt′,gT(z,t)=e-QTt+QTA∫0te-QF(t-t′)dt′+QFT∫0te-QF(t-t′)qp[gF(z,t′)]qd[gC(z,t′)]dt′,
(40)gC(z,t)=e-λt+λ∫0te-λ(t-t′)gF(z,t′)dt′.
Taking into account the second formula in (11), one has(41)qd[gC(z,t′)]=1-q1(d)+q1(d)gC(z,t′).
The notation applied does not show that the delayed neutrons also participate in the process.


### 4.1.1. Expectations

The expectations of detected number ofF type particles can be easily calculated by using the formulae:
(42)E{N(t)∣F}=[∂gF(z,t)∂z]z=1=n1(F)(t),E{N(t)∣T}=[∂gT(z,t)∂z]z=1=n1(T)(t),E{N(t)∣C}=[∂gC(z,t)∂z]z=1=n1(C)(t).After simple considerations one obtains the Laplace transforms of the expectations given by(43)ñ1(F)(s)=QFD(s+QT)(s+λ)sN(s),ñ1(T)(s)=QFDQFT[q1(p)(s+λ)+q1(d)λ]sN(s),ñ1(C)(s)=QFDλ(s+QT)sN(s),
where 𝒩(s) is defined by (20). Clearly, these expressions are the Laplace transforms of the following integrals:
(44)n1(F)(t)=QFD∫0tm1(F)(t′∣F)dt′,n1(T)(t)=QFD∫0tm1(F)(t′∣T)dt′,n1(C)(t)=QFD∫0tm1(F)(t′∣C)dt′.
By using the well-known Tauberian theorem [16], from (43) one obtains immediately the expectations of the number of F type particles detected in time interval [0,∞]. They are given by
(45)n1(F)(∞)=QFDQTλω1ω2ω3,n1(T)(∞)=QFDλq1(p)+q1(d)ω1ω2ω3n1(C)(∞)=QFDQTλω1ω2ω3,,
where(46)ω1ω2ω3=λ[QFQT-QFTQTF(q1(p)+q1(d))].
It is worth to note that the total number or detected fast neutrons is the same whether the starting particle is a fast neutron or a delayed neutron precursor. However, if the starting particle is a thermal neutron, then, as expected, one obtains a different expectation for the number of detected fast neutrons.


### 4.1.2. Second Factorial Moments and the Variances

For the characterization of the detecting process, one needs the variances of the number of the fast neutrons, counted in the time interval[0,t], for the three different types of starting particles. In order to determine the variances, one has to calculate the second factorial moments. From equations (9) and (40), it follows that
(47)n2(F)(t)=QFT∫0te-QF(t-t′)n2(T)(t′)dt′,n2(T)(t)=QTF∫0te-QF(t-t′)[q1(p)n2(F)(t′)+q1(d)n2(C)(t′)]dt′+QTF∫0te-QF(t-t′)×{q2[n1(F)(t′)]2+2q1(p)q1(d)n1(F)(t′)n1(C)(t′)}dt′,n2(C)(t)=λ∫0te-λ(t-t′)n2(F)(t′)dt′.
Introducing the notation(48)BFC(t′)=q2[n1(F)(t′)]2+2q1(p)q1(d)n1(F)(t′)n1(C)(t′),
and performing a Laplace transformation on (47), after the usual algebraic manipulations including the application of the convolution theorem, one obtains the solutions as
(49)n2(F)(t)=QTF∫0tm1(T)(t-t′∣F)BFC(t′)dt′,(50)n2(T)(t)=QTF∫0tm1(T)(t-t′∣T)BFC(t′)dt′,(51)n2(C)(t)=QFT∫0tm1(T)(t-t′∣C)BFC(t′)dt′,With the help of these expressions, the variances of the number of fast neutron detections for the three starting particle types are given by(52)D2{N(t)∣F}=n2(F)(t)+n1(F)(t,∣F)[1-n1(F)(t,∣F)],D2{N(t)∣T}=n2(T)(t)+n1(F)(t,∣T)[1-n1(F)(t,∣T)],D2{N(t)∣C}=n2(C)(t)+n1(F)(t,∣C)[1-n1(F)(t,∣C)].


## 4.1.1. Expectations

The expectations of detected number ofF type particles can be easily calculated by using the formulae:
(42)E{N(t)∣F}=[∂gF(z,t)∂z]z=1=n1(F)(t),E{N(t)∣T}=[∂gT(z,t)∂z]z=1=n1(T)(t),E{N(t)∣C}=[∂gC(z,t)∂z]z=1=n1(C)(t).After simple considerations one obtains the Laplace transforms of the expectations given by(43)ñ1(F)(s)=QFD(s+QT)(s+λ)sN(s),ñ1(T)(s)=QFDQFT[q1(p)(s+λ)+q1(d)λ]sN(s),ñ1(C)(s)=QFDλ(s+QT)sN(s),
where 𝒩(s) is defined by (20). Clearly, these expressions are the Laplace transforms of the following integrals:
(44)n1(F)(t)=QFD∫0tm1(F)(t′∣F)dt′,n1(T)(t)=QFD∫0tm1(F)(t′∣T)dt′,n1(C)(t)=QFD∫0tm1(F)(t′∣C)dt′.
By using the well-known Tauberian theorem [16], from (43) one obtains immediately the expectations of the number of F type particles detected in time interval [0,∞]. They are given by
(45)n1(F)(∞)=QFDQTλω1ω2ω3,n1(T)(∞)=QFDλq1(p)+q1(d)ω1ω2ω3n1(C)(∞)=QFDQTλω1ω2ω3,,
where(46)ω1ω2ω3=λ[QFQT-QFTQTF(q1(p)+q1(d))].
It is worth to note that the total number or detected fast neutrons is the same whether the starting particle is a fast neutron or a delayed neutron precursor. However, if the starting particle is a thermal neutron, then, as expected, one obtains a different expectation for the number of detected fast neutrons.


## 4.1.2. Second Factorial Moments and the Variances

For the characterization of the detecting process, one needs the variances of the number of the fast neutrons, counted in the time interval[0,t], for the three different types of starting particles. In order to determine the variances, one has to calculate the second factorial moments. From equations (9) and (40), it follows that
(47)n2(F)(t)=QFT∫0te-QF(t-t′)n2(T)(t′)dt′,n2(T)(t)=QTF∫0te-QF(t-t′)[q1(p)n2(F)(t′)+q1(d)n2(C)(t′)]dt′+QTF∫0te-QF(t-t′)×{q2[n1(F)(t′)]2+2q1(p)q1(d)n1(F)(t′)n1(C)(t′)}dt′,n2(C)(t)=λ∫0te-λ(t-t′)n2(F)(t′)dt′.
Introducing the notation(48)BFC(t′)=q2[n1(F)(t′)]2+2q1(p)q1(d)n1(F)(t′)n1(C)(t′),
and performing a Laplace transformation on (47), after the usual algebraic manipulations including the application of the convolution theorem, one obtains the solutions as
(49)n2(F)(t)=QTF∫0tm1(T)(t-t′∣F)BFC(t′)dt′,(50)n2(T)(t)=QTF∫0tm1(T)(t-t′∣T)BFC(t′)dt′,(51)n2(C)(t)=QFT∫0tm1(T)(t-t′∣C)BFC(t′)dt′,With the help of these expressions, the variances of the number of fast neutron detections for the three starting particle types are given by(52)D2{N(t)∣F}=n2(F)(t)+n1(F)(t,∣F)[1-n1(F)(t,∣F)],D2{N(t)∣T}=n2(T)(t)+n1(F)(t,∣T)[1-n1(F)(t,∣T)],D2{N(t)∣C}=n2(C)(t)+n1(F)(t,∣C)[1-n1(F)(t,∣C)].


## 4.2. Detection in the Time Interval[t-u,t],t>u

Since the detection in the time interval[0,t-u) is excluded and only in the time interval [t-u,t] is permitted, one can write that
(53)P{N(t,u)=n∣SF}=pF(n,t,u)=∑j+k+l=n∑n1=0∞∑n2=0∞∑n3=0∞p(n1,n2,n3,t-u∣F)×Z(1)(j,u∣n1)Z(2)(k,u∣n2)Z(3)(l,u∣n3),
provided that at the time moment t=0 one fast neutron was in the system. It is easy to prove that the probabilities Z(1)(j,un1),Z(2)(k,un2), and Z(3)(ℓ,u∣n3) are given by the following formulae:
(54)Z(1)(j,u∣n1)=∑j1+⋯+jn1=j∏i=1n1pF(ji,u),Z(2)(k,u∣n2)=∑k1+⋯+kn2=k∏i=1n2pT(ki,u),Z(3)(l,u∣n3)=∑l1+⋯+ln3=l∏i=1n3pC(li,u),
where pF(ji,u),pT(ki,u), and pC(ℓi,u) are defined by (37). From (54) one can immediately see that
(55)∑j=0∞Z(1)(j,u∣n1)zj=[gF(z,u)]n1,∑k=0∞Z(2)(k,u∣n2)zk=[gT(z,u)]n2,∑l=0∞Z(3)(l,u∣n3)zl=[gC(z,u)]n3.
Thus, the generating function(56)gF(z,t,u)=∑n=0∞pF(n,t,u)zn,t>u
can be written in the form:(57)gF(z,t,u)=g[gF(z,u),gT(z,u),gC(z,u),t-u∣SF].
Since, as it is seen from(8)(58)g(z1,z2,z3,0∣F)=z1,
from (57), one obtains(59)gF(z,u,u)=g[gF(z,u),gT(z,u),gC(z,u),0∣SF]=gF(z,u),
which corresponds to the condition(60)lim⁡t↓ugF(z,t,u)=gF(z,u).
Expressions similar to (57) can be obtained for the cases when the starting particle is a fast neutron or a delayed neutron precursor.


### 4.2.1. Expectation and Second Factorial Moment

For later use, let us determine the expectation:(61)n1(F)(t,u)=[∂gF(z,t,u)∂z]z=1
of the number of fast neutron detections in the time interval [t-u,t], where t>u, provided that the starting particle was also a fast neutron. By using (57) one obtains
(62)n1(F)(t,u)=m1(F)(t-u∣F)n1(F)(u)+m1(T)(t-u∣F)n1(T)(u)+m1(C)(t-u∣F)n1(C)(u),
where n1(F)(u),n1(T)(u), and n1(C)(u) are defined by (44).The calculation of the second factorial moment:(63)n2(F)(t,u)=[∂2gF(z,t,u)∂z2]z=1
is straightforward, but rather tedious. One finds
(64)n2(F)(t,u)=m1(F)(t-u∣F)n2(F)(u)+m1(T)(t-u∣F)×n2(T)(u)+m1(C)(t-u∣F)n2(C)(u)+m2(F)×(t-u∣F)[n1(F)(u)]2+m2(T)(t-u∣F)×[n1(T)(u)]2+m2(C)(t-u∣F)[n1(C)(u)]2+2m2(FT)(t-u∣F)n1(F)(u)n1(T)(u)+2m2(FC)×n1(F)(u)n1(C)(u)+2m2(TC)(t-u∣F)n1(T)(u)n1(C)(u).
In this expression, the second factorial moments n2(F)(u),n2(T)(u), and n2(C)(u) have been determined already by (49), (50), and (51) respectively hence they are already known. With the substitution of all the known functions, an explicit expression can be obtained which will contain three exponentials with the known exponents. The coefficients multiplying the exponents would take to much space to display and hence will not be shown here.


## 4.2.1. Expectation and Second Factorial Moment

For later use, let us determine the expectation:(61)n1(F)(t,u)=[∂gF(z,t,u)∂z]z=1
of the number of fast neutron detections in the time interval [t-u,t], where t>u, provided that the starting particle was also a fast neutron. By using (57) one obtains
(62)n1(F)(t,u)=m1(F)(t-u∣F)n1(F)(u)+m1(T)(t-u∣F)n1(T)(u)+m1(C)(t-u∣F)n1(C)(u),
where n1(F)(u),n1(T)(u), and n1(C)(u) are defined by (44).The calculation of the second factorial moment:(63)n2(F)(t,u)=[∂2gF(z,t,u)∂z2]z=1
is straightforward, but rather tedious. One finds
(64)n2(F)(t,u)=m1(F)(t-u∣F)n2(F)(u)+m1(T)(t-u∣F)×n2(T)(u)+m1(C)(t-u∣F)n2(C)(u)+m2(F)×(t-u∣F)[n1(F)(u)]2+m2(T)(t-u∣F)×[n1(T)(u)]2+m2(C)(t-u∣F)[n1(C)(u)]2+2m2(FT)(t-u∣F)n1(F)(u)n1(T)(u)+2m2(FC)×n1(F)(u)n1(C)(u)+2m2(TC)(t-u∣F)n1(T)(u)n1(C)(u).
In this expression, the second factorial moments n2(F)(u),n2(T)(u), and n2(C)(u) have been determined already by (49), (50), and (51) respectively hence they are already known. With the substitution of all the known functions, an explicit expression can be obtained which will contain three exponentials with the known exponents. The coefficients multiplying the exponents would take to much space to display and hence will not be shown here.


## 5. Process with Randomly Injected Particles

The last step in the derivation of the variance to mean or Feynman-alpha formula is to calculate the first two factorial moments of the detected fast neutrons, induced in a subcritical reactor by a stationary source of fast neutrons. Suppose that at the time instantt=0 there are no particles present in the system, but as time passes fast neutrons appear randomly with a given intensity and initiate branching processes independently of one another. The theory of injection of particles is expounded in the book by Pázsit and Pál [2] for particles of one type. The generalization for particles of three types is straightforward.


### 5.1. Joint Distribution of the Numbers of Particles

In this work we assume that the source events constitute a Poisson point process, that is, that the random time interval between two consecutive injections of fast neutrons follows an exponential distribution with parametersF. This corresponds to the case of an ADS driven by a DD or DT neutron generator in continuous mode. Neutron sources of future ADS will operate with spallation sources and/or in pulsed mode, which have a non-Poisson character. However, the generalisation of the treatment below to non-Poisson processes has been already done in other context (see, e.g., [2, 4, 5]), and the treatment presented in this paper can also be extended to the case of non-Poisson sources in a straightforward way.For the case of a Poisson source, it can be easily shown that the generating function of the probabilityP(n1,n2,n3,t;F) of the event:(65){n1(t)=n1,n2(t)=n2,n3(t)=n3},
provided that at time t=0 there were no particles present in the system, is given by(66)G(z1,z2,z3,t;F)=exp⁡{sF∫0t[g(z1,z2,z3,t′∣F)-1]}.
In general, if the injected particles are of type i, and if the intensity of the injection of particle type i is si, then one has(67)G(z1,z2,z3,t;Si)=exp⁡{si∫0t[g(z1,z2,z3,t′∣Si)-1]}.


#### 5.1.1. Calculation of the Expectations and the Variances

By using the logarithm of the generating function, one can write the formulae:(68)M1(F)(t;F)=[∂ln⁡G(z1,z2,z3,t;F)∂z1]z1=z2=z3=1=sF∫0tm1(F)(t′∣F)dt′,M1(T)(t;F)=[∂ln⁡G(z1,z2,z3,t;F)∂z2]z1=z2=z3=1=sF∫0tm1(T)(t′∣F)dt′,M1(C)(t;F)=[∂ln⁡G(z1,z2,z3,t;F)∂z3]z1=z2=z3=1=sF∫0tm1(C)(t′∣F)dt′,
giving the expectations of the numbers of fast and thermal neutrons and delayed neutron precursors, respectively, at the time moment t≥0, provided that the type of the particles injected into the system was fast neutrons. Similar formulae can be derived for the other six expectations.In the subcritical state, that is, ifω1,ω2, and ω3, are positive real numbers, then the process is asymptotically stationary, consequently one has
(69)M1(F)(∞;F)=sFQTλNst,M1(T)(∞;F)=sFQFTλNst,M1(C)(∞;F)=sFq1(d)QFTQTFNst,
where(70)Nst=ω1ω2ω3=λ[QFQT-(q1(p)+q1(d))QFTQTF].
The rest of the stationary expectations are given by the following formulae:
(71)M1(F)(∞;T)=sT(q1(p)+q1(d))QTFλNst,M1(T)(∞;T)=sTQFλNst,M1(C)(∞;T)=sTq1(d)QFQTFNst,M1(F)(∞;C)=sCQTλNst,M1(T)(∞;C)=sCQFTλNst,M1(C)(∞;C)=sCq1(d)QFQT-q1(p)QFTQTFNst.For the determination of variances of the number of fast neutrons, generated by the three different possible starting particles, one obtains the expressions:(72)VF(t;Si)=[∂2ln⁡G(z1,z2,z3,t;Si)∂z12]z1=z2=z3=1+[∂ln⁡G(z1,z2,z3,t;Si)∂z1]z1=z2=z3=1=si∫0t[m2(F)(t′∣Si)+m1(F)(t′∣Si)]dt′.


### 5.2. Distribution of the Number of Fast Neutrons Detected in a Given Time Interval

If the injection of fast neutrons into the medium is performed by a stationary source emitting particles according to a Poisson process defined by the intensity parametersF, then applying the procedure described in [2], one can prove that generating function of the probability mass function:(73)P{N(t,u)=n∣n1=n2=n3=0}=PF(n,t)
is given by
(74)GF(z,t,u)=GF(z,u)exp⁡{sF∫ut[gF(z,t′,u)-1]dt′}ift>u,GF(z,t)=exp⁡{sF∫0t[gF(z,t′)-1]dt′}ift≤u,
where(75)gF(z,t′,u)=g[gF(z,u),gT(z,u),gC(z,u),t′-u∣F].By using the method described in [2, pages 85-86], one can prove that the improper integral:(76)lim⁡t→∞∫ut[gF(z,t′,u)-1]dt′
exists, consequently the stationary generating function:(77)lim⁡t→∞GF(z,t,u)=GF(st)(z,u)
also exists, if the three roots of the characteristic function (18) are nonnegative. In this case the random function N(t,u) converges in distribution to a random function N(st)(u), if t→∞. Hence, 𝔾F(st)(z,u) is the generating function of N(st)(u), that is,(78)GF(st)(z,u)=E{zN(st)(u)}.
It is useful to rewrite the stationary generating function in the following form:(79)ln⁡GF(st)(z,u)=sF∫0u[gF(z,t)-1]dt+sF∫u∞[gF(z,t,u)-1]dt.
Thus one finds(80)ln⁡GF(st)(z,u)=sF∫0u[gF(z,t)-1]dt+sF∫0∞{g[gF(z,u),gT(z,u),gC(z,u),t∣SF]-1}dt.
If the injected particles are thermal neutrons, then the corresponding formula reads as
(81)ln⁡GT(st)(z,u)=sT∫0u[gT(z,t)-1]dt+sT∫0∞{g[gF(z,u),gT(z,u),gC(z,u),t∣ST]-1}dt.
In the continuation, we will only deal with the case of fast neutron injection.


#### 5.2.1. Calculation of the Expectation and the Variance to Mean in Stationary Case

In the stationary case the expectation of the number of the detection of fast neutrons in the time intervalu is given by(82)NF(st)(u)=[∂ln⁡GF(st)(z,u)∂z]z=1=sF∫0un1(F)(t′)dt′+sF∫u∞n1(F)(t′,u)dt′.
By using (79), one can write
(83)∫u∞n1(F)(t′,u)dt′=n1(F)(u)∫0∞m1(F)(t∣F)dt+n1(T)(u)∫0∞m1(T)(t∣F)dt+n1(C)(u)∫0∞m1(C)(t∣F)dt,
where
(84)∫0∞m1(F)(t∣F)dt=λQTNst,∫0∞m1(T)(t∣F)dt=λQFTNst,∫0∞m1(C)(t∣F)dt=q1(d)QFTQTFNst.
After a lenghty algebra, one obtains(85)NF(st)(u)=sFQFDQTQFQT-QFTQTF(q1(p)+q1(d))u,
which does not contain the delayed neutron precursor decay constant λ.By using the logarithmic generation functionln𝔾F(st)(z,u), one can write down the ratio of the variance to mean of the number of F types particles detected in the time interval u in the form:(86)D2{N(st)(u)}E{N(st)(u)}=VF(st)(u)NF(st)(u)=1+1NF(st)(u)[∂2ln⁡GF(st)(z,u)∂z2]z=1=1+Y(u),
where (87)[∂2ln⁡GF(st)(z,u)∂z2]z=1=sF∫0un2(F)(t′)dt′+sF∫u∞n2(F)(t′,u)dt′,
and n2(F)(t′,u) is given by (64). The Feynman Y(u) function can be written in the traditional form as(88)Y(u)=∑i=13Yi(1-1-e-ωiuωiu).
The coefficients Yi are rather involved functions of the roots ωi as well as the various reaction intensities and the first and second factorial moments of the number of fast neutrons and the delayed neutron precursors per fission. Due to their very extensive character, they are not given here.Figure5 shows a quantitative illustration of the dependence of the variance to mean of the number of fast neutron detections on the detection time u. The left-hand side figure shows the linear-linear, while the right-had side one shows a log-linear plot (for the calculations the following parameter values were used: QFD=1/10, QFA=7/30, QFT=2/3, QTA=14/10, QTF=3/5, q1(p)=3, q1(d)=0.02, λ=0.1).Dependence of the ratio of the variance to mean of the number ofF types particle detections on the detection time u. The left-hand side figure shows the linear-linear, while the right-hand side one the log-linear plot.
(a)(b)In view of the fact that there are three decay constantsωi,i=1,2,3 one could expect a more complicated structure for the variance to mean than what is seen in the left-hand side figure, with for example, two or three different plateaus. Such is the case with the traditional (one-group) Feynman-alpha formula with delayed neutrons, which displays two different plateaus.However, the Feynman-alpha formula will show the different plateaus if the decay constants are sufficiently different and differ by orders of magnitude from each other. Even the traditional (one-group) Feynman formula demonstrates this, since even if six different delayed neutron precursor types are accounted for, there are only two plateaus distinguishable: one corresponding to the prompt neutrons and only one more for all the six delayed neutron groups [17]. In the present two-group treatment, there will be two decay constants associated with the prompt neutrons, and one with the delayed neutrons. As discussed already in [14], quantitatively these two prompt neutron decay constants are not separated sufficiently from each other to make the two decay constants easily observable in the lin-lin plot, for example, by displaying two plateaus. The decay constant corresponding to the delayed neutrons, on the other hand, deviates sufficiently from the prompt decay constants. The plateau corresponding to the delayed neutrons is visible on the plot with logarithmic scale on the time axis (right-hand side figure).


## 5.1. Joint Distribution of the Numbers of Particles

In this work we assume that the source events constitute a Poisson point process, that is, that the random time interval between two consecutive injections of fast neutrons follows an exponential distribution with parametersF. This corresponds to the case of an ADS driven by a DD or DT neutron generator in continuous mode. Neutron sources of future ADS will operate with spallation sources and/or in pulsed mode, which have a non-Poisson character. However, the generalisation of the treatment below to non-Poisson processes has been already done in other context (see, e.g., [2, 4, 5]), and the treatment presented in this paper can also be extended to the case of non-Poisson sources in a straightforward way.For the case of a Poisson source, it can be easily shown that the generating function of the probabilityP(n1,n2,n3,t;F) of the event:(65){n1(t)=n1,n2(t)=n2,n3(t)=n3},
provided that at time t=0 there were no particles present in the system, is given by(66)G(z1,z2,z3,t;F)=exp⁡{sF∫0t[g(z1,z2,z3,t′∣F)-1]}.
In general, if the injected particles are of type i, and if the intensity of the injection of particle type i is si, then one has(67)G(z1,z2,z3,t;Si)=exp⁡{si∫0t[g(z1,z2,z3,t′∣Si)-1]}.


### 5.1.1. Calculation of the Expectations and the Variances

By using the logarithm of the generating function, one can write the formulae:(68)M1(F)(t;F)=[∂ln⁡G(z1,z2,z3,t;F)∂z1]z1=z2=z3=1=sF∫0tm1(F)(t′∣F)dt′,M1(T)(t;F)=[∂ln⁡G(z1,z2,z3,t;F)∂z2]z1=z2=z3=1=sF∫0tm1(T)(t′∣F)dt′,M1(C)(t;F)=[∂ln⁡G(z1,z2,z3,t;F)∂z3]z1=z2=z3=1=sF∫0tm1(C)(t′∣F)dt′,
giving the expectations of the numbers of fast and thermal neutrons and delayed neutron precursors, respectively, at the time moment t≥0, provided that the type of the particles injected into the system was fast neutrons. Similar formulae can be derived for the other six expectations.In the subcritical state, that is, ifω1,ω2, and ω3, are positive real numbers, then the process is asymptotically stationary, consequently one has
(69)M1(F)(∞;F)=sFQTλNst,M1(T)(∞;F)=sFQFTλNst,M1(C)(∞;F)=sFq1(d)QFTQTFNst,
where(70)Nst=ω1ω2ω3=λ[QFQT-(q1(p)+q1(d))QFTQTF].
The rest of the stationary expectations are given by the following formulae:
(71)M1(F)(∞;T)=sT(q1(p)+q1(d))QTFλNst,M1(T)(∞;T)=sTQFλNst,M1(C)(∞;T)=sTq1(d)QFQTFNst,M1(F)(∞;C)=sCQTλNst,M1(T)(∞;C)=sCQFTλNst,M1(C)(∞;C)=sCq1(d)QFQT-q1(p)QFTQTFNst.For the determination of variances of the number of fast neutrons, generated by the three different possible starting particles, one obtains the expressions:(72)VF(t;Si)=[∂2ln⁡G(z1,z2,z3,t;Si)∂z12]z1=z2=z3=1+[∂ln⁡G(z1,z2,z3,t;Si)∂z1]z1=z2=z3=1=si∫0t[m2(F)(t′∣Si)+m1(F)(t′∣Si)]dt′.


## 5.1.1. Calculation of the Expectations and the Variances

By using the logarithm of the generating function, one can write the formulae:(68)M1(F)(t;F)=[∂ln⁡G(z1,z2,z3,t;F)∂z1]z1=z2=z3=1=sF∫0tm1(F)(t′∣F)dt′,M1(T)(t;F)=[∂ln⁡G(z1,z2,z3,t;F)∂z2]z1=z2=z3=1=sF∫0tm1(T)(t′∣F)dt′,M1(C)(t;F)=[∂ln⁡G(z1,z2,z3,t;F)∂z3]z1=z2=z3=1=sF∫0tm1(C)(t′∣F)dt′,
giving the expectations of the numbers of fast and thermal neutrons and delayed neutron precursors, respectively, at the time moment t≥0, provided that the type of the particles injected into the system was fast neutrons. Similar formulae can be derived for the other six expectations.In the subcritical state, that is, ifω1,ω2, and ω3, are positive real numbers, then the process is asymptotically stationary, consequently one has
(69)M1(F)(∞;F)=sFQTλNst,M1(T)(∞;F)=sFQFTλNst,M1(C)(∞;F)=sFq1(d)QFTQTFNst,
where(70)Nst=ω1ω2ω3=λ[QFQT-(q1(p)+q1(d))QFTQTF].
The rest of the stationary expectations are given by the following formulae:
(71)M1(F)(∞;T)=sT(q1(p)+q1(d))QTFλNst,M1(T)(∞;T)=sTQFλNst,M1(C)(∞;T)=sTq1(d)QFQTFNst,M1(F)(∞;C)=sCQTλNst,M1(T)(∞;C)=sCQFTλNst,M1(C)(∞;C)=sCq1(d)QFQT-q1(p)QFTQTFNst.For the determination of variances of the number of fast neutrons, generated by the three different possible starting particles, one obtains the expressions:(72)VF(t;Si)=[∂2ln⁡G(z1,z2,z3,t;Si)∂z12]z1=z2=z3=1+[∂ln⁡G(z1,z2,z3,t;Si)∂z1]z1=z2=z3=1=si∫0t[m2(F)(t′∣Si)+m1(F)(t′∣Si)]dt′.


## 5.2. Distribution of the Number of Fast Neutrons Detected in a Given Time Interval

If the injection of fast neutrons into the medium is performed by a stationary source emitting particles according to a Poisson process defined by the intensity parametersF, then applying the procedure described in [2], one can prove that generating function of the probability mass function:(73)P{N(t,u)=n∣n1=n2=n3=0}=PF(n,t)
is given by
(74)GF(z,t,u)=GF(z,u)exp⁡{sF∫ut[gF(z,t′,u)-1]dt′}ift>u,GF(z,t)=exp⁡{sF∫0t[gF(z,t′)-1]dt′}ift≤u,
where(75)gF(z,t′,u)=g[gF(z,u),gT(z,u),gC(z,u),t′-u∣F].By using the method described in [2, pages 85-86], one can prove that the improper integral:(76)lim⁡t→∞∫ut[gF(z,t′,u)-1]dt′
exists, consequently the stationary generating function:(77)lim⁡t→∞GF(z,t,u)=GF(st)(z,u)
also exists, if the three roots of the characteristic function (18) are nonnegative. In this case the random function N(t,u) converges in distribution to a random function N(st)(u), if t→∞. Hence, 𝔾F(st)(z,u) is the generating function of N(st)(u), that is,(78)GF(st)(z,u)=E{zN(st)(u)}.
It is useful to rewrite the stationary generating function in the following form:(79)ln⁡GF(st)(z,u)=sF∫0u[gF(z,t)-1]dt+sF∫u∞[gF(z,t,u)-1]dt.
Thus one finds(80)ln⁡GF(st)(z,u)=sF∫0u[gF(z,t)-1]dt+sF∫0∞{g[gF(z,u),gT(z,u),gC(z,u),t∣SF]-1}dt.
If the injected particles are thermal neutrons, then the corresponding formula reads as
(81)ln⁡GT(st)(z,u)=sT∫0u[gT(z,t)-1]dt+sT∫0∞{g[gF(z,u),gT(z,u),gC(z,u),t∣ST]-1}dt.
In the continuation, we will only deal with the case of fast neutron injection.


### 5.2.1. Calculation of the Expectation and the Variance to Mean in Stationary Case

In the stationary case the expectation of the number of the detection of fast neutrons in the time intervalu is given by(82)NF(st)(u)=[∂ln⁡GF(st)(z,u)∂z]z=1=sF∫0un1(F)(t′)dt′+sF∫u∞n1(F)(t′,u)dt′.
By using (79), one can write
(83)∫u∞n1(F)(t′,u)dt′=n1(F)(u)∫0∞m1(F)(t∣F)dt+n1(T)(u)∫0∞m1(T)(t∣F)dt+n1(C)(u)∫0∞m1(C)(t∣F)dt,
where
(84)∫0∞m1(F)(t∣F)dt=λQTNst,∫0∞m1(T)(t∣F)dt=λQFTNst,∫0∞m1(C)(t∣F)dt=q1(d)QFTQTFNst.
After a lenghty algebra, one obtains(85)NF(st)(u)=sFQFDQTQFQT-QFTQTF(q1(p)+q1(d))u,
which does not contain the delayed neutron precursor decay constant λ.By using the logarithmic generation functionln𝔾F(st)(z,u), one can write down the ratio of the variance to mean of the number of F types particles detected in the time interval u in the form:(86)D2{N(st)(u)}E{N(st)(u)}=VF(st)(u)NF(st)(u)=1+1NF(st)(u)[∂2ln⁡GF(st)(z,u)∂z2]z=1=1+Y(u),
where (87)[∂2ln⁡GF(st)(z,u)∂z2]z=1=sF∫0un2(F)(t′)dt′+sF∫u∞n2(F)(t′,u)dt′,
and n2(F)(t′,u) is given by (64). The Feynman Y(u) function can be written in the traditional form as(88)Y(u)=∑i=13Yi(1-1-e-ωiuωiu).
The coefficients Yi are rather involved functions of the roots ωi as well as the various reaction intensities and the first and second factorial moments of the number of fast neutrons and the delayed neutron precursors per fission. Due to their very extensive character, they are not given here.Figure5 shows a quantitative illustration of the dependence of the variance to mean of the number of fast neutron detections on the detection time u. The left-hand side figure shows the linear-linear, while the right-had side one shows a log-linear plot (for the calculations the following parameter values were used: QFD=1/10, QFA=7/30, QFT=2/3, QTA=14/10, QTF=3/5, q1(p)=3, q1(d)=0.02, λ=0.1).Dependence of the ratio of the variance to mean of the number ofF types particle detections on the detection time u. The left-hand side figure shows the linear-linear, while the right-hand side one the log-linear plot.
(a)(b)In view of the fact that there are three decay constantsωi,i=1,2,3 one could expect a more complicated structure for the variance to mean than what is seen in the left-hand side figure, with for example, two or three different plateaus. Such is the case with the traditional (one-group) Feynman-alpha formula with delayed neutrons, which displays two different plateaus.However, the Feynman-alpha formula will show the different plateaus if the decay constants are sufficiently different and differ by orders of magnitude from each other. Even the traditional (one-group) Feynman formula demonstrates this, since even if six different delayed neutron precursor types are accounted for, there are only two plateaus distinguishable: one corresponding to the prompt neutrons and only one more for all the six delayed neutron groups [17]. In the present two-group treatment, there will be two decay constants associated with the prompt neutrons, and one with the delayed neutrons. As discussed already in [14], quantitatively these two prompt neutron decay constants are not separated sufficiently from each other to make the two decay constants easily observable in the lin-lin plot, for example, by displaying two plateaus. The decay constant corresponding to the delayed neutrons, on the other hand, deviates sufficiently from the prompt decay constants. The plateau corresponding to the delayed neutrons is visible on the plot with logarithmic scale on the time axis (right-hand side figure).


## 5.2.1. Calculation of the Expectation and the Variance to Mean in Stationary Case

In the stationary case the expectation of the number of the detection of fast neutrons in the time intervalu is given by(82)NF(st)(u)=[∂ln⁡GF(st)(z,u)∂z]z=1=sF∫0un1(F)(t′)dt′+sF∫u∞n1(F)(t′,u)dt′.
By using (79), one can write
(83)∫u∞n1(F)(t′,u)dt′=n1(F)(u)∫0∞m1(F)(t∣F)dt+n1(T)(u)∫0∞m1(T)(t∣F)dt+n1(C)(u)∫0∞m1(C)(t∣F)dt,
where
(84)∫0∞m1(F)(t∣F)dt=λQTNst,∫0∞m1(T)(t∣F)dt=λQFTNst,∫0∞m1(C)(t∣F)dt=q1(d)QFTQTFNst.
After a lenghty algebra, one obtains(85)NF(st)(u)=sFQFDQTQFQT-QFTQTF(q1(p)+q1(d))u,
which does not contain the delayed neutron precursor decay constant λ.By using the logarithmic generation functionln𝔾F(st)(z,u), one can write down the ratio of the variance to mean of the number of F types particles detected in the time interval u in the form:(86)D2{N(st)(u)}E{N(st)(u)}=VF(st)(u)NF(st)(u)=1+1NF(st)(u)[∂2ln⁡GF(st)(z,u)∂z2]z=1=1+Y(u),
where (87)[∂2ln⁡GF(st)(z,u)∂z2]z=1=sF∫0un2(F)(t′)dt′+sF∫u∞n2(F)(t′,u)dt′,
and n2(F)(t′,u) is given by (64). The Feynman Y(u) function can be written in the traditional form as(88)Y(u)=∑i=13Yi(1-1-e-ωiuωiu).
The coefficients Yi are rather involved functions of the roots ωi as well as the various reaction intensities and the first and second factorial moments of the number of fast neutrons and the delayed neutron precursors per fission. Due to their very extensive character, they are not given here.Figure5 shows a quantitative illustration of the dependence of the variance to mean of the number of fast neutron detections on the detection time u. The left-hand side figure shows the linear-linear, while the right-had side one shows a log-linear plot (for the calculations the following parameter values were used: QFD=1/10, QFA=7/30, QFT=2/3, QTA=14/10, QTF=3/5, q1(p)=3, q1(d)=0.02, λ=0.1).Dependence of the ratio of the variance to mean of the number ofF types particle detections on the detection time u. The left-hand side figure shows the linear-linear, while the right-hand side one the log-linear plot.
(a)(b)In view of the fact that there are three decay constantsωi,i=1,2,3 one could expect a more complicated structure for the variance to mean than what is seen in the left-hand side figure, with for example, two or three different plateaus. Such is the case with the traditional (one-group) Feynman-alpha formula with delayed neutrons, which displays two different plateaus.However, the Feynman-alpha formula will show the different plateaus if the decay constants are sufficiently different and differ by orders of magnitude from each other. Even the traditional (one-group) Feynman formula demonstrates this, since even if six different delayed neutron precursor types are accounted for, there are only two plateaus distinguishable: one corresponding to the prompt neutrons and only one more for all the six delayed neutron groups [17]. In the present two-group treatment, there will be two decay constants associated with the prompt neutrons, and one with the delayed neutrons. As discussed already in [14], quantitatively these two prompt neutron decay constants are not separated sufficiently from each other to make the two decay constants easily observable in the lin-lin plot, for example, by displaying two plateaus. The decay constant corresponding to the delayed neutrons, on the other hand, deviates sufficiently from the prompt decay constants. The plateau corresponding to the delayed neutrons is visible on the plot with logarithmic scale on the time axis (right-hand side figure).


## 6. Conclusions

The variance to mean formula was derived in a two-group treatment with one group of delayed neutrons with the use of the backward master equation technique. The temporal behaviour of both the first and second factorial moment of the detected particles is determined by three exponentials. The various factorial moments can be fully determined analytically; however, the expressions in most cases are too lengthy to write them out. Qualitatively, the form of the solutions is the same as in the traditional case, but the two decay constants associated with the prompt response of the system, are not distinguishable in the plots. The relationship between the decay constants and the subcritical reactivity of the system will be investigated in future communications.



---

*Source: 620808-2012-05-16.xml*