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.

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 [79], 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 [1012]. 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. Figure 1 illustrates the possible reactions induced by the fast and thermal neutrons and the delayed neutron precursors, respectively.

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 notations F, 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 𝐷𝜈 for the latter. In the literature on neutron fluctuations, the transition probabilities (more correctly, transition intensities) are usually denoted by 𝜆𝑖, where the subscript 𝑖 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) where 𝑄𝐹𝐴 and 𝑄𝐹𝑇 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 by 𝑓(𝑘,) the probability that a thermal neutron produces 𝑘0 fast neutrons and 0 delayed neutron precursors, that is, particles of type C. Suppose that𝑓(𝑘,)=𝑓𝑝(𝑘)𝑓𝑑(),(2) that is, the numbers fast neutrons and delayed neutron precursors created in one reaction are independent. Further, let 𝑓𝑑() be given by 𝑓𝑑()=1𝑞1(𝑑)𝑞,if=0,1(𝑑),if=1,0,if>1,(3) where 𝑞1(𝑑)1. The probability that a delayed neutron precursor produced at time 𝑡=0 decays to a fast neutron during the time interval not larger than 𝑡0 is given by𝑇𝑑=1𝑒𝜆𝑡,𝜆0.(4) 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 functions 𝐧𝐹(𝑡),𝐧𝑇(𝑡), and 𝐧𝐶(𝑡), giving the numbers of particles of types F, T, and C, respectively, at the time moment 𝑡0. Define the probabilities:𝒫𝐧𝐹(𝑡)=𝑛1,𝐧𝑇(𝑡)=𝑛2,𝐧𝐶(𝑡)=𝑛3𝑆𝑗𝑛=𝑝1,𝑛2,𝑛3,𝑡𝑆𝑗,(5) where the conditions 𝑆𝑗, 𝑗={1,2,3},{𝑆𝑗}={𝐹,𝑇,𝐶}, indicating the type of particle starting the process, are defined as 𝑆1𝐧𝐹=𝐹(0)=1,𝐧𝑇(0)=0,𝐧𝐶,𝑆(0)=02𝐧𝑇=𝐹(0)=0,𝐧𝑇(0)=1,𝐧𝐶(,𝑆0)=03𝐧𝐶=𝐹(0)=0,𝐧𝑇(0)=0,𝐧𝐶,(0)=1(6) respectively. Introduce also the corresponding generating functions: 𝑔𝑧1,𝑧2,𝑧3,𝑡𝑆𝑗=𝑛1=0𝑛2=0𝑛3=0𝑝𝑛1,𝑛2,𝑛3,𝑡𝑆𝑗𝑧𝑛11𝑧𝑛22𝑧𝑛33.(7) 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,𝑡), and the summing up of the probabilities of the mutually exclusive events generated by the respective reactions, one can write that 𝑔𝑧1,𝑧2,𝑧3,𝑡𝐹=𝑒𝑄𝐹𝑡𝑧1+𝑄𝐹𝐴𝑡0𝑒𝑄𝐹(𝑡𝑡)𝑑𝑡+𝑄𝐹𝑇𝑡0𝑒𝑄𝐹(𝑡𝑡)𝑔𝑧1,𝑧2,𝑧3,𝑡𝑇𝑑𝑡,𝑔𝑧(8)1,𝑧2,𝑧3,𝑡𝑇=𝑒𝑄𝑇𝑡𝑧2+𝑄𝑇𝐴𝑡0𝑒𝑄𝑇(𝑡𝑡)𝑑𝑡+𝑄𝑇𝐹𝑡0𝑒𝑄𝑇(𝑡𝑡)𝑞𝑝𝑔𝑧1,𝑧2,𝑧3,𝑡𝐹×𝑞𝑑𝑔𝑧1,𝑧2,𝑧3,𝑡𝐶𝑑𝑡,(9) as well as𝑔𝑧1,𝑧2,𝑧3,𝑡𝐶=𝑒𝜆𝑡𝑧3+𝜆𝑡0𝑒𝜆(𝑡𝑡)𝑔𝑧1,𝑧2,𝑧3,𝑡𝐹𝑑𝑡,(10) where𝑞𝑝(𝑧)=𝑘=0𝑓𝑝(𝑘)𝑧𝑘,𝑞𝑑(𝑧)==0𝑓𝑑()=1𝑞1(𝑑)+𝑞1(𝑑)𝑧.(11) For later use, we note that𝑑𝑗𝑞𝑝(𝑧)𝑑𝑧𝑗𝑧=1=𝑞𝑗(𝑝),𝑑𝑗𝑞𝑑(𝑧)𝑑𝑧𝑗𝑧=1=𝑞𝑗(𝑑)=𝛿𝑗1𝑞1(𝑑)(12) 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𝑞1(𝑝)=𝜈𝑝=(1𝛽)𝜈,𝑞1(𝑑)=𝜈𝑑=𝛽𝜈,(13) with 𝜈=𝜈𝑝+𝜈𝑑, and 𝛽=𝜈𝑑/(𝜈𝑝+𝜈𝑑).

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) 𝑚1(𝐹)(𝑡𝑆𝑗) of the number of fast neutrons, induced by one starting fast neutron, thermal neutron and delayed neutron precursor, respectively, one obtains the equations: 𝑚1(𝐹)𝑧(𝑡𝐹)=𝜕𝑔1,𝑧2,𝑧3,𝑡𝐹𝜕𝑧1𝑧1=𝑧2=𝑧3=1=𝑒𝑄𝐹𝑡+𝑄𝐹𝑇𝑡0𝑒𝑄𝐹(𝑡𝑡)𝑚1(𝐹)𝑡𝑚𝑇𝑑𝑡,1(𝐹)𝑧(𝑡𝑇)=𝜕𝑔1,𝑧2,𝑧3,𝑡𝑇𝜕𝑧1𝑧1=𝑧2=𝑧3=1=𝑄𝑇𝐹𝑡0𝑒𝑄𝑇(𝑡𝑡)×𝑞1(𝑝)𝑚1(𝐹)𝑡𝐹+𝑞1(𝑑)𝑚1(𝐹)𝑡𝐶𝑑𝑡,𝑚1(𝐹)(𝑡𝐶)=𝜕𝑔(𝑧1,𝑧2,𝑧3,𝑡𝐶)𝜕𝑧1𝑧1=𝑧2=𝑧3=1=𝜆𝑡0𝑒𝜆(𝑡𝑡)𝑚1(𝐹)𝑡𝐹𝑑𝑡.(14) Equations for the quantities 𝑚1(𝑇)(𝑡𝑆𝑗) and 𝑚1(𝐶)(𝑡𝑆𝑗) 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:𝑚(𝑆𝑗)1𝑠𝑆𝑖=0𝑒𝑠𝑡𝑚(𝑆𝑗)1𝑡𝑆𝑖𝑑𝑡,(15) then, from (14) one obtains 𝑚1(𝐹)1(𝑠𝐹)=𝑠+𝑄𝐹+𝑄𝐹𝑇𝑠+𝑄𝐹𝑚1(𝐹)(𝑠𝑇),𝑚1(𝐹)𝑄(𝑠𝑇)=𝑇𝐹𝑠+𝑄𝑇𝑞1(𝑝)𝑚1(𝐹)(𝑠𝐹)+𝑞1(𝑑)𝑚1(𝐹),(𝑠𝐶)𝑚1(𝐹)𝜆(𝑠𝐶)=𝑠+𝜆𝑚1(𝐹)(𝑠𝐹).(16) After elementary algebra one obtains𝑚1(𝐹)(𝑠𝐹)=𝑠+𝑄𝑇(𝑠+𝜆)𝒩(𝑠),(17) where𝒩(𝑠)=𝑠+𝑄𝐹𝑠+𝑄𝑇(𝑠+𝜆)𝑄𝐹𝑇𝑄𝑇𝐹𝑞1(𝑝)(𝑠+𝜆)+𝑞1(𝑑)𝜆(18) is a third-order polynomial of 𝑠. It can be proven that the roots of the equation: 𝒩(𝑠)=𝑠3+𝑄𝐹+𝑄𝑇𝑠+𝜆2+𝑄𝐹𝑄𝑇𝑄+𝜆𝐹+𝑄𝑇𝑞1(𝑝)𝑄𝐹𝑇𝑄𝑇𝐹𝑠𝑄+𝜆𝐹𝑄𝑇𝑄𝐹𝑇𝑄𝑇𝐹𝑞1(𝑝)+𝑞1(𝑑)=0(19) are all real. Hence, introducing the notation 𝑠𝑖=𝜔𝑖, 𝑖=1,2,3, one can write that𝒩(𝑠)=𝑠+𝜔1𝑠+𝜔2𝑠+𝜔3.(20)

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 𝒩(𝑠) on 𝑠 is shown (for the calculations, the following parameter values were used: 𝑄𝐹𝐴=1/3, 𝑄𝐹𝑇=2/3, 𝑄𝐹=1, 𝑄𝑇𝐴=14/10, 𝑄𝑇𝐹=3/5, 𝑄𝑇=2, q1(𝑝)=3, 𝑞1(𝑑)=0.02, and 𝜆=0.1). The right hand side figure is an enlargement of the 𝒩(𝑠) in the interval (0.4,0).

The algebraic solutions for the Laplace transforms 𝑚1(𝑇)(𝑠𝑆𝑗) and 𝑚1(𝐶)(𝑠𝑆𝑗) can be obtained in a similar manner. Here we only list these solutions, which for 𝑚1(𝑇)(𝑠𝑆𝑗) are given by𝑚1(𝑇)𝑄(𝑠𝐹)=𝐹𝑇(𝑠+𝜆),𝒩(𝑠)𝑚1(𝑇)(𝑠𝑇)=𝑠+𝑄𝐹(𝑠+𝜆),𝒩(𝑠)𝑚1(𝑇)𝑄(𝑠𝐶)=𝐹𝑇𝜆𝒩.(𝑠)(21) and for 𝑚1(𝐶)(𝑠𝑆𝑗) as𝑚1(𝐶)(𝑠𝐹)=𝑞1(𝑑)𝑄𝐹𝑇𝑄𝑇𝐹,𝒩(𝑠)𝑚1(𝐶)(𝑠𝑇)=𝑞1(𝑑)𝑄𝑇𝐹𝑠+𝑄𝐹,𝒩(𝑠)𝑚1(𝐶)(𝑠𝐶)=𝑞1(𝑑)𝑠+𝑄𝐹𝑠+𝑄𝑇𝑞1(𝑝)𝑄𝐹𝑇𝑄𝑇𝐹.𝒩(𝑠)(22)

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, of 𝑒𝜔1𝑡,𝑒𝜔2𝑡 and 𝑒𝜔3𝑡. As an illustration, we give the expectation 𝑚1(𝐹)(𝑡𝐹) of the number of fast neutrons, generated by one initial fast neutron injected into the system. One obtains 𝑚1(𝐹)𝜔(𝑡𝐹)=21𝑄𝑇𝜔+𝜆1+𝑄𝑇𝜆𝜔1𝜔2𝜔1𝜔3𝑒𝜔1𝑡𝜔22𝑄𝑇𝜔+𝜆2+𝑄𝑇𝜆𝜔1𝜔2𝜔2𝜔3𝑒𝜔2𝑡+𝜔23𝑄𝑇𝜔+𝜆3+𝑄𝑇𝜆𝜔1𝜔3𝜔2𝜔3𝑒𝜔3𝑡.(23) 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: Q𝐹𝐴=1/3, 𝑄𝐹𝑇=2/3, 𝑄𝐹=1, 𝑄𝑇𝐴=14/10, 𝑄𝑇𝐹=3/5, 𝑄𝑇=2, 𝑞1(𝑝)=3, 𝑞1(𝑑)=0.02, and 𝜆=0.1).

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 time 𝑡0, induced by an initial fast neutron, is given by the formula: 𝐃2𝐧1(𝑡)𝐹=𝑚2(𝐹)(𝑡𝐹)+𝑚1(𝐹)(𝑡𝐹)1𝑚1(𝐹)(𝑡𝐹).(24) 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: 𝑚2(𝐹)𝜕(𝑡𝐹)=2𝑔𝑧1,𝑧2,𝑧3,𝑡𝐹𝜕𝑧21𝑧1=𝑧2=𝑧3=1,𝑚2(𝐹)𝜕(𝑡𝑇)=2𝑔𝑧1,𝑧2,𝑧3,𝑡𝑇𝜕𝑧21𝑧1=𝑧2=𝑧3=1,𝑚2(𝐹)𝜕(𝑡𝐶)=2𝑔𝑧1,𝑧2,𝑧3,𝑡𝐶𝜕𝑧21𝑧1=𝑧2=𝑧3=1,(25) one can derive the following equations:𝑚2(𝐹)(𝑡𝐹)=𝑄𝐹𝑇𝑡0𝑒𝑄𝐹(𝑡𝑡)𝑚2(𝐹)𝑡𝑚𝑇𝑑𝑡,(26)2(𝐹)𝑄(𝑡𝑇)=𝑇𝐹𝑡0𝑒𝑄𝑇(𝑡𝑡)×𝑞1(𝑝)𝑚2(𝐹)𝑡𝐹+𝑞2𝑚1(𝐹)𝑡𝐹2+2𝑞1(𝑝)𝑞1(𝑑)𝑚1(𝐹)𝑡𝐹×𝑚1(𝐹)𝑡𝐶+𝑞1(𝑑)𝑚2(𝐹)𝑡𝐶𝑑𝑡,𝑚(27)2(𝐹)(𝑡𝐶)=𝜆𝑡0𝑒𝜆(𝑡𝑡)𝑚2(𝐹)(𝑡𝐹)𝑑𝑡.(28)

Analogous definitions can be introduced for 𝑚2(𝑇)(𝑡𝑆𝑗), and 𝑚2(𝐶)(𝑡𝑆𝑗) 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 𝑚2(𝑇)(𝑡𝑆𝑗) and 𝑚2(𝐶)(𝑡𝑆𝑗) as𝐴𝐹𝑡=𝑞2𝑚1(𝐹)𝑡𝐹2+2𝑞1(𝑝)𝑞1(𝑑)𝑚1(𝐹)𝑡𝑚𝐹1(𝐹)𝑡,𝐴𝐶𝑇𝑡=𝑞2𝑚1(𝑇)𝑡𝐹2+2𝑞1(𝑝)𝑞1(𝑑)𝑚1(𝑇)𝑡𝑚𝐹1(𝑇)𝑡,𝐴𝐶𝐶𝑡=𝑞2𝑚1(𝐶)𝑡𝐹2+2𝑞1(𝑝)𝑞1(𝑑)𝑚1(𝐶)𝑡𝑚𝐹1(𝐶)𝑡.𝐶(29) 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𝑚2(𝐹)𝑄(𝑠𝐹)=𝐹𝑇𝑠+𝑄𝐹𝑚2(𝐹)(𝑠𝑇),𝑚2(𝐹)𝑄(𝑠𝑇)=𝑇𝐹𝑠+𝑄𝑇,×𝑞1(𝑝)𝑚2(𝐹)(𝑠𝐹)+𝑞1(𝑑)𝑚2(𝐹)𝐴(𝑠𝐶)+𝐹,(𝑠)𝑚2(𝐹)𝜆(𝑠𝐶)=𝑠+𝜆𝑚2(𝐹)(𝑠𝐹),(30) and similar equations for 𝑚2(𝑇)(𝑠𝑆𝑗) and 𝑚2(𝐶)(𝑠𝑆𝑗). After simple algebra the solutions can be written in the following form: 𝑚2(𝐹)(𝑠𝐹)=𝑄𝑇𝐹𝑄𝐹𝑇(𝑠+𝜆)𝐴𝒩(𝑠)𝐹(𝑠)=𝑄𝑇𝐹𝑚1(𝑇)𝐴(𝑠𝐹)𝐹(𝑠),𝑚2(𝐹)(𝑠𝑇)=𝑄𝑇𝐹𝑠+𝑄𝑇(𝑠+𝜆)𝐴𝒩(𝑠)𝐹(𝑠)=𝑄𝑇𝐹𝑚1(𝑇)(𝐴𝑠𝑇)𝐹(𝑠),𝑚2(𝐹)(𝑠𝐶)=𝑄𝑇𝐹𝑄𝐹𝑇𝜆𝐴𝒩(𝑠)𝐹(𝑠)=𝑄𝑇𝐹𝑚1(𝑇)𝐴(𝑠𝐶)𝐹(𝑠).(31) 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: 𝑚2(𝐹)(𝑡𝐹)=𝑄𝑇𝐹𝑡0𝑚1(𝑇)𝑡𝑡𝐴𝐹𝐹𝑡𝑚𝑑𝑡,2(𝐹)(𝑡𝑇)=𝑄𝑇𝐹𝑡0𝑚1(𝑇)𝑡𝑡𝐴𝑇𝐹𝑡𝑚𝑑𝑡,2(𝐹)(𝑡𝐶)=𝑄𝑇𝐹𝑡0𝑚1(𝑇)𝑡𝑡𝐴𝐶𝐹𝑡𝑑𝑡,(32) 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 𝐴𝐹(𝑡), 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 functions 𝑒𝜔1𝑡,𝑒𝜔2𝑡, and 𝑒𝜔3𝑡 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 instant 𝑡=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: 𝐧𝐂𝐨𝐯1(𝑡),𝐧2(𝑡)𝑆𝐹=𝑚2(𝐹𝑇)(𝑡𝐹)𝑚1(𝐹)(𝑡𝐹)𝑚1(𝑇)𝐧(𝑡𝐹),𝐂𝐨𝐯1(𝑡),𝐧3(𝑡)𝑆𝐹=𝑚2(𝐹𝐶)(𝑡𝐹)𝑚1(𝐹)(𝑡𝐹)𝑚1(𝐶)𝐧(𝑡𝐹),𝐂𝐨𝐯2(𝑡),𝐧3(𝑡)𝑆𝐹=𝑚2(𝑇𝐶)(𝑡𝐹)𝑚1(𝑇)(𝑡𝐹)𝑚1(𝐶)(𝑡𝐹).(33) 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:𝑚2(𝐹𝑇)𝜕(𝑡𝐹)=2𝑔𝑧1,𝑧2,𝑧3,𝑡𝐹𝜕𝑧1𝜕𝑧2𝑧1=𝑧2=𝑧3=1,𝑚2(𝐹𝐶)𝜕(𝑡𝐹)=2𝑔𝑧1,𝑧2,𝑧3,𝑡𝐹𝜕𝑧1𝑧3𝑧1=𝑧2=𝑧3=1,𝑚2(𝑇𝐶)𝜕(𝑡𝐹)=2𝑔𝑧1,𝑧2,𝑧3,𝑡𝐹𝜕𝑧2𝑧3𝑧1=𝑧2=𝑧3=1.(34)

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 𝐴𝐹𝑇𝑡=𝑞2𝑚1(𝐹)𝑡𝑚𝐹1(𝑇)𝑡𝐹+𝑞1(𝑝)𝑞1(𝑑)𝑚1(𝐹)𝑡𝑚𝐹1(𝑇)𝑡𝐶+𝑚1(𝑇)𝑡𝑚𝐹1(𝐹)𝑡,𝐶(35) and similarly for 𝐴𝐹𝐶(𝑡) and 𝐴𝑇𝐶(𝑡) 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: 𝑚2(𝐹𝑇)(𝑡𝐹)=𝑄𝑇𝐹𝑡0𝑚1(𝑇)𝑡𝑡𝐴𝐹𝐹𝑇𝑡𝑚𝑑𝑡,2(𝐹𝐶)(𝑡𝐹)=𝑄𝑇𝐹𝑡0𝑚1(𝑇)𝑡𝑡𝐴𝐹𝐹𝐶𝑡𝑚𝑑𝑡,2(𝑇𝐶)(𝑡𝐹)=𝑄𝑇𝐹𝑡0𝑚1(𝑇)𝑡𝑡𝐴𝐹𝑇𝐶𝑡𝑑𝑡.(36)

By using the formulae (36), one can immediately calculate the covariances (33).

Figure 4 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.

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 𝑡=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 by 𝑄𝐹𝐷𝑑𝑡+𝑜(𝑑𝑡) the probability that a fast neutron is detected in the time interval (𝑡,𝑡+𝑑𝑡). Obviously, the detected neutron is also absorbed, that is, removed from the branching process. Let 𝐍(𝑡,𝑢) denote the number of fast neutrons detected in the time interval [𝑡𝑢,𝑡]. If 𝑡<𝑢, then 𝐍(𝑡) stands for the number of fast neutron detections in the time interval [0,𝑡].

4.1. Detection in Time Interval [0,𝑡]

Define the probabilities: 𝒫{𝐍(𝑡)=𝑛𝐹}=𝑝𝐹𝒫(𝑛,𝑡),{𝐍(𝑡)=𝑛𝑇}=𝑝𝑇(𝑛,𝑡),𝒫{𝐍(𝑡)=𝑛𝐶}=𝑝𝐶(𝑛,𝑡),(37) and introduce the generating functions: 𝐄𝑧𝐍(𝑡)=𝐹𝑛=0𝑝𝐹(𝑛,𝑡)𝑧𝑛=𝑔𝐹(𝐄𝑧𝑧,𝑡),𝐍(𝑡)=𝑇𝑛=0𝑝𝑇(𝑛,𝑡)𝑧𝑛=𝑔𝑇𝐄𝑧(𝑧,𝑡),𝐍(𝑡)=𝐶𝑛=0𝑝𝐶(𝑛,𝑡)𝑧𝑛=𝑔𝐶(𝑧,𝑡).(38) The equations determining the generating function can be easily written down. One obtains 𝑔𝐹(𝑧,𝑡)=𝑒𝑄𝐹𝑡+𝑄𝐹𝐴𝑡0𝑒𝑄𝐹(𝑡𝑡)𝑑𝑡+𝑧𝑄𝐹𝐷𝑡0𝑒𝑄𝐹(𝑡t)𝑑𝑡+𝑄𝐹𝑇𝑡0𝑒𝑄𝐹(𝑡𝑡)𝑔𝑇𝑔(𝑧,𝑡)𝑑𝑡,𝑇(𝑧,𝑡)=𝑒𝑄𝑇𝑡+𝑄𝑇𝐴𝑡0𝑒𝑄𝐹(𝑡𝑡)𝑑𝑡+𝑄𝐹𝑇𝑡0𝑒𝑄𝐹(𝑡𝑡)𝑞𝑝𝑔𝐹𝑧,𝑡𝑞𝑑𝑔𝐶𝑧,𝑡𝑑𝑡,(39)𝑔𝐶(𝑧,𝑡)=𝑒𝜆𝑡+𝜆𝑡0𝑒𝜆(𝑡𝑡)𝑔𝐹𝑧,𝑡𝑑𝑡.(40) Taking into account the second formula in (11), one has𝑞𝑑𝑔𝐶𝑧,𝑡=1𝑞1(𝑑)+𝑞1(𝑑)𝑔𝐶𝑧,𝑡.(41) The notation applied does not show that the delayed neutrons also participate in the process.

4.1.1. Expectations

The expectations of detected number of F type particles can be easily calculated by using the formulae: 𝐄{𝐍(𝑡)𝐹}=𝜕𝑔𝐹(𝑧,𝑡)𝜕𝑧𝑧=1=𝑛1(𝐹)(𝑡),𝐄{𝐍(𝑡)𝑇}=𝜕𝑔𝑇(𝑧,𝑡)𝜕𝑧𝑧=1=𝑛1(𝑇)(𝑡),𝐄{𝐍(𝑡)𝐶}=𝜕𝑔𝐶(𝑧,𝑡)𝜕𝑧𝑧=1=𝑛1(𝐶)(𝑡).(42)

After simple considerations one obtains the Laplace transforms of the expectations given by ̃𝑛1(𝐹)(𝑠)=𝑄𝐹𝐷𝑠+𝑄𝑇(𝑠+𝜆),𝑠𝒩(𝑠)̃𝑛1(𝑇)(𝑠)=𝑄𝐹𝐷𝑄𝐹𝑇𝑞1(𝑝)(𝑠+𝜆)+𝑞1(𝑑)𝜆,𝑠𝒩(𝑠)̃𝑛1(𝐶)(𝑠)=𝑄𝐹𝐷𝜆𝑠+𝑄𝑇,𝑠𝒩(𝑠)(43) where 𝒩(𝑠) is defined by (20). Clearly, these expressions are the Laplace transforms of the following integrals: 𝑛1(𝐹)(𝑡)=𝑄𝐹𝐷𝑡0𝑚1(𝐹)𝑡𝑛𝐹𝑑𝑡,1(𝑇)(𝑡)=𝑄𝐹𝐷𝑡0𝑚1(𝐹)𝑡𝑛𝑇𝑑𝑡,1(𝐶)(𝑡)=𝑄𝐹𝐷𝑡0𝑚1(𝐹)𝑡𝐶𝑑𝑡.(44) 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 𝑛1(𝐹)()=𝑄𝐹𝐷𝑄𝑇𝜆𝜔1𝜔2𝜔3,𝑛1(𝑇)()=𝑄𝐹𝐷𝜆𝑞1(𝑝)+𝑞1(𝑑)𝜔1𝜔2𝜔3𝑛1(𝐶)()=𝑄𝐹𝐷𝑄𝑇𝜆𝜔1𝜔2𝜔3,,(45) where𝜔1𝜔2𝜔3𝑄=𝜆𝐹𝑄𝑇𝑄𝐹𝑇𝑄𝑇𝐹𝑞1(𝑝)+𝑞1(𝑑).(46) 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,𝑡], 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 𝑛2(𝐹)(𝑡)=𝑄𝐹𝑇𝑡0𝑒𝑄𝐹(𝑡𝑡)𝑛2(𝑇)𝑡𝑑𝑡,𝑛2(𝑇)(𝑡)=𝑄𝑇𝐹𝑡0𝑒𝑄𝐹(𝑡𝑡)𝑞1(𝑝)𝑛2(𝐹)𝑡+𝑞1(𝑑)𝑛2(𝐶)(𝑡)𝑑𝑡+𝑄𝑇𝐹𝑡0𝑒𝑄𝐹(𝑡𝑡)×𝑞2𝑛1(𝐹)(𝑡)2+2𝑞1(𝑝)𝑞1(𝑑)𝑛1(𝐹)𝑡𝑛1(𝐶)𝑡𝑑𝑡,𝑛2(𝐶)(𝑡)=𝜆𝑡0𝑒𝜆(𝑡𝑡)𝑛2(𝐹)𝑡𝑑𝑡.(47) Introducing the notation𝐵𝐹𝐶𝑡=𝑞2𝑛1(𝐹)𝑡2+2𝑞1(𝑝)𝑞1(𝑑)𝑛1(𝐹)𝑡𝑛1(𝐶)𝑡,(48) and performing a Laplace transformation on (47), after the usual algebraic manipulations including the application of the convolution theorem, one obtains the solutions as 𝑛2(𝐹)(𝑡)=𝑄𝑇𝐹𝑡0𝑚1(𝑇)𝑡𝑡𝐵𝐹𝐹𝐶𝑡𝑑𝑡,𝑛(49)2(𝑇)(𝑡)=𝑄𝑇𝐹𝑡0𝑚1(𝑇)𝑡𝑡𝐵𝑇𝐹𝐶𝑡𝑑𝑡𝑛,(50)2(𝐶)(𝑡)=𝑄𝐹𝑇𝑡0𝑚1(𝑇)𝑡𝑡𝐵𝐶𝐹𝐶𝑡𝑑𝑡,(51)

With the help of these expressions, the variances of the number of fast neutron detections for the three starting particle types are given by 𝐃2{𝐍(𝑡)𝐹}=𝑛2(𝐹)(𝑡)+𝑛1(𝐹)(𝑡,𝐹)1𝑛1(𝐹),𝐃(𝑡,𝐹)2{𝐍(𝑡)𝑇}=𝑛2(𝑇)(𝑡)+𝑛1(𝐹)(𝑡,𝑇)1𝑛1(𝐹),𝐃(𝑡,𝑇)2{𝐍(𝑡)𝐶}=𝑛2(𝐶)(𝑡)+𝑛1(𝐹)(𝑡,𝐶)1𝑛1(𝐹).(𝑡,𝐶)(52)

4.2. Detection in the Time Interval [𝑡𝑢,𝑡],𝑡>𝑢

Since the detection in the time interval [0,𝑡𝑢) is excluded and only in the time interval [𝑡𝑢,𝑡] is permitted, one can write that 𝒫𝐍(𝑡,𝑢)=𝑛𝑆𝐹=𝗉𝐹=(𝑛,𝑡,𝑢)𝑗+𝑘+=𝑛𝑛1=0𝑛2=0𝑛3=0𝑝𝑛1,𝑛2,𝑛3,𝑡𝑢𝐹×𝑍(1)𝑗,𝑢𝑛1𝑍(2)𝑘,𝑢𝑛2𝑍(3),𝑢𝑛3,(53) provided that at the time moment 𝑡=0 one fast neutron was in the system. It is easy to prove that the probabilities 𝑍(1)(𝑗,𝑢𝑛1),𝑍(2)(𝑘,𝑢𝑛2), and 𝑍(3)(,𝑢𝑛3) are given by the following formulae: 𝑍(1)𝑗,𝑢𝑛1=𝑗1++𝑗𝑛1𝑛=𝑗1𝑖=1𝑝𝐹𝑗𝑖,𝑍,𝑢(2)𝑘,𝑢𝑛2=𝑘1++𝑘𝑛2𝑛=𝑘2𝑖=1𝑝𝑇𝑘𝑖,𝑍,𝑢(3),𝑢𝑛3=1++𝑛3𝑛=3𝑖=1𝑝𝐶𝑖,,𝑢(54) where 𝑝𝐹(𝑗𝑖,𝑢),𝑝𝑇(𝑘𝑖,𝑢), and 𝑝𝐶(𝑖,𝑢) are defined by (37). From (54) one can immediately see that 𝑗=0𝑍(1)𝑗,𝑢𝑛1𝑧𝑗=𝑔𝐹(𝑧,𝑢)𝑛1,𝑘=0𝑍(2)𝑘,𝑢𝑛2𝑧𝑘=𝑔𝑇(𝑧,𝑢)𝑛2,=0𝑍(3),𝑢𝑛3𝑧=𝑔𝐶(𝑧,𝑢)𝑛3.(55) Thus, the generating function𝗀𝐹(𝑧,𝑡,𝑢)=𝑛=0𝗉𝐹(𝑛,𝑡,𝑢)𝑧𝑛,𝑡>𝑢(56) can be written in the form:𝗀𝐹𝑔(𝑧,𝑡,𝑢)=𝑔𝐹(𝑧,𝑢),𝑔𝑇(𝑧,𝑢),𝑔𝐶(𝑧,𝑢),𝑡𝑢𝑆𝐹.(57) Since, as it is seen from(8)𝑔𝑧1,𝑧2,𝑧3,0𝐹=𝑧1,(58) from (57), one obtains𝗀𝐹𝑔(𝑧,𝑢,𝑢)=𝑔𝐹(𝑧,𝑢),𝑔𝑇(𝑧,𝑢),𝑔𝐶(𝑧,𝑢),0𝑆𝐹=𝑔𝐹(𝑧,𝑢),(59) which corresponds to the conditionlim𝑡𝑢𝗀𝐹(𝑧,𝑡,𝑢)=𝑔𝐹(𝑧,𝑢).(60) 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:𝗇1(𝐹)(𝑡,𝑢)=𝜕𝗀𝐹(𝑧,𝑡,𝑢)𝜕𝑧𝑧=1(61) of the number of fast neutron detections in the time interval [𝑡𝑢,𝑡], where 𝑡>𝑢, provided that the starting particle was also a fast neutron. By using (57) one obtains 𝗇1(𝐹)(𝑡,𝑢)=𝑚1(𝐹)(𝑡𝑢𝐹)𝑛1(𝐹)(𝑢)+𝑚1(𝑇)(𝑡𝑢𝐹)𝑛1(𝑇)(𝑢)+𝑚1(𝐶)(𝑡𝑢𝐹)𝑛1(𝐶)(𝑢),(62) where 𝑛1(𝐹)(𝑢),𝑛1(𝑇)(𝑢), and 𝑛1(𝐶)(𝑢) are defined by (44).

The calculation of the second factorial moment:𝗇2(𝐹)𝜕(𝑡,𝑢)=2g𝐹(𝑧,𝑡,𝑢)𝜕𝑧2𝑧=1(63) is straightforward, but rather tedious. One finds 𝑛2(𝐹)(𝑡,𝑢)=𝑚1(𝐹)(𝑡𝑢𝐹)𝑛2(𝐹)(𝑢)+𝑚1(𝑇)(𝑡𝑢𝐹)×𝑛2(𝑇)(𝑢)+𝑚1(𝐶)(𝑡𝑢𝐹)𝑛2(𝐶)(𝑢)+𝑚2(𝐹)𝑛×(𝑡𝑢𝐹)1(𝐹)(𝑢)2+𝑚2(𝑇)×𝑛(𝑡𝑢𝐹)1(𝑇)(𝑢)2+𝑚2(𝐶)(𝑛𝑡𝑢𝐹)1(𝐶)(𝑢)2+2𝑚2(𝐹𝑇)(𝑡𝑢𝐹)𝑛1(𝐹)(𝑢)𝑛1(𝑇)(𝑢)+2𝑚2(𝐹𝐶)×𝑛1(𝐹)(𝑢)𝑛1(𝐶)(𝑢)+2𝑚2(𝑇𝐶)(𝑡𝑢𝐹)𝑛1(𝑇)(𝑢)𝑛1(𝐶)(𝑢).(64) In this expression, the second factorial moments 𝑛2(𝐹)(𝑢),𝑛2(𝑇)(𝑢), and 𝑛2(𝐶)(𝑢) 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 instant 𝑡=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 parameter 𝑠𝐹. 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 probability 𝑃(𝑛1,𝑛2,𝑛3,𝑡;𝐹) of the event:𝐧1(𝑡)=𝑛1,𝐧2(𝑡)=𝑛2,𝐧3(𝑡)=𝑛3,(65) provided that at time 𝑡=0 there were no particles present in the system, is given by𝐺𝑧1,𝑧2,𝑧3𝑠,𝑡;𝐹=exp𝐹t0𝑔𝑧1,𝑧2,𝑧3,𝑡.𝐹1(66) In general, if the injected particles are of type 𝑖, and if the intensity of the injection of particle type 𝑖 is 𝑠𝑖, then one has𝐺𝑧1,𝑧2,𝑧3,𝑡;𝑆𝑖𝑠=exp𝑖𝑡0𝑔𝑧1,𝑧2,𝑧3,𝑡𝑆𝑖.1(67)

5.1.1. Calculation of the Expectations and the Variances

By using the logarithm of the generating function, one can write the formulae: 𝑀1(𝐹)𝑧(𝑡;𝐹)=𝜕ln𝐺1,𝑧2,𝑧3,𝑡;𝐹𝜕𝑧1𝑧1=𝑧2=𝑧3=1=𝑠𝐹𝑡0𝑚1(𝐹)𝑡𝑀𝐹𝑑𝑡,1(𝑇)𝑧(𝑡;𝐹)=𝜕ln𝐺1,𝑧2,𝑧3,𝑡;𝐹𝜕𝑧2𝑧1=𝑧2=𝑧3=1=𝑠𝐹𝑡0𝑚1(𝑇)𝑡𝑀𝐹𝑑𝑡,1(𝐶)𝑧(𝑡;𝐹)=𝜕ln𝐺1,𝑧2,𝑧3,𝑡;𝐹𝜕𝑧3𝑧1=𝑧2=𝑧3=1=𝑠𝐹𝑡0𝑚1(𝐶)𝑡𝐹𝑑𝑡,(68) giving the expectations of the numbers of fast and thermal neutrons and delayed neutron precursors, respectively, at the time moment 𝑡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 𝑀1(𝐹)(;𝐹)=𝑠𝐹𝑄𝑇𝜆𝒩𝑠𝑡,𝑀1(𝑇)(;𝐹)=𝑠𝐹𝑄𝐹𝑇𝜆𝒩𝑠𝑡,𝑀1(𝐶)(;𝐹)=𝑠𝐹𝑞1(𝑑)𝑄𝐹𝑇𝑄𝑇𝐹𝒩𝑠𝑡,(69) where𝒩𝑠𝑡=𝜔1𝜔2𝜔3𝑄=𝜆𝐹𝑄𝑇𝑞1(𝑝)+𝑞1(𝑑)𝑄𝐹𝑇𝑄𝑇𝐹.(70) The rest of the stationary expectations are given by the following formulae: 𝑀1(𝐹)(;𝑇)=𝑠𝑇𝑞1(𝑝)+𝑞1(𝑑)𝑄𝑇𝐹𝜆𝒩𝑠𝑡,𝑀1(𝑇)(;𝑇)=𝑠𝑇𝑄𝐹𝜆𝒩𝑠𝑡,𝑀1(𝐶)(;𝑇)=𝑠𝑇𝑞1(𝑑)𝑄𝐹𝑄𝑇𝐹𝒩𝑠𝑡,𝑀1(𝐹)(;𝐶)=𝑠𝐶𝑄𝑇𝜆𝒩𝑠𝑡,𝑀1(𝑇)(;𝐶)=𝑠𝐶𝑄𝐹𝑇𝜆𝒩𝑠𝑡,𝑀1(𝐶)(;𝐶)=𝑠𝐶𝑞1(𝑑)𝑄𝐹𝑄𝑇𝑞1(𝑝)𝑄𝐹𝑇𝑄𝑇𝐹𝒩𝑠𝑡.(71)

For the determination of variances of the number of fast neutrons, generated by the three different possible starting particles, one obtains the expressions: 𝑉𝐹𝑡;𝑆𝑖=𝜕2𝑧ln𝐺1,𝑧2,𝑧3,𝑡;𝑆𝑖𝜕𝑧21𝑧1=𝑧2=𝑧3=1+𝑧𝜕ln𝐺1,𝑧2,𝑧3,𝑡;𝑆𝑖𝜕𝑧1𝑧1=𝑧2=𝑧3=1=𝑠𝑖𝑡0𝑚2(𝐹)𝑡𝑆𝑖+𝑚1(𝐹)𝑡𝑆𝑖𝑑𝑡.(72)

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 parameter 𝑠𝐹, then applying the procedure described in [2], one can prove that generating function of the probability mass function:𝒫𝐍(𝑡,𝑢)=𝑛𝐧1=𝐧2=𝐧3=0=𝑃𝐹(𝑛,𝑡)(73) is given by 𝔾𝐹(𝑧,𝑡,𝑢)=𝔾𝐹(𝑠𝑧,𝑢)exp𝐹𝑡𝑢g𝐹𝑧,𝑡,𝑢1𝑑𝑡𝔾if𝑡>𝑢,𝐹𝑠(𝑧,𝑡)=exp𝐹𝑡0𝑔𝐹𝑧,𝑡1𝑑𝑡if𝑡𝑢,(74) whereg𝐹𝑧,𝑡𝑔,𝑢=𝑔𝐹(𝑧,𝑢),𝑔𝑇(𝑧,𝑢),𝑔𝐶(𝑧,𝑢),𝑡𝑢𝐹.(75)

By using the method described in [2, pages 85-86], one can prove that the improper integral:lim𝑡𝑡𝑢𝑔𝐹𝑧,𝑡,𝑢1𝑑𝑡(76) exists, consequently the stationary generating function:lim𝑡𝔾𝐹(𝑧,𝑡,𝑢)=𝔾𝐹(𝑠𝑡)(𝑧,𝑢)(77) also exists, if the three roots of the characteristic function (18) are nonnegative. In this case the random function 𝐍(𝑡,𝑢) converges in distribution to a random function 𝐍(𝑠𝑡)(𝑢), if 𝑡. Hence, 𝔾𝐹(𝑠𝑡)(𝑧,𝑢) is the generating function of 𝐍(𝑠𝑡)(𝑢), that is,𝔾𝐹(𝑠𝑡)𝑧(𝑧,𝑢)=𝐄𝐍(𝑠𝑡)(𝑢).(78) It is useful to rewrite the stationary generating function in the following form:ln𝔾𝐹(𝑠𝑡)(𝑧,𝑢)=𝑠𝐹𝑢0𝑔𝐹(𝑧,𝑡)1𝑑𝑡+𝑠𝐹𝑢g𝐹(𝑧,𝑡,𝑢)1𝑑𝑡.(79) Thus one findsln𝔾𝐹(𝑠𝑡)(𝑧,𝑢)=𝑠𝐹𝑢0𝑔𝐹(𝑧,𝑡)1𝑑𝑡+𝑠𝐹0𝑔𝑔𝐹(𝑧,𝑢),𝑔𝑇(𝑧,𝑢),𝑔𝐶(𝑧,𝑢),𝑡𝑆𝐹1𝑑𝑡.(80) If the injected particles are thermal neutrons, then the corresponding formula reads as ln𝔾𝑇(𝑠𝑡)(𝑧,𝑢)=𝑠𝑇𝑢0𝑔𝑇(𝑧,𝑡)1𝑑𝑡+𝑠𝑇0𝑔𝑔𝐹(𝑧,𝑢),𝑔𝑇(𝑧,𝑢),𝑔𝐶(𝑧,𝑢),𝑡𝑆𝑇1𝑑𝑡.(81) 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 interval 𝑢 is given by𝐹(𝑠𝑡)(𝑢)=𝜕ln𝔾𝐹(𝑠𝑡)(𝑧,𝑢)𝜕𝑧𝑧=1=𝑠𝐹𝑢0𝑛1(𝐹)𝑡𝑑𝑡+𝑠𝐹𝑢𝗇1(𝐹)𝑡,𝑢𝑑𝑡.(82) By using (79), one can write 𝑢𝗇1(𝐹)𝑡,𝑢𝑑𝑡=𝑛1(𝐹)(𝑢)0𝑚1(𝐹)(𝑡𝐹)𝑑𝑡+𝑛1(𝑇)(𝑢)0𝑚1(𝑇)(𝑡𝐹)𝑑𝑡+𝑛1(𝐶)(𝑢)0𝑚1(𝐶)(𝑡𝐹)𝑑𝑡,(83) where 0𝑚1(𝐹)(𝑡𝐹)𝑑𝑡=𝜆𝑄𝑇𝒩𝑠𝑡,0𝑚1(𝑇)(𝑡𝐹)𝑑𝑡=𝜆𝑄𝐹𝑇𝒩𝑠𝑡,0𝑚1(𝐶)(𝑡𝐹)𝑑𝑡=𝑞1(𝑑)𝑄𝐹𝑇𝑄𝑇𝐹𝒩𝑠𝑡.(84) After a lenghty algebra, one obtains𝐹(𝑠𝑡)(𝑢)=𝑠𝐹𝑄𝐹𝐷𝑄𝑇𝑄𝐹𝑄𝑇𝑄𝐹𝑇𝑄𝑇𝐹𝑞1(𝑝)+𝑞1(𝑑)𝑢,(85) which does not contain the delayed neutron precursor decay constant 𝜆.

By using the logarithmic generation function ln𝔾𝐹(𝑠𝑡)(𝑧,𝑢), one can write down the ratio of the variance to mean of the number of F types particles detected in the time interval 𝑢 in the form:𝐃2𝐍(𝑠𝑡)(𝑢)𝐄𝐍(𝑠𝑡)=𝕍(𝑢)𝐹(𝑠𝑡)(𝑢)𝐹(𝑠𝑡)1(𝑢)=1+𝐹(𝑠𝑡)(𝜕𝑢)2ln𝔾𝐹(𝑠𝑡)(𝑧,𝑢)𝜕𝑧2𝑧=1=1+𝑌(𝑢),(86) where 𝜕2ln𝔾𝐹(𝑠𝑡)(𝑧,𝑢)𝜕𝑧2𝑧=1=𝑠𝐹𝑢0𝑛2(𝐹)𝑡𝑑𝑡+𝑠𝐹𝑢𝗇2(𝐹)𝑡,𝑢𝑑𝑡,(87) and 𝑛2(𝐹)(𝑡,𝑢) is given by (64). The Feynman 𝑌(𝑢) function can be written in the traditional form as𝑌(𝑢)=3𝑖=1𝑌𝑖11𝑒𝜔𝑖𝑢𝜔𝑖𝑢.(88) The coefficients 𝑌𝑖 are rather involved functions of the roots 𝜔𝑖 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.

Figure 5 shows a quantitative illustration of the dependence of the variance to mean of the number of fast neutron detections on the detection time 𝑢. 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: 𝑄𝐹𝐷=1/10, 𝑄𝐹𝐴=7/30, 𝑄𝐹𝑇=2/3, 𝑄𝑇𝐴=14/10, 𝑄𝑇𝐹=3/5, 𝑞1(𝑝)=3, 𝑞1(𝑑)=0.02, 𝜆=0.1).

In view of the fact that there are three decay constants 𝜔𝑖,𝑖=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.

Acknowledgments

This work was supported by the Hungarian Academy of Sciences and the FP7 EU Collaborative Research Project FREYA, Grant Agreement no. FP7-269665.