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

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.


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][8][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 2 Science and Technology of Nuclear Installations emphasis was on the spatial effects as being the reason for the appearance of the multiple alpha modes [10][11][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 energydependent 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.

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 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 3 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 where Q FA and Q FT are the intensities of the absorption and thermalization of fast neutrons, respectively. Similar, selfobvious 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 f (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 that is, the numbers fast neutrons and delayed neutron precursors created in one reaction are independent. Further, let f d ( ) be given by where q (d) 1 ≤ 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 With this all quantities that are needed to formulate the problem are defined.

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 n F (t), n T (t), and n C (t), giving the numbers of particles of types F, T, and C, respectively, at the time moment t ≥ 0. Define the probabilities: P n F (t) = n 1 , n T (t) = n 2 , n C (t) = n 3 | S j = p n 1 , n 2 , n 3 , t | S j , (5) where the conditions S j , j = {1, 2, 3}, {S j } = {F, T, C}, indicating the type of particle starting the process, are defined as respectively. Introduce also the corresponding generating functions: 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 as well as where For later use, we note that Science and Technology of Nuclear Installations 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 with ν = ν p + ν d , and β = ν d /(ν p + ν d ).

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) m (F) 1 (t | S j ) of the number of fast neutrons, induced by one starting fast neutron, thermal neutron and delayed neutron precursor, respectively, one obtains the equations: Equations for the quantities m (T) 1 (t | S j ) and m (C) 1 (t | S j ) 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: then, from (14) one obtains After elementary algebra one obtains where is a third-order polynomial of s. It can be proven that the roots of the equation: are all real. Hence, introducing the notation 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 N (s) on s is shown (for the calculations, the following parameter values were used: The algebraic solutions for the Laplace transforms m (T) 1 (s | S j ) and m (C) 1 (s | S j ) can be obtained in a similar manner. Here we only list these solutions, which for m (T) 1 (s | S j ) are given by and for m (C) 1 (s | S j ) as 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, Science and Technology of Nuclear Installations of e −ω1t , e −ω2t and e −ω3t . As an illustration, we give the expectation m (F) 1 (t | F) of the number of fast neutrons, generated by one initial fast neutron injected into the system. One obtains 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: .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.

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 t ≥ 0, induced by an initial fast neutron, is given by the formula: Similar expressions can be derived for the other 8 variances.

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: one can derive the following equations: Analogous definitions can be introduced for m (T) 2 (t | S j ), and m (C) 2 (t | S j ) 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)   6 Science and Technology of Nuclear Installations Expectations F Expectations T Figure 3: Dependence of the expectations of the numbers of particles of types F, T, and C on the time, assuming that the starting particle was either F or T. and in the corresponding equations for m (T) 2 (t | S j ) and 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 Science and Technology of Nuclear Installations 7 and similar equations for m (T) 2 (s | S j ) and m (C) 2 (s | S j ). After simple algebra the solutions can be written in the following form: 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: 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 singleparticle 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 A F (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 functions e −ω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.

Covariances of the Numbers of Particles of Different Types.
Although not needed explicitly for the calculation of the twogroup 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 t = 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: 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.

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: 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 and similarly for A FC (t ) and A TC (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: 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 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.

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 Q FD dt + 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) Time Interval [0, t]. Define the probabilities:

Detection in
and introduce the generating functions: The equations determining the generating function can be easily written down. One obtains Taking into account the second formula in (11), one has The notation applied does not show that the delayed neutrons also participate in the process.

Expectations.
The expectations of detected number of F type particles can be easily calculated by using the formulae: Science and Technology of Nuclear Installations  Figure 4: Dependence of the covariances between the numbers of particles of types F 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.
After simple considerations one obtains the Laplace transforms of the expectations given by , where N (s) is defined by (20). Clearly, these expressions are the Laplace transforms of the following integrals: 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 where 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.

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 Introducing the notation and performing a Laplace transformation on (47), after the usual algebraic manipulations including the application of the convolution theorem, one obtains the solutions as With the help of these expressions, the variances of the number of fast neutron detections for the three starting particle types are given by

Detection in the Time Interval
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, un 1 ), Z (2) (k, un 2 ), and Z (3) ( , u | n 3 ) are given by the following formulae: where p F ( j i , u), p T (k i , u), and p C ( i , u) are defined by (37). From (54) one can immediately see that Thus, the generating function can be written in the form: Since, as it is seen from(8) from (57), one obtains which corresponds to the condition Expressions similar to (57) can be obtained for the cases when the starting particle is a fast neutron or a delayed neutron precursor.

Expectation and Second Factorial Moment.
For later use, let us determine the expectation: Science and Technology of Nuclear Installations 11 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 where n (F) 1 (u), n (T) 1 (u), and n (C) 1 (u) are defined by (44). The calculation of the second factorial moment: is straightforward, but rather tedious. One finds In this expression, the second factorial moments n (F) 2 (u), n (T) 2 (u), and n (C) 2 (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.

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 t = 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.

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 s F . 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 P(n 1 , n 2 , n 3 , t; F) of the event: provided that at time t = 0 there were no particles present in the system, is given by In general, if the injected particles are of type i, and if the intensity of the injection of particle type i is s i , then one has (67)

Calculation of the Expectations and the Variances.
By using the logarithm of the generating function, one can write the formulae: 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 where The rest of the stationary expectations are given by the following formulae: For the determination of variances of the number of fast neutrons, generated by the three different possible starting particles, one obtains the expressions:

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 s F , then applying the procedure described in [2], one can prove that generating function of the probability mass function: is given by where By using the method described in [2, pages 85-86], one can prove that the improper integral: exists, consequently the stationary generating function: 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 It is useful to rewrite the stationary generating function in the following form: Thus one finds If the injected particles are thermal neutrons, then the corresponding formula reads as In the continuation, we will only deal with the case of fast neutron injection. Variance-to-mean ratio

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 u is given by By using (79), one can write After a lenghty algebra, one obtains which does not contain the delayed neutron precursor decay constant λ.
By using the logarithmic generation function ln G (st) F (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: where and n (F) 2 (t , u) is given by (64). The Feynman Y (u) function can be written in the traditional form as The coefficients Y i 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. 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 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: Q FD = 1/10, Q FA = 7/30, Q FT = 2/3, Q TA = 14/10, Q TF = 3/5, q (p) 14 Science and Technology of Nuclear Installations 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) Feynmanalpha 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).

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.