REACTIVITY MEASUREMENTS IN ACCELERATOR DRIVEN SYSTEMS

Ch10-I045064.tex

27/6/2007

C H A P T E R

9: 42

Page 259

T E N

Reactivity Measurements in Accelerator Driven Systems

Contents 10.1 10.2 10.3 10.4

Steady Spallation Source Pulsed Poisson Source with Finite Pulse Width Pulsed Compound Poisson Source with Finite Width Periodic Instantaneous Pulses

260 264 281 283

The traditional methods of reactivity measurement with fluctuation analysis, discussed in the previous chapter, are all based on the use of a steady extraneous source with Poisson statistics. In the measurements a radioactive source is used, such as an Am–Be, Sb–Be or Pu–Be source with a constant intensity and simple Poisson statistics. The study of the statistical fluctuations in subcritical source driven systems became interesting again recently due to the interest in so-called accelerator driven systems, (ADS) [74–76, 83]. An ADS is a subcritical reactor, driven by a strong external source, which utilises e.g. the spallation reaction for neutron generation. The advantages of such a system are that it has excellent operational safety properties (a criticality accident is practically impossible), can utilise fertile nuclides as fuel such as 232 Th or 238 U much easier than traditional reactors (hereby having access to more abundant fuel reserves as the current reactors), and finally, on the long run and with proper design, it would produce less high-level nuclear waste. As a matter of fact, an ADS can incinerate more waste than it produces, so it can also be used primarily for transmutation of nuclear waste with energy production as a by-product. Due to the faster neutron spectrum and hence shorter neutron generation times, as well as the smaller fraction of delayed neutrons in most minor actinides whose burning is one purpose of the ADS, such a system has to be run at a relatively deep subcriticality, such as with a keff = 0.95 or lower. Hence to achieve sufficiently high power, a rather intensive extraneous source is needed. The dominating current ADS concept is based on a source utilising spallation. A spallation source will differ from a traditional one on at least two accounts. First, it emits a large, random number of neutrons on each spallation event, hence it will have a so-called compound Poisson statistics. Second, most likely it will run in a pulsed mode, i.e. it will be non-stationary. Both of these features deviate from those of the traditional sources. Hence the statistics of the neutron numbers and the detection events will also be different in subcritical systems driven with such a source. The branching processes generated in a subcritical system, driven with a spallation source will be investigated in this chapter. First, the Feynman- and Rossi-alpha formulae will be derived by assuming a steady source with compound Poisson statistics. Then the same formulae will be derived for periodic pulsed sources with a Poisson distribution of neutrons in the pulse. Both deterministic and random pulsing will be considered. Just as in the preceding chapter, the treatment will be confined to infinite homogeneous systems and detectors and to one energy group. This will provide some, albeit not perfect, transparency to the derivations and results, yet being suitable for comparison with experimental results and evaluation of measurements. Interested readers will find space- and energy- dependent treatments in [28, 64–66, 77, 78].

Neutron fluctuations ISBN-13: 978-0-08-045064-3

© 2008 Elsevier Ltd. All rights reserved.

259

Ch10-I045064.tex

27/6/2007

9: 42

Page 260

260

Imre Pázsit & Lénárd Pál

10.1 Steady Spallation Source In a spallation source, a large and random number of neutrons (typically several tens of neutrons) are generated by each incoming projectile, usually a high-energy proton. In a thick target, a whole shower of spallation reactions will take place. It will be assumed that the target is small enough such that all neutrons in one reaction can be assumed as being born simultaneously. Such sources will be called multiple emission sources. These were already described in Part I, Section 3.2, for the case without delayed neutrons. Again, it would be possible to extend the formalism of Chapter 3 to include delayed neutrons. However, similarly to the previous chapter, the traditional reactor physics derivation will be given here instead. It can be mentioned by passing that emission of multiple neutrons from a source takes place even with sources based on spontaneous fission, such as a 252 Cf source. Likewise, in a core containing spent fuel, there is an intrinsic source from the higher actinides that were generated during the previous fuel cycle and which contain isotopes with spontaneous fission, such as even mass number Pu nuclides. So the case of sources with compound Poisson statistics is interesting even outside the ADS area, and it has been already discussed in the literature [77,78]. However, the effects of the relatively low multiplicity of the source neutrons in case of a 252 Cf source or from other transuranic elements in spent fuel is negligible, hence this question was not given much attention before. As it will be seen here, the large neutron multiplicity of spallation neutrons (or even that of the neutrons in the pulse of a neutron generator) has a much more significant effect, which warrants a study of those effects in more detail. The temporal distribution of spallation events can be assumed exponential, leading to Poisson statistics of the number of spallation events in a time interval with a constant intensity S. The distribution of the number of particles emitted per event will be denoted by pq (n) satisfying the equation ∞


pq (n) = 1,

(10.1)

n=0

and its generating function as r(z) =

∞


pq (n)zn .

(10.2)

n=0

Only the first two moments of the above distribution will enter the relevant final formulae, and these are r1 = q =




npq (n)

and

r2 = q(q − 1) =

n




n(n − 1)pq (n).

(10.3)

n

For consistence of notation it is practical to introduce the Diven factor of the source: Dq =

q(q − 1) r2 = 2. 2 q r1

(10.4)

10.1.1 Feynman-alpha with a steady spallation source The assumptions of the model will also be the same as before: an infinite homogeneous system, infinite detector, one group of delayed neutrons, and one-group theory. Similarly to the derivation of the Rossi-alpha formula in the previous chapter, in the remaining part of Section 10.1 the solution will be given only by the backward formalism. The derivation used follows closely that of [79]. The derivation based on the forward equations is found in [80]. Since the difference compared to the traditional case lies entirely in the different statistics of the source, the distributions induced by a single particle will not be affected. Hence it is only the formula connecting the single-particle-induced and source-induced distributions which has to be re-calculated; the rest of the

Ch10-I045064.tex

27/6/2007

9: 42

Page 261

261

Reactivity Measurements in Accelerator Driven Systems

calculations will go along a line very similar to that given in Section 9.3. With the same arguments as with (9.59), one can write, in first order of dt: P(N , C, Z, T , t) = (1 − S dt)P(N , C, Z, T , t − dt)

P(N − n, C − c, Z − z, T , t) pq (k) Ak (n, c, z, T , t), + S dt n,c,z

(10.5)

k

where, analogously to (9.77), the function Ak (n, c, z, T , t) is defined as k 




Ak (n, c, z, T , t) =

p(nj , cj , zj , T , t).

(10.6)

n1 +···+nk =n j=1 c1 +···+ck =c z1 +···+zk =z

With the usual steps, the solution of the equation connecting the corresponding generating functions of the single-particle-induced distribution p


xn yc vz p(n, c, z, T , t) g(x, y, v, T , t) = n

c

z

and that of the source-induced distribution P,


xN yC vZ P(N , C, Z, T , t) G(x, y, v, T , t) = N

C

Z

which were introduced in (9.60) and (9.61), reads as   t  G(x, y, v, T , t) = exp S {r[g(x, y, v, T , t  )] − 1}dt  .

(10.7)

0

This is formally identical with the solution (3.24) of Section 3.2. The form of the relationship is not affected by the presence of delayed neutrons or detectors in the system; those only affect p(n, c, z, T , t), but not the source properties. It can also be shown that the asymptotic value, i.e. the limit of the generating function, limt→∞ G(x, y, v, T , t) exists if the system is subcritical. The various asymptotic factorial moments, corresponding to (9.65)–(9.70) are now easy to calculate. The first moments read as  ∞ r1 S
, (10.8) n(t)dt = N = r1 S −ρ 0  ∞ z(t, T )dt = λf NT . (10.9) Z(T ) = r1 S 0

The second moments, and notably the modified variances are obtained as  ∞  ∞ qnn (t)dt + r2 S n2 (t)dt µNN = N 0

and

 µZZ (T ) = 0

∞

(10.10)

0

 qzz (t, T )dt + r2 S

∞

z2 (t, T )dt,

(10.11)

0

where the second factorial moments qnn (t) and qzz (t, T ) of the single-particle-induced distributions are the same as before, given by (9.95) and (9.103). An inspection of (9.103) shows that the expression of qzz (t, T ) contains two different quadratic forms of z(t, T ), out of which one having the same form as the last term on

Ch10-I045064.tex

27/6/2007

9: 42

Page 262

262

Imre Pázsit & Lénárd Pál

the right-hand side of (10.11), which expresses the effect of the multiple source. Hence it can be envisaged that the dependence of µZZ on the measurement time T will be the same as in the traditional case of a simple Poisson source, but its amplitude enhanced with a factor, depending on the various parameters of the system. Since Z(T ) has also the same dependence on the measurement time as in the traditional case, the structure of the Feynman-alpha formula with a steady spallation source (or any multiple emission source, i.e. a 252 Cf source) will be the same as in the traditional case. Evaluating the integrals in (10.11) shows that this is indeed the case. Again, after performing the integrals, by re-denoting the measurement time T as t, the Feynman-alpha formula is now obtained as follows: λ2f µZZ (t) = [W (αp )fp (t) − W (αd )fd (t)], Z(t) (α + λ)(αp − αd )

(10.12)

where the parameters αp , αd and the functions fp (t) and fd (t) are the same as in (9.113) and (9.112), respectively, whereas the function W (s) is modified to   λ2 ν λ2 W (s) = 1 − 2 νp (νp − 1) + r2 (−ρ) − 2 2 νp νd . (10.13) s s r1 To see the effect of the extra term arising from the source multiplicity, it is worth re-writing the formula in the standard Feynman Y form µZZ (t) = Yp fp (t) + Yd fd (t) Z(t) and consider the prompt term Yp fp (t) ≡ Y1 f1 (t) only. Noting as before that αp  λ, the last term of (10.13) can be neglected, and one obtains

    Dνp Dνp r 1 Dq 1 − e −αt 1 − e −αt (1 + δ) 1 − 1+ (−ρ) 1− = . (10.14) Y1 f1 (t) = ν Dνp (ρ − β)2 αt (ρ − β)2 αt Here the parameter δ=

r1 Dq (−ρ) νDνp

(10.15)

was introduced, following Pázsit and Yamane [80]. In a similar manner, the effect of the source multiplicity can be quantified in the case of six delayed neutron groups with a generalisation of (9.120): ⎡ ⎤ 6 2 ν ν 
λ j j 0 ⎦. (10.16) Yi = 2λ2f ωi ⎣ν0 (ν0 − 1)(1 + δ) − 2 2 − λ2 s i j j=1 It is thus seen that for multiple emission sources, the variance to mean formula changes only very slightly. The time-dependence of the formula, which is the basis of the reactivity determination, is the same as with a simple Poisson source. Another similarity is that the source strength disappears from the expression for the variance to mean even with multiple emission sources. It is also obvious that the presence of the multiple emission is beneficial for the application of the Feynman method for the measurement of the reactivity, since it increases the amplitude of the useful part of the variance, i.e. the factor Y which stands for the deviation from the Poisson statistics in the count rate. The physical reason for the fact that the time-dependence does not change but the amplitude increases is simple. The time-dependence is determined by the rate by which the correlations between the numbers of the neutrons (or corresponding detector counts) born in the same chain, will die out. This depends only on the properties of the system, hence the presence of a steady multiple emission source will not change it. The amplitude of the variance to mean, on the other hand, depends also on the number of neutrons having a

Ch10-I045064.tex

27/6/2007

9: 42

Page 263

263

Reactivity Measurements in Accelerator Driven Systems

4.0 Am–Be Cf-252 Random–PNG Fitted curves

Y value

3.0

2.0

1.0

0.0 0.00

0.01

0.02 0.03 Gate width (s)

0.04

Figure 10.1 Measured Y values in an experiment with an Am–Be source, a neutron generator. keff = 0.9874 (from [81]).

0.05

252

Cf source and a randomly pulsed

common ancestor. In the traditional case all those neutrons are generated in the chain, started by the individual source particles. The number of such neutrons will increase when several neutrons enter the system from the source simultaneously, hence the increase of the amplitude of the time-dependent part of the variance to mean. The presence of the ‘enhancement factor’ δ > 0 is beneficial for the determination of the reactivity, which can be important for deeply subcritical systems. A quantitative estimate of the parameters determining the magnitude of δ yields that the Diven factors show a very little variation whether it regards fission or spallation, i.e. the factor Dq /Dνp will be in the order of unity. For a 252 Cf source in a core loaded with enriched uranium, the factor r1 /ν = q/ν is also in the order of unity, hence for a keff = 0.95, the correction represented by δ is below 10%. The situation is the same for an intrinsic source represented by spent fuel. For spallation sources, on the other hand, the situation is different. For spallation with protons in the GeV range and a thick Pb target, r1 ≈ 40 neutrons per spallation event, and accordingly, with keff = 0.95, δ ≈ 3. For larger subcriticalities the value increases further, and for keff = 0.7 one has δ ≈ 27. Hence for deeply subcritical systems the large source multiplicity becomes a very significant factor in measuring the reactivity of the system. The enhancing effect of the factor δ was confirmed in an experiment performed by Kitamura et al. [86]. In that experiment, a multiple emission source with a large δ value was achieved in form of narrow pulses of a D–T neutron generator which were triggered in a random fashion. The random triggering was induced by detecting gamma-rays from a 60 Co source, and using the detections as trigger signals. The average detection rate was approximately 108 pulses per second. In addition, an Am–Be, and a 252 Cf source were also used. The high δ value of the randomly pulsed neutron generator was attained by the large mean number of neutrons in the pulse. The results are shown in Fig. 10.1 for a measurement in a system with keff = 0.9874 (ρ ≈ −0.013). It is seen that the Feynman-Y curves for the multiple emission 252 Cf source are indistinguishable from those of the Am–Be source, agreeing with the numerical estimates of δ for a 252 Cf source and the level of subcriticality, as indicated above. However, the Y function corresponding to the true randomly pulsed neutron generator has a markedly larger amplitude than the other two curves, showing the effect of large source multiplicity.

10.1.2 Rossi-alpha with a steady spallation source Based on the foregoing, the derivation of the Rossi-alpha formula can be made quite simply. All one needs is the generalisation of the formula (9.145) to the case of multiple source emission, similarly as in (10.7). One obtains Gst (x1 , y1 , v1 ; x2 , y2 , v2 , τ) = lim G(x1 , y1 , v1 , t; x2 , y2 , v2 , t + τ) t→∞   ∞     = exp S {r[g(x1 , y1 , v1 , t ; x2 , y2 , v2 , t + τ)] − 1}dt . 0

(10.17)

Ch10-I045064.tex

27/6/2007

9: 42

Page 264

264

Imre Pázsit & Lénárd Pál

Again, like in the previous chapter, notation on the infinitesimal measuring time intervals dt at t and t + τ will be omitted. From this, one obtains for the covariance the equation  ∞  ∞ qzz (t, τ)dt + r2 S z(t)z(t + τ)dt, (10.18) CZZ (τ) = N 0

0

where qzz (t, τ) was defined before in (9.160). The effect of the multiple emission source is contained in the last term on the right-hand side. Evaluating the integral yields the results that, similarly to the case of the Feynman-alpha formula, the time-dependence of the Rossi-alpha remains the same with multiple emission sources as in the traditional case, whereas the amplitudes are modified the same way as in the case of the Feynman-alpha in the previous section. Concretely, one has λ2f dτ CZZ (τ) Pr (τ)dτ = = [W (αp )fp (τ) − W (αd )fd (τ)], Z dt 2(α + λ)(αp − αd )

(10.19)

where the functions W (s) are equal to those in (10.13) and the functions fi (τ) with (9.166). For the case of six groups of delayed neutron precursors can be made the same way as in the case of the Feynman-alpha formula. One obtains ⎡ ⎤ 6 6 6 2 ν ν 

λ 0 j 1
j ⎦ αi e −αi τ , Pr (τ) = Yi αi e −αi τ = λ2f ωi ⎣ν0 (ν0 − 1)(1 + δ) − 2 (10.20) 2 − λ2 2 i=0 s i j i=0 j=1 where the factors ωi are defined in (9.119).

10.2 Pulsed Poisson Source with Finite Pulse Width Based on the available technology of particle accelerators, it seems likely that a future ADS will be driven by a pulsed spallation source. Likewise, in several current model ADS experiments, a traditional pulsed neutron generator was used [83]. The case of the pulsed source is hence of practical interest, and in this section the Feynman- and Rossi-alpha formulae corresponding to a pulsed source we will be derived. In this section, delayed neutrons will be neglected in the derivation throughout. This is mostly motivated by the fact that in cores containing minor actinides, or even in traditional cores utilising fast fission, the delayed neutron fraction is much smaller than in reactors based on thermal fission of 235 U. The results with the inclusion of one group of delayed neutrons will be touched upon at the end of the section. At first we shall treat the case of a traditional neutron generator emitting finite width pulses with a simple Poisson distribution of neutrons in the pulse. Then the case of the finite width spallation source will be discussed. The derivations will mostly refer to the Feynman-alpha method. The derivation of the Rossi-alpha formula will be given without details, since the derivation can be inferred from the case of the Feynman formula. For very narrow pulses, it is a good approximation that all neutrons in the pulse are emitted simultaneously, as is indicated by the example in Fig. 10.1 [81]. The modelling of narrow pulses was already mentioned and discussed, at the level of the particle distribution, in Section 3. The case of narrow pulses can be described much simpler by assuming strongly periodic processes for the injection, i.e. periodic instantaneous injection, such as with the methods of Section 3. This methodology was applied by Degweker [60, 65] for the treatment of neutron fluctuations and for the pulsed Feynman-alpha method in accelerator driven systems. The case of periodic instantaneous injection will be treated at the end of this chapter. In the pulsed mode with finite width pulses, the distribution of the neutron injection can be considered as an inhomogeneous Poisson process, i.e. one with a time-dependent intensity S(t). Such cases were already treated in Chapter 3, however, those calculations concerned the distribution of the particle number, and not the number of detections in an interval. The extension of those methods to the distribution of the detection process would be relatively complicated. In the forthcoming a treatment will be used that builds on the backward

Ch10-I045064.tex

27/6/2007

9: 42

Page 265

265

Reactivity Measurements in Accelerator Driven Systems

equation based formalism of the preceding chapter, and by which both deterministic and random pulsing as well as arbitrary pulse shapes can be treated in a general framework [87]. As was seen in the preceding section and the previous chapter, the backward approach relies on the determination of the single-particle-induced distribution and the way how the single-particle and sourceinduced distributions are related to each other. Similarly to the case of multiple emission sources, the pulsed characteristics of the source will not alter the single-particle-induced distributions, hence again the difference as compared to the steady sources will lie at the expression that connects the single-particle-induced distribution with those due to an external source, switched on at t = 0. At the level of the generating functions, this formula with an inhomogeneous Poisson source reads as   t    (10.21) S(t )[g(x, v, T , t − t ) − 1]dt . G(x, v, T , t) = exp 0

Here the generating functions g(x, v, T , t) and G(x, v, T , t) are the same as those defined in (9.60) and (9.61), except that the auxiliary variable y, corresponding to the delayed neutrons, is missing. The formula needs to be considered in its asymptotic form when t → ∞. In the traditional case of constant source, this can be made in (10.21) by extending the upper limit of the integration to infinity. However, with a periodically stationary source, this limit has to be taken in a particular way which will be discussed soon below, hence the substitution is not denoted here.

10.2.1 Source properties and pulsing methods A pulsed neutron generator produces a periodic sequence (‘train’) of pulses.1 This means that the number of injected source neutrons, as well as the number of neutrons in the system and the number of detected neutrons, will not be stationary stochastic processes in the general sense, rather they will be periodically stationary. In addition to periodic stationarity, all moments will be oscillating quantities, in contrast to the smooth (non-oscillatory) behaviour of the moments in the case of a steady source. There exists a way of making the process stationary, although except the first moment, the moments cannot be made non-oscillatory. Namely, one can perform a Feynman- (or Rossi-) alpha experiment in two different ways: one can open the counting gate either synchronised with the pulsing, or randomly. The first case is referred to as ‘deterministic pulsing’, whereas the second as ‘stochastic pulsing’. In this latter case the process will be stationary. In reality of course it is not the pulsing of the generator which is random, rather the starting point of the measurements. In fact, since measurements are made nowadays with a high resolution time-to-digital converter (as opposed to a shift register), such that the individual detection times are registered together with the neutron pulse trigger, a given measurement can be evaluated both by the deterministic or stochastic pulsing method. The periodic pulsed character of the source is formulated mathematically as follows. Assume that the pulse repetition period is T0 , and that the pulse train consists of a periodic sum of the same pulse shape f (t) such that it is non-vanishing only within [0, T0 ), i.e.  ≥ 0 for 0 ≤ t ≤ T0 , (10.22) f (t) = 0 otherwise. Here f (t) is assumed to be normalised either by its integral or by the maximum value, whereas the total intensity of the source will be described by the help of an intensity factor S0 . Hence the time-dependent intensity of the source is given as S(t) = S0

∞
n=0

f (t − nT0 ) ≡ S0

∞


fn (t).

(10.23)

n=0

Various forms of the pulse shape f (t) will be specified later. For the discussion of the solution technique that follows below, f (t) can be left unspecified. This will exemplify one advantage of the technique used, namely 1 True

stochastic pulsing has also been achieved, see [81] and Fig. 10.1, but at high performance accelerators only the periodic pulsing is practical.

Ch10-I045064.tex

27/6/2007

9: 42

Page 266

266

Imre Pázsit & Lénárd Pál

that it can be applied to various pulse shapes with roughly the same effort, and that a substantial part of the steps does not depend on the pulse shape. This means that new pulse shapes can be treated with a rather moderate extra effort, once the problem was fully solved for one particular pulse shape. The case of random pulsing is then described by introducing a random variable ξ, taking non-negative real values ξ ∈ [0, T0 ) with a density p(ξ). Hence, P{ξ ≤ ξ ≤ ξ + dξ} ≡ p(ξ)dξ

(10.24)

is the probability that the random variable ξ takes a value between ξ and ξ + dξ. The random variable ξ describes the random starting point of the measurement within a pulse period T0 . The corresponding random source is then denoted as S(t|ξ = ξ) ≡ S(t|ξ) = S(t − ξ),

(10.25)

and with (10.25) one has S(t|ξ) = S0

∞


f (t − nT0 − ξ).

(10.26)

n=0

Such a source whose intensity itself is a random process is called a doubly random Poisson process, or Cox process. The deterministic and stochastic pulsing methods can actually be treated in a common framework by specifying the corresponding probability densities of the random variable ξ. For the stochastic pulsing, one will have  1/T0 for ξ ≤ T0 , p(ξ) = (10.27) 0 for ξ > T0 , whereas for the deterministic case it is given as p(ξ) = δ(ξ).

(10.28)

Expectations with respect to the random variable ξ will be denoted by brackets in this section in the definitions; wherever practical, bracket-free notations will be introduced for the first two factorial moments, as usual. Hence one has  T0 S(t|ξ) = S(t|ξ)p(ξ)dξ (10.29) 0

and S(t  |ξ)S(t  |ξ) ≡ K (t , t  ) =



T0

S(t  |ξ)S(t  |ξ)p(ξ)dξ.

(10.30)

0

It is easily confirmed that the covariance CovS (t  , t  ) = S(t  |ξ)S(t  |ξ) − S(t  |ξ)S(t  |ξ) = K (t  , t  ) − S(t  |ξ)S(t  |ξ)

(10.31)

is identically zero if and only if the density of the random variable ξ is a Dirac-delta function as in (10.28), i.e. for the case of a deterministic pulsing. For any other distribution p(ξ), CovS (t  , t  ) = 0, and in particular  t t CovS (t  , t  )dt  dt  > 0. (10.32) 0

0

This fact will have a bearing also on the variance to mean of the neutrons and the detections in a system driven with a randomly pulsed source. Along the same lines, the generating function
G(x, t|ξ) = P{N(t) = N |ξ = ξ}xN (10.33) N

Ch10-I045064.tex

27/6/2007

9: 42

Page 267

267

Reactivity Measurements in Accelerator Driven Systems

of the probability that N particles will be emitted by the source, i.e. by a doubly stochastic Poisson process during the time period [0, t), (which is the same as the number of particles at time t in an infinite system without absorption and multiplications), is given as    t S(t|ξ)dt . (10.34) G(x, t|ξ) = exp (x − 1) 0

The expected value of this generating function with respect to ξ is defined as     t  T0 G(x, t|ξ)p(ξ)dξ = exp (x − 1) S(t|ξ)dt . G(x, t|ξ) = 0

(10.35)

0

As an illustration let us consider now the case of stochastic pulsing with square pulses of width W < T0 , i.e. with f (t) = H (t) − H (t − W ) = (t, W ),

(10.36)

and hence S(t|ξ) = S0

∞


(t − nT0 − ξ, W ),

t ≥ 0.

(10.37)

n=−∞

For the first moment of the number of neutrons emitted by the random source between (0, t), from (32) one obtains     t  T0  1 S0 W d = S(t  |ξ)dξdt  = t, (10.38) G(x, t|ξ) N (t|ξ) ≡ N (t) = dx T0 0 0 T0 x=1 where, in line with the previous notations, the expectation sign was dropped. It is seen that the expectation N (t) is a smooth function of time, as expected. The second moment of the number of neutrons emitted by the stochastic source can be calculated from    2  t  t  T0  t t  d 1      N (t)(N (t) − 1) = = G(x, t|ξ) S(t |ξ)S(t |ξ)dt dt dξ = K (t  , t  )dt  dt  , dx2 T0 0 0 0 0 0 x=1 (10.39)   where the kernel K (t , t ) was introduced in (10.30). Then the variance of the number of neutrons in the system at time t is given from  t t 2 (t) = N (t)(N (t) − 1) − N 2 (t) + N (t) = N (t) + CovS (t  , t  )dt  dt  . (10.40) σN 0

0

By using the concrete form of the stochastic pulse train (10.37) in (10.31), the integrals in (10.40) can be easily evaluated by Laplace transform methods. One notices that K (t  , t  ) = K (|t  − t  |) ≡ K (u) is a symmetric function of its arguments. Hence,  t 0

t

K (t  , t  )dt  dt  = 2

0

 t 0

t

K (t  − t  )dt  dt  .

The inner integral of the above can be calculated by first performing a Laplace transform    t K˜ (s) K˜1 (s)    L K (t − t )dt = = , s s(1 − e −sT0 ) 0 where K˜1 (s) ≡



T0

e 0

−su

S2 K (u)du = 0 T0



(10.41)

0

 e −sW − 1 e −sT0 (e sW − 1) W (1 − e −sT0 ) . + + s2 s2 s

(10.42)

Ch10-I045064.tex

27/6/2007

9: 42

Page 268

268

Imre Pázsit & Lénárd Pál

The Laplace inverse of (10.42) is easily obtained by calculating the residues at s = 0 and s = ±2nπi/T0 ; n = 0, 1, 2 . . . (K˜1 (s) does not have a singularity at s = 0), which will in effect give the solution in form of a Fourier series. The remaining integral in (10.41) can be performed term by term on the series. The result, multiplied by 2, yields the second moment quantity N (t)(N (t) − 1). This finally leads to     ∞ 2 (t) σN nπt 2 2S0 T03
−4 nπW 2 sin =1+ sin n ≥ 1. N (t) Wtπ4 n=1 T0 T0

(10.43)

Equations (10.40) and (10.43) show the over-Poisson variance of the number of neutrons emitted by a stochastically pulsed source. In general, doubly stochastic processes have over-Poisson variances. The over-Poisson character can be explained by the positive integral of the covariance function, and hence, ultimately, by the fact that stochastic pulsing introduces correlations between the numbers of source particles at different times.

10.2.2 Calculation of the factorial moments for arbitrary pulse shapes and pulsing methods For the calculation of the factorial moments of the detector counts in a subcritical reactor with a subcritical source, either deterministically or stochastically pulsed, one has to use (10.21) generalised to the case of the pulsed source, denoted as in (10.25):  t (10.44) S(t  |ξ){g(x, v, t − t  ) − 1}dt  . G(x, v, t|ξ) = exp 0

The factorial moments have then to be taken from the expectation 

T0

G(x, v, t|ξ) =

G(x, v, t|ξ)p(ξ)dξ,

(10.45)

0

where, in the calculations of the factorial moments, first the derivatives with respect to x or z will be taken, after which the expectation w.r.t. ξ can be performed. The pulsing type is described by the probability density p(ξ) of ξ; furthermore, the pulse shape f (t) is only needed when the integrals w.r.t. ξ are performed. Thus, first the generic results will be derived, with the dependence on the realisation of the random variable ξ retained, from which then the deterministic and stochastic cases can be obtained for various pulse shapes. Again use can be made of the fact that the single-particle-induced distributions, i.e. the factorial moments of g(x, v, t) of (10.44) are known from the foregoing. Hence one can turn immediately to the calculation of the moments of G(x, v, t) of equation (10.45). The expressions for the first moments are as follows    t ∂G(x, v, t|ξ)  N (t) ≡ N (t|ξ) = = S(t  |ξ)n(t − t  )dt  ,  ∂x 0 x=v=1 

where

 N (t|ξ) =

t

S(t  |ξ)n(t − t  )dt  .

(10.46)

(10.47)

0

Likewise,  Z(t, T ) = Z(t, T |ξ) =

   t  t ∂G(x, v, t|ξ)     = S(t |ξ)z(t − t , T )dt = λ (t  , T )N (t − t  )dt  , d  ∂v 0 0 x=v=1 (10.48)

Ch10-I045064.tex

27/6/2007

9: 42

Page 269

269

Reactivity Measurements in Accelerator Driven Systems

with

 Z(t, T |ξ) =

t









S(t |ξ)z(t − t , T )dt = λd

0

t

(t  , T )N (t − t  |ξ)dt  .

(10.49)

0

The reason for the fact that the source-induced detector count can be written in two different ways is that it is a double convolution between three different functions. By using previous definitions and expressions, one out of two convolutions can be re-denoted as an autonomous function/notation, and this can be done by two different ways. This dual representation gives a possibility to choose the form that suits the actual computations better, which fact will be utilised below. The second factorial moment of the source-induced distribution    2 ∂ G(x, v, t|ξ)  (10.50) MZZ (t, T ) ≡ Z(Z − 1) =  ∂v2 x=v=1 can also be expressed in two equivalent ways, i.e.  t  t    2 S(t |ξ)mzz (t − t , T )dt + Z (t, T |ξ) = qzz (t  , T )N (t − t  , ξ)dt  + Z 2 (t, T |ξ), (10.51) MZZ (t, T |ξ) = 0

0

where mzz (t, T ) is the second factorial moment of the single-particle-induced neutron counts, defined in (9.54) and qzz is given in (9.103). For the Feynman-alpha formula one needs the modified variance, which can be written as µZZ (t, T ) = MZZ (t, T ) − Z 2 (t, T ). (10.52) From the above this leads to  µZZ (t, T ) = µZZ (t, T |ξ) =  =

t

S(t  |ξ)mzz (t − t  , T )dt  + Z 2 (t, T |ξ) − Z 2 (t, T )

0 t

qzz (t  , T )N (t − t  )dt  + Z 2 (t, T |ξ) − Z 2 (t, T ).

(10.53)

0

Here, it is obvious that due to the expectation with respect to the random variable ξ, the last two terms of (10.53) will be different. This is not the case for a steady source or for the deterministically pulsed case, where no such expectation has to be taken, and hence those two terms cancel each other. In that case, with neglecting notation on ξ, the convolution can be rearranged to  t  t    S(t )mzz (t − t )dt = qzz (t  , T )N (t − t  )dt  , (10.54) µZZ (t, T ) = 0

0

with qzz (t, T ) being the source term, given by   qzz (t, T ) = λf νp (νp − 1) z2 (t, T ),

(10.55)

which is obtained from (9.103) by neglecting the delayed neutrons. This is what is used in principle in the deterministic method.

10.2.3 General calculation of the variance to mean with arbitrary pulse shapes and pulsing methods As was seen in the previous section, the task is to calculate the kernels N (t|ξ), Z(t, T |ξ), and with their help the second order quantities MZZ (t, T |ξ) and/or µZZ (t, T |ξ). These need to be obtained for the case t → ∞ and then the expectation value of these quantities w.r.t. ξ needs to be taken. The formal solution for these can be given for arbitrary pulse shapes, with the pulse shape only affecting certain terms of a Fourier expansion. In this section the pulse will be described by its ‘mother function’ f (t), as defined in (10.22).

Ch10-I045064.tex

27/6/2007

9: 42

Page 270

270

Imre Pázsit & Lénárd Pál

The starting point is the Laplace transform of (10.47), which reads as ˜ N˜ (s|ξ) = S(s|ξ)˜ n(s).

(10.56)

By virtue of (10.22)–(10.26) and the fact that without delayed neutrons one has n(t) = e −αt , Equation (10.56) can be written as N˜ (s|ξ) =

˜ S0 e −sξ f˜ (s) S(s|ξ) = . s+α (s + α)(1 − e −sT0 )

Here, f˜ (s) is the Laplace transform of the truncated function f (t), which can be calculated as  T0  ∞ f (t)e −st dt = f (t)e −st dt, f˜ (s) =

(10.57)

(10.58)

0

0

whereas the term (1 − e −sT0 ) in the denominator arises from the summing up of the geometric series, given rise by the shifted integrals for the terms f (t − nT0 ). Since the function f˜ (s) does not have any singularities of its own, for times t > T0 , the inverse transform of (10.57) is entirely determined by the residues of the function e −sξ e st f˜ (s) . (s + α)(1 − e −sT0 )

(10.59)

The actual form of the pulse shape only affects the numerical values of the residues, but not the residue structure, which is as follows: 1. A pole at s = −α, which describes the transient after switching on the source at t = 0. 2. A pole at s = 0, the corresponding residue gives the asymptotic mean value of the oscillating function N (t). 3. An infinite number of complex conjugate roots on the imaginary axis at the values s = ±iωn , yielding harmonic functions in the time domain, representing a Fourier series expansion of the oscillating part of N (t) in the form sin (ωn t) and cos(ωn t). Here the notation ωn ≡

2nπ , T0

n = 1, 2, · · ·

(10.60)

was introduced. The term arising from the residue at s = −α will not give a contribution to the asymptotic value, but we keep it here for completeness. The residues at s = 0 and s = −α are given as Res s=0

f˜ (s = 0) e −sξ e st f˜ (s) = ≡ a0 , −sT (s + α)(1 − e 0 ) αT0

e −sξ e st f˜ (s) f˜ (s = −α) −α(t−ξ) ≡ a−α e −α(t−ξ) , = − αT e −sT 0 s=−α (s + α)(1 − e ) e 0 −1 Res

(10.61)

(10.62)

whereas Res{iωn } + Res{−iωn } can be written as Res

s=+iωn

e −sξ e st f˜ (s) e −sξ e st f˜ (s) + Res = an sin (ωn t) + bn cos (ωn t) cos (ωn ξ) −sT (s + α)(1 − e 0 ) s=−iωn (s + α)(1 − e −sT0 ) + bn sin (ωn t) − an cos (ωn t) sin (ωn ξ). (10.63)

Ch10-I045064.tex

27/6/2007

9: 42

Page 271

271

Reactivity Measurements in Accelerator Driven Systems

Here the parameters an and bn are defined as     ˜ ˜ 2 ωn f (iωn ) − α f (iωn ) an = T0 ωn2 + α2 and

(10.64)

    ˜ (iωn ) + ωn  f˜ (iωn ) α f 2 . bn = T0 ωn2 + α2

(10.65)

Then, N (t|ξ) is given as N (t|ξ) = Na (t|ξ) + S0 a−α e −α(t−ξ) = S0 [a−α e −α(t−ξ) + a0 ] + S0

∞


[an sin(ωn t) + bn cos(ωn t)] cos(ωn ξ)

n=1

+ S0

∞


[bn sin(ωn t) − an cos(ωn t)] sin(ωn ξ),

(10.66)

n=1

where Na (t|ξ) stands for the periodically asymptotic value of N (t|ξ). The actual value of the parameters a−α , a0 , an and bn will depend on the pulse form. Formally, however, all subsequent formulae, including the variance to mean can be calculated from the above expression. Equation (10.66) shows that in the stochastic case, where averaging w.r.t. ξ incurs integrations of the trigonometric functions over their period, all oscillating terms disappear in N (t), as is expected. In the deterministic case the oscillating term due to the cos(ωn ξ) remains present. For the calculation of the asymptotic value of the detector counts Z(T | ξ) = lim Z(t, T |ξ), t→∞

the last equality of (10.49) is the more suitable, by simply integrating (10.66) term by term. The execution of the limit t → ∞ requires, however, some attention. If the start of the measuring interval, t − T , always coincides with the beginning of a pulse, i.e. at time KT0 with some integer K , then of necessity one has t = KT0 + T , and it follows that letting t → ∞ equals to letting K → ∞. This procedure is not necessary for the stochastic pulsing, where the asymptotic number of neutrons becomes constant after the averaging w.r.t. ξ, but for the sake of uniform treatment, the limit will be calculated for both cases according to the above. With that, one has:  KT0 +T  KT0 +T  T (t  , T )N (t − t  |)dt  = lim λd N (t|ξ)dt = λd Na (t|ξ)dt. Z(T |ξ) = lim λd K →∞

K →∞

0

KT0

0

(10.67) The last step above results from the periodic character of the asymptotic form Na (t|ξ). Performing the integral will yield the result   ∞
Cn (T ) cos (ωn ξ) + Sn (T ) sin (ωn ξ) , (10.68) Z(T |ξ) = S0 λd a0 T + n=1

where Cn (T ) ≡

an {1 − cos(ωn T )} + bn sin(ωn T ) , ωn

(10.69)

Sn (T ) ≡

bn {1 − cos(ωn T )} − an sin(ωn T ) . ωn

(10.70)

Ch10-I045064.tex

27/6/2007

9: 42

Page 272

272

Imre Pázsit & Lénárd Pál

This is the general expression for arbitrary pulse shapes and pulsing techniques. The case of the deterministic pulsing is then obtained by substituting ξ = 0, whereas the stochastic pulsing is obtained by integrating (10.68) between zero and T0 w.r.t. ξ which leads to the disappearance of all oscillating parts. The calculation of the second order moments goes on similar lines. One finds it again that, due to the simple explicit expression for N (t|ξ), it is more convenient to use the last identity of (10.51) or (10.54). We shall use here the last equality of (10.51) where, from the concrete form of Z(t, T ) and (10.55) one has ⎧ 2 λ λf ν(ν − 1) ⎪ ⎪ (1 − e −αt )2 , 0 < t ≤ T, ⎨ d 2 α qzz (T , t) ≡ (10.71) 2 ⎪ ⎪ ⎩ λd λf ν(ν − 1) (e αT −1 )2 e −2αt , T < t. α2 As (10.71) shows, one needs to evaluate the integrals over trigonometric functions multiplied by the zeroth, first and second power of e −αt . These integrals can be jointly evaluated in a formal manner. The execution of the limit t → ∞ is performed the same way as in connection with (10.67). After lengthy but straightforward calculations, the result is given in the following form:    S0 λ2d λf ν(ν − 1) 1 − e −αT a T 1 − 0 α2 αT  ∞
−ωn an + 2αbn + an An (T ) + bn Bn (T ) + (1 − e −αT )2 cos(ωn ξ) ωn2 + (2α)2 n=1  ∞ 
−αT 2 ωn bn + 2αan + bn An (T ) − an Bn (T ) − (1 − e sin(ωn ξ) + Z 2 (T |ξ) − Z 2 (T ). ) 2 (2α)2 ω n n=1 (10.72)

µZZ (T |ξ) =

Here the following notations were introduced: An (T ) ≡ pn (0, T ) − 2pn (α, T ) + pn (2α, T ), with pn (α, T ) ≡ e

−αT



T

e αt sin(ωn t)dt =

0

α sin(ωn T ) + ωn {e −αT − cos(ωn T )} , ωn2 + α2

(10.73)

(10.74)

and Bn (T ) ≡ qn (0, T ) − 2qn (α, T ) + qn (2α, T ), with qn (α, T ) ≡ e −αT



T

e αt cos(ωn t)dt =

0

ωn sin(ωn T ) − α{e −αT − cos(ωn T )} . ωn2 + α2

(10.75)

(10.76)

Equations (10.68)–(10.70) and (10.72)–(10.76) serve as the formal solution for the Feynman-alpha formulae with arbitrary pulse shapes and arbitrary pulsing statistics. In the following they will be used to derive the concrete cases of square and Gaussian pulses and deterministic and stochastic pulsing methods.

10.2.4 Treatment of the pulse shapes and pulsing methods Square pulses The sequence of square pulses is described by f (t) = (t, W ),

(10.77)

Ch10-I045064.tex

27/6/2007

9: 42

Page 273

273

Reactivity Measurements in Accelerator Driven Systems

where the square window function was defined in (9.79). From this one has 1 − e −sW f˜ (s) = . s

(10.78)

Putting this into (10.61), (10.64) and (10.65) will yield the following results for a0 , an and bn : W , αT0

(10.79)

an =

ωn sin(ωn W ) + α{1 − cos(ωn W )} , nπ(ωn2 + α2 )

(10.80)

bn =

α sin(ωn W ) − ωn {1 − cos(ωn W )} . nπ(ωn2 + α2 )

(10.81)

a0 =

and

Gaussian pulses The reason for treating Gaussian-like pulse shapes is that there is experimental evidence that neutron generator pulses may have a bell-shaped form of time-dependence, which can be approximated, within the region they differ from zero, with a Gaussian fit. This was the case for instance with the GENEPI pulsed neutron generator, used in the European research program MUSE with the fast research reactor core MASURCA [83]. The pulse shape of the generator is shown in Fig. 10.2. Such a pulse can be approximated with a Gaussian, whose overwhelming area is embedded into a square window of width W if W ≈ 2 × FWHM (full width at half maximum). This is mostly for the purpose of easy comparison with the case of square pulses. Then one has   2 t − W2 f (t) = exp − (t − W/2, W ), (10.82) 2σ 2 where the window function (t, T ) is the same as before. The normalisation of the pulse shape f (t) was chosen such that the maximum value of the pulse function is unity, i.e. the same as with the square shape, just for practical purposes. With this, f˜ (s) is given as      ∞  W 2 2 √ −sW +(σs)2 − t − W2 − t − W2 −st −st ≈ = 2πσe 2 . e exp e exp (10.83) f˜ (s) = 2 2 2σ 2σ 0 −∞

Count rate (Arbitrary unit)

0.25 Measurement Fitting

0.2 0.15 FWHM

0.1 0.05 0

0

2

4

6

8

10

Time (µs)

Figure 10.2 Time shape of the neutron pulse of GENEPI, obtained by detection of associated alpha-particles of the T(d,n)4 He reaction in a silicon detector. The Gaussian fit gives a FWHM of 400 ns. (from [84]).

Ch10-I045064.tex

27/6/2007

274

9: 42

Page 274

Imre Pázsit & Lénárd Pál

By choosing σ = W /8, the cut-off values of f (t) become f (0) = f (W ) ≈ 3.4 · 10−4 , hence the error of the approximation in (10.83) is indeed negligible. Substituting (10.83) into (10.61), (10.64) and (10.65) yields √ 2πσ a0 = , (10.84) αT0   (ωn σ)2 √ exp − !ω  !ω " " 2 2πσ 2 n n an = cos α sin W + ω W , (10.85) n T0 (ωn2 + α2 ) 2 2   (ωn σ)2 √ exp −  !ω " " !ω 2 2πσ 2 n n bn = α cos W − ω W . (10.86) sin n T0 (ωn2 + α2 ) 2 2 The above two sets of parameters a0 , an and bn can now be used in the formulae with deterministic and stochastic pulsing, derived below.

Deterministic pulsing This case is equivalent with substituting ξ = 0 into the relevant formulae (10.68) and (10.72). All expressions containing the sin(ωn ξ) terms disappear, and likewise the last two terms of (10.72) cancel out. The results for the mean and the modified variance will become as follows (for simplicity, the expectation value sign, indicating the averaging w.r.t. ξ will be omitted):   ∞
an {1 − cos(ωn T )} + bn sin(ωn T ) , (10.87) Z(T ) = S0 λd a0 T + ωn n=1 and    S0 λ2d λf ν(ν − 1) 1 − e −αT µZZ (T ) = a0 T 1 − α2 αT   ∞ ∞

−ωn an + 2αbn −αT 2 an An (T ) + bn Bn (T ) + (1 − e + ) . ωn2 + (2α)2 n=1 n=1 (10.88) The particular cases of square and Gaussian pulses can be then obtained by substituting (10.79)–(10.81) or (10.84)–(10.86), respectively, into (10.87) and (10.88). As seen, the source strength S0 cancels out from the formula Y (T ) = µZ (T )/Z, similarly as in the case with continuous sources. Since both Z(T ) and µZ (T ) are asymptotically linearly dependent on the measurement time T , the Y (T ) curve goes into saturation, similarly as in the case with continuous sources. The formulae above also show that the oscillating deviations from the continuous Feynman-curve tend to zero asymptotically with the time gate T . The relative weight of the oscillations depends on the relationship between the pulse angular frequency ω = 2π/T0 and the prompt neutron time constant α. For the case ω  α, i.e. high pulse repetition frequency, the oscillations are small. Furthermore, with increasing pulse width W the relative magnitude of the oscillations decreases. For short pulses with low pulse frequency, ω ≤ α, the deviations from the smooth Y (T ) become quite significant. The magnitude of the oscillations depends rather weakly on the level of subcriticality, through the parameter α in the formulae, and it increases with increasing subcriticalities. The above qualitative analysis is illustrated quantitatively on Fig. 10.3, which shows deterministic Feynman Y (T ) functions for various pulse frequencies, pulse forms and pulse widths. As Fig. 10.3d shows, for low frequency pulsing, the pulse shape does have a certain influence on the fine structure of the Y (T ) function.

Ch10-I045064.tex

27/6/2007

9: 42

Page 275

275

Reactivity Measurements in Accelerator Driven Systems

W
0.001 ms, T0
2 ms

1

1

0.8

0.8 Y(T )

Y(T )

W
1 ms, T0
2 ms

0.6 0.4 Pulsed Steady

0.2 0

0

2

4

(a)

6

8

10

0.4

0

12

T (ms)

Y(T )

1

0.6 Gaussian Square Steady

0.4 0.2 4

(c)

6

4

8

10

6

8

12

0.6 Gaussian Square Steady

0.4 0.2 0

12

0

2

(d)

T (ms)

10

W
0.5 ms, T0
5 ms, s
0.125 ms

0.8

2

2

T (ms)

1

0

0

(b)

0.8

0

Pulsed Steady

0.2

W
0.5 ms, T0
2 ms, s
0.125 ms

Y(T )

0.6

4

6

8

10

12

T (ms)

Figure 10.3 Theoretical Feynman Y-curves for deterministic pulsing, with square and Gaussian pulses.

Stochastic pulsing In that case, as (10.27) shows, one needs to integrate (10.68) and (10.72) over the pulse period. Then all harmonic functions disappear in (10.68), and also the ones seen explicitly in (10.72). The only nontrivial case is the integration of the second last term in (10.72), Z 2 (T | ξ), which contains second order products of the sine and cosine functions, i.e. requires the evaluation of 1 T0



T0



0

∞


2 {Cn (T ) cos (ωn ξ) + Sn (T ) sin (ωn ξ)}

dξ.

(10.89)

n=1

However, substituting the definitions of Cn and Sn , and making use of the orthogonality relationships between the trigonometric functions, this integral is readily shown to be equal to ∞
C 2 (T ) + S 2 (T ) n

n=1

n

2

=2

∞ 2
a + b2 n

n=1

n

ωn2

sin2

!ω

n

2

" T .

(10.90)

Because the expected value of the detector counts is an especially simple function for the case of stochastic pulsing, i.e. Z(T ) = S0 λd a0 T , (10.91) it is practical to give directly the Y (T ) function. It will have the relatively simple form   λd λf ν(ν − 1) 1 − e −αT 1 − Y (T ) = α2 αT ∞ 2S0 λd
an2 + bn2 2 ! ωn " + T . sin a0 T n=1 ωn2 2

(10.92)

Ch10-I045064.tex

27/6/2007

9: 42

Page 276

276

Imre Pázsit & Lénárd Pál

α
300 s1

1

1

0.8

0.8

0.6 0.4

0.6 0.4

Stochastic Steady

0.2 0

α
400 s1

1.2

Y(T )

Y(T )

1.2

0

10

(a)

20 30 T (ms)

40

Stochastic Steady

0.2 0

50

0

10

(b)

α
500 s1 1

1

0.8

0.8

Y(T )

Y(T )

50

1.2

0.6 0.4

0.6 0.4

Stochastic Steady

0.2

(c)

40

α
600 s1

1.2

0

20 30 T (ms)

0

10

20

30

T (ms)

40

Stochastic Steady

0.2 0

50 (d)

0

10

20

30

40

50

T (ms)

Figure 10.4 Theoretical Feynman Y-curves for stochastic pulsing, with square pulses. W = 1 ms,T0 = 15 ms.

Again, the particular cases of square and Gaussian pulses are obtained by substituting (10.79)–(10.81) or (10.84)– (10.86), respectively, into (10.92). As the above shows, one significant difference between the deterministic and stochastic cases concerns the relative weight of the oscillating part to the smooth part. Most notably, in the stochastic pulsing case, the oscillating part is linear in the source strength, i.e. the source strength does not disappear from the variance to mean. This is a clear consequence of the ‘randomising’ of the pulse, which leads to a qualitatively different dispersion of the source neutrons, and in particular to the non-zero value of the auto-covariance of the source functions, as discussed earlier. This property is then transferred to the statistics of the neutron chain and that of the detector counts. It is interesting that the Diven factor of fission, which is present in all formulae of the continuous source and deterministic pulsing, is absent from the oscillating part of the Feynman formula of the stochastic pulsing. This can be expressed as if the oscillating part is controlled by the source statistics instead of the statistics of the multiplication of the fission chain. This also means that with a strong source, the relative oscillations become large. There are indications in recent experiments that this indeed occurs [85]. A second feature is that the smooth part of (10.92) is proportional to 1/α, whereas the oscillating part, through the factor 1/a0 , is proportional to α. This means that the relative weight of the oscillating part will increase faster with increasing subcriticalities than in the case of deterministic pulsing. Again, this can be said to be due to the fact that the oscillating part is influenced mostly by the source properties, whose significance increases in deep subcritical systems. Some representative calculated values of the stochastic pulsed Feynman-alpha curves are shown in Fig. 10.4, to illustrate the above qualitative analysis. The above formulae and their application to the evaluation of experiments was verified in measurements. Among others they were applied in the EU-supported project MUSE [66, 85]. Another dedicated series of measurements were made at the KUCA reactor (Kyoto University Critical Assembly) by using a pulsed neutron generator. The measurements were made with various subcriticality levels, by various pulse repetition frequencies (or periods), and pulse widths. Details of these measurements are found in [82, 86].The reference reactivities and corresponding alpha values were determined by the rod drop and the pulsed neutron technique, respectively. The data could be evaluated both by the deterministic and the stochastic method. The results of the deterministic method for the four different subcriticality levels are shown in Fig. 10.5, and those for the

Ch10-I045064.tex

27/6/2007

9: 42

Page 277

277

Reactivity Measurements in Accelerator Driven Systems

Deterministic method: 0.65$

Y value

3.0

α (Ref.)
266  2 s α (Fitted)
241  1 s1

2.0 1.0

Deterministic method: 1.30$ 2.5

1

2.0 Y value

4.0

0.0 0.00 0.01 0.02 0.03 0.04 0.05 (a) Gate widthT (s)

1.5 1.0 0.5

Measured Fitted curve

α (Ref.)
494  3 s1 α (Fitted)
500  2 s1

1.0 0.8 0.6 0.4 0.2

0.8

Measured Fitted curve

0.0 0.00 0.01 0.02 0.03 0.04 0.05 Gate widthT (s) (c)

Figure 10.5

Deterministic method: 2.72$ 1.0

Y value

Y value

1.2

Measured Fitted curve

0.0 0.00 0.01 0.02 0.03 0.04 0.05 (b) Gate widthT (s)

Deterministic method: 2.07$ 1.4

α (Ref.)
369  3 s1 α (Fitted)
355  1 s1

α (Ref.)
598  4 s1 α (Fitted)
607  4 s1

0.6 0.4 0.2

Measured Fitted curve

0.0 0.00 0.01 0.02 0.03 0.04 0.05 Gate widthT (s) (d)

Measured and fitted results for the deterministic pulsing (from [87]).

stochastic pulsing on Fig. 10.6. The experimentally determined Y values are shown as symbols together with the experimental errors, and the fitted curves are shown as continuous lines. The reference alpha values and the ones obtained from the curve fitting are shown in all four sub-figures. The fitting was made by assuming a square pulse form, based on experimental evidence with the neutron generator. The agreement between the experimental and the fitted curves is generally good. Especially for the stochastic pulsing the agreement between the fitted curves and alpha parameters with the measurements and the reference alpha values is excellent. There are larger deviations in the case of the deterministic method. As it will be seen, these deviations remain also when the measurement is modelled with periodic, instantaneous (infinitely narrow) pulses. This indicates that some assumptions of the model do not perfectly match the physical situation. One explanation was given in [66], showing that the simultaneous existence of an intrinsic, continuous source together with the pulsed source changes the shape of the Feynman-Y curve such that the relative amplitude of the oscillating part becomes smaller, as in the measurements above. In the stochastic pulsing method, the averaging procedure corresponding to the random distribution of the start of the measurement eliminates the effect of such unaccounted phenomena, which explains why this latter method gives better results. In addition, in the stochastic pulsing method an extra free parameter occurs, the source intensity, which multiplies only the oscillatory part, and hence adjusts the relative weight of the smooth and the oscillatory part, making it possible to compensate for the effect of an inherent source. No similar possibility exists in the deterministic pulsing method, rather one has to modify the formula by adding a traditional Feynman-Y curve. At any rate, the results here and other publications suggest, that the stochastic method appears to be more suitable for evaluating measurements than the deterministic pulsing.

10.2.5 Rossi-alpha with pulsed Poisson source Based on previous material, the Rossi-alpha formula for pulsed sources can be easily derived. As we have seen in the previous chapter, the difference between the Feynman- and Rossi-alpha formulae lies in the difference in the equation connecting the single-particle-induced and source-induced equations, so that instead of one-point distributions, one has to handle two-point ones. Then, as was seen in the previous section, the

Ch10-I045064.tex

27/6/2007

9: 42

Page 278

278

Imre Pázsit & Lénárd Pál

Stochastic method: 0.65$ 4.0

Stochastic method: 1.30$ 2.5

1

α (Ref.)
266  2 s

1

α (Fitted)
253  1 s

2.0 1.0 0.0 0.00 0.01 0.02 0.03 0.04 0.05 Gate width T (s)

1.5 1.0 0.5

Measured Fitted curve

(a)

Gate width T (s)

(b)

Stochastic method: 2.72$

α (Ref.)
494  3 s

α (Fitted)
495  3 s1

1

α (Fitted)
601  4 s

1.0 0.8 0.6 0.4

Measured Fitted curve

0.2

(c)

0.6 0.4

0.0 0.00 0.01 0.02 0.03 0.04 0.05 (d)

Gate width T (s)

Measured Fitted curve

0.2

0.0 0.00 0.01 0.02 0.03 0.04 0.05

Figure 10.6

α (Ref.)
598  4 s1

0.8 Y value

Y value

1.0

1

1.2

Measured Fitted curve

0.0 0.00 0.01 0.02 0.03 0.04 0.05

Stochastic method: 2.07$ 1.4

α (Fitted)
373  2 s1

2.0 Y value

Y value

3.0

α (Ref.)
369  3 s1

Gate width T (s)

Measured and fitted results for the stochastic pulsing (from [87]).

formula has to be modified in order to account for the time-dependence of the intensity of the source in the pulsed case. On the other hand, the single-particle-induced distributions do not change. In practice, one can therefore nearly ‘taylor together’ the solution from already known elements. For brevity, for the Rossi-alpha method only the case of the stochastic pulsing will be treated. Partly because experience shows that evaluating the measurements with the method of stochastic pulsing yields more consistent results than the deterministic pulsing; and partly because it is more similar to the ‘spirit’ of the original Rossi-alpha measurement where the counting gate is opened at a detection event in the first detector, which happens randomly in relation to the pulsing. Hence we start with an equation, corresponding to (10.17), but such that the time-dependent source intensity is accounted for, and simplified in the sense that the delayed neutrons are neglected: Gst (x1 , v1 ; x2 , v2 , τ|ξ) = lim G(x1 , v1 , t; x2 , v2 , t + τ|ξ) t→∞   t     = lim exp S(t |ξ)[g(x1 , v1 , t − t ; x2 , v2 , t + τ − t ) − 1]dt . (10.93) t→∞

0

The factorial moments that are sought have to be calculated from the expectation of the above generating function with respect to ξ, i.e.  Gst (x1 , v1 ; x2 , v2 , τ) = lim G(x1 , v1 , t; x2 , v2 , t + τ | ξ) = lim t→∞

t→∞ 0

T0

G(x1 , v1 , t; x2 , v2 , t + τ|ξ)p(ξ)dξ.

(10.94) Again, as in the previous cases of simple Poisson source and a compound Poisson source, the notations on the infinitesimal measurement times dt1 ≡ dt and dt2 ≡ dτ are neglected. The dependence of the distributions on the random variable ξ is denoted. The notations are otherwise standard and do not need explanations. The first moments regarding the number of neutrons N (t|ξ) or the detections Z(t, dt|ξ) at t or t + τ are the same as in the previous section, with straightforward alterations of the detection time from T to dt. The

Ch10-I045064.tex

27/6/2007

9: 42

Page 279

279

Reactivity Measurements in Accelerator Driven Systems

covariance function of the detector counts in the infinitesimal intervals dt at t and dτ at t + τ can be calculated from (10.93). Similarly to the modified variance of (10.53), one finds that due to the time-dependent stochastic source, the covariance will contain additional terms as compared to the steady source. The result is  t qzz (t, τ)N (t − t  |ξ)dt  CZZ (t, τ) = CZZ (t, τ|ξ)dξ = 0

+ Z(t, dt|ξ)Z(t + τ, dτ|ξ) − Z(t, dt|ξ)Z(t + τ, dτ|ξ).

(10.95)

Here, N (t|ξ) is given in (10.66), the single-particle-induced second factorial moment source term qzz (t, τ) is the same as in (9.160), and Z(t, dt|ξ) reads as Z(t, dt|ξ) =

λd S0 f˜ (s = −α) e −α(t−ξ) λd S0 f˜ (s = 0) dt − dt αT0 α e αT0 − 1 + λ d S0

∞


{Cn (t, dt) cos (ωn ξ) + Sn (t, dt) sin (ωn ξ)} ,

(10.96)

n=1

where now, due to the infinitesimal length of the measurement time, the coefficients Cn (t, dt) and Sn (t, dt) are given as Cn (t, dt) = [an sin(ωn t) + bn cos(ωn t)]dt, n = 1, 2, . . . (10.97) and Sn (t, dt) = [bn sin(ωn t) − an cos(ωn t)]dt, n = 1, 2, . . . , (10.98) whereas the coefficients an and bn are the same as in (10.64) and (10.65). As mentioned there, they depend on the actual pulse shape. The calculation goes along the same lines as in the previous section. The only term whose calculation requires some effort is the second last term of (10.95). However, due to the orthogonality of the trigonometric functions over the period T0 , over which the averaging in ξ is to be performed, this goes without problems. The final result, corresponding to the asymptotic case, will be written here in the alternative form of the Rossi-alpha formula, based on the function R(τ), as defined in (9.163) for the traditional case. For the pulsed source case treated here it is expressed as R(τ)dτ ≡ lim

t→∞

Z(t, dt|ξ)Z(t + τ, dτ|ξ) CZZ (τ) = + Zdτ. Z(t, dt|ξ) Zdt

A short calculation yields R(τ)dτ =

λd λf ν(ν − 1) −ατ λs S0 f˜ (s = 0) e dτ + dτ 2α αT0 ∞

+

λd S0 αT0
2 (an + bn2 ) cos(ωn τ)dτ. 2f˜ (s = 0)

(10.99)

n=1

Equation (10.99) shows that, similarly to the Feynman-alpha formula for randomly pulsed sources, the pulsed Rossi-alpha formula consists of two terms equivalent to (with an amplitude scaling) the Rossi-alpha formula of the traditional case, and an additional non-negative oscillating term. This last term shows periodical, undamped oscillations with the pulse period T0 . This is a definite difference compared to the pulsed Feynmanalpha formula, in which the oscillating part is not periodic, rather the oscillations have a decaying amplitude. The periodic character of the oscillatory part in the pulsed Rossi-alpha formula makes it very suitable for evaluating measurements. Consideration of the various pulse shapes goes exactly as in the previous sections, and the corresponding coefficients an and bn have already been calculated. For square pulses, i.e. f (t) = (t, W ),

Ch10-I045064.tex

27/6/2007

9: 42

Page 280

280

Imre Pázsit & Lénárd Pál

25 000 Equation (10.101)

R(τ) (s1)

20 000

First two terms of Equation (10.101)

15 000 10 000 5000 0 0.000

0.002

0.004 0.006 Time interval (s)

0.008

0.010

Figure 10.7 An example of a stochastically pulsed Rossi-alpha curve for square pulses,T0 = 2 ms,W = 10−3 ms.

these are given by (10.80) and (10.81), yielding an2 + bn2 =

! " 4 2 ωn W , sin n2 π2 (α2 + ωn2 ) 2

n = 1, 2, . . . .

(10.100)

Hence the Rossi-alpha formula for square pulses becomes R(τ) =

λd λf ν(ν − 1) −ατ λd S0 W + e 2α αT0

∞ ! " 2λd S0 αT0
1 2 ωn + sin W cos(ωn τ). π2 W n=1 n2 (α2 + ωn2 ) 2

(10.101)

For Gaussian pulse shapes, an and bn are given by (10.85) and (10.86), thus an2 + bn2 = and hence

8πσ 2 2 2 e −ωn σ , + ωn2 )

T02 (α2

n = 1, 2, . . . ,

√ λd λf ν(ν − 1) −ατ λd S0 2πσ R(τ) = e dτ + 2α αT0 √ ∞
1 2λd S0 α 2πσ 2 2 + e −ωn σ cos(ωn τ) . 2 2 T0 α + ωn n=1

(10.102)

(10.103)

Figure 10.7 shows a plot of the pulsed Rossi-alpha formula R(τ) for square pulses, i.e. (10.101). The parameters used correspond to a typical thermal reactor system [86]. The curve demonstrates the periodic oscillations around the traditional Rossi-alpha term, i.e. the first two terms of (10.101). The fact that the pulsed Rossi-alpha formula consists of a traditional Rossi-alpha expression and a periodic part, can be utilised in an interesting way. Normally, one would try to fit a curve to the expression (10.101). This is, however, a somewhat complicated task, in view of the fact that (10.101) contains a series expansion in which several terms need to be kept for good accuracy, since for narrow pulses the Rossi-alpha curve has sharp edges, as is seen in Fig. 10.7. This is a definite difference compared to the stochastically pulsed Feynman method, which remains smooth (edge-free) even for very narrow pulses. However, one can avoid fitting the highly oscillating terms by determining them from the experiment itself, taking the oscillatory part from the section of the covariance curve for large delay times t.

Ch10-I045064.tex

27/6/2007

9: 42

Page 281

281

Reactivity Measurements in Accelerator Driven Systems

0.08

R(τ) (100 µs1)

Measured 0.06

0.04

0.02 0.00

0.02

0.04 0.06 Time interval (s)

0.08

0.10

Figure 10.8 Result of a pulsed Rossi-alpha experiment (−ρ = 0.65$,T0 = 20 ms). From Kitamura et al. [86]. 0.05 Measured Fitted curve

R(τ) (100 µs1)

0.04 0.03 0.02 0.01 0.00 0.01 0.00

0.02

0.04

0.06

0.08

Time interval (s)

Figure 10.9 The non-oscillating component of Rossi-alpha curve (−ρ = 0.65$, T0 = 20 ms). From Kitamura et al. [86].

This has been demonstrated in an experiment performed in the KUCA using a D–T pulsed neutron source [86]. Figure 10.8 shows an experimental result of the pulsed Rossi-alpha method. In Fig. 10.9, the traditional Rossi-alpha curve is shown, which was obtained by subtracting the section of the experimental curve ranging from 0.08 s to 0.10 s from the curve ranging from 0.00 s to 0.08 s in Fig. 10.8. This part of the curve contains only the periodic term. The result of a single exponential fitting is also shown in Fig. 10.9. The resulting decay constant, i.e. the α ≈ 263 ± 1 s−1 showed a good agreement with that obtained by the pulsed neutron source method performed by using the same core configuration. More examples are given in [82, 86]. Further applications of the pulsed Rossi-alpha method in experiments are found in [85].

10.3 Pulsed Compound Poisson Source with Finite Width The above results can be extended to more complicated models with the same formalism. One extension is to assume that the statistics of the source neutron injection within the pulse follows a compound Poisson distribution. Such a case would correspond to a spallation based accelerator with a finite width of proton pulses. Another extension is to generalise the description to the inclusion of delayed neutrons in the branching process. The first of these two extensions will be described briefly below. The case of delayed neutrons is not considered here; such results can be found in [88].

Ch10-I045064.tex

27/6/2007

9: 42

Page 282

282

Imre Pázsit & Lénárd Pál

The extension of the results to a compound inhomogeneous Poisson source is straightforward, based on the treatment of the previous cases. From combining (10.7) and (10.44), one has  G(x, v, t|ξ) = exp

t





S(t |ξ){r[g(x, v, t − t )] − 1}dt



,

(10.104)

0

# where as before, r(z) = n pq (n) zn is the generating function of the source probability distribution pq (n). The calculations can be performed with basically the same techniques as in the foregoing, although they are more involved. Here we only quote the final results,2 using the same notations as before.

Feynman-alpha, deterministic pulsing The Feynman Y -term is obtained as   f˜ (s = 0) 1 − e −αT Y (T ) =Y0 T (1 + δ) 1 − αT0 αT    ∞
b n ωn + Y0 δ An (T ) an 1 + 1 − an α n=1 + Y0

∞
n=1



bn



a n ωn 1+ 1+ bn α

 δ Bn (T )

   1 ωn2 (1 −ω 1 + 1 + a − δ) + 2αb δ ∞ n n n
2 α2 −αT 2 (1 − e ) , + Y0 (2α)2 + ωn2 n=1

(10.105)

where Y0 =

λd λf ν(ν − 1) λd r1 S0 . α2 Z(T )

Here the parameter δ, defined in (10.15), δ=

r 1 Dq (−ρ) νDνp

appears again. Just as in the case of the steady spallation source, the part of the variance to the mean that exceeds unity is amplified by the appearance of the factor δ in the first term on the right hand side. This term contains a time-dependent factor which is the same as that of the traditional method. However, the amplitude of the oscillating part is also amplified, through the occurrence of the factor r1 = q. Hence, what regards the applicability of the pulsed Feynman-alpha for the evaluation of measurements, the presence of the multiple source may be slightly less advantageous, due to the amplification of the oscillating terms. In addition, since the factor δ is proportional to the absolute value of the reactivity, for systems near to critical, the oscillating terms are amplified more than the smooth part. 2 These

calculations were performed by Y. Kitamura (unpublished work).

Ch10-I045064.tex

27/6/2007

9: 42

Page 283

283

Reactivity Measurements in Accelerator Driven Systems

Feynman-alpha, stochastic pulsing As expected, the Feynman-alpha formula for the stochastic pulsing takes a much simpler form even in the case of a multiple emission source. The result is   λd λf ν(ν − 1) 1 − e −αT (1 + δ) 1 − α2 αT ∞ ! λd r1 S0 αT0
an2 + bn2 ωn " +2 sin T . 2 f˜ (s = 0)T n=1 ωn2

Y (T ) =

(10.106)

Here, the first term on the right-hand side, corresponding to the traditional formula, is modified exactly the same way as in the case of the steady spallation source. Similarly to the case of deterministic pulsing, the oscillatory part is amplified too, with the factor r1 .Thus for large source multiplicities r1 and slight subcriticalities |ρ|, the oscillatory parts will be amplified more by the source multiplicity than the smooth part.

Rossi-alpha, random pulsing As is quite predictable at this point, the formula (10.99) gets modified into Pr (τ)dτ =

λd λf ν(ν − 1) (1 + δ)e −ατ dτ 2α ∞ λd r1 S0 αT0
2 (an + bn2 ) cos(ωn τ)dτ. + 2f˜ (s = 0) n=1

(10.107)

The same comments are valid regarding the effect of the source multiplicity as in the case of the Feynman-alpha with random pulsing.

10.4 Periodic Instantaneous Pulses In this section the case of infinitely narrow pulses, i.e. instantaneous injection, will be treated.The Feynmanalpha formula will be calculated for both deterministic pulsing and random pulsing. For the Rossi-alpha, due to the nature of the measurement (the measurement gate is opened at random through detection of a neutron in the system), only the random pulsing case is interesting. These cases will be treated after an elegant and effective formalism used by Degweker in two publications [60, 65]. The formalism is basically the same as the one developed in Section 3.2 for the distribution of the number of particles, but extended to the case of distribution of the detections.

10.4.1 Feynman-alpha with deterministic pulsing It is supposed that pulses are emitted into the system at times t = nT0 , with n running through the integers from −∞ to ∞. The measurement starts at t = +0, and lasts over a gate length T , i.e. its start coincides with the arrival of a pulse (hence this case is also called ‘pulse-triggered Feynman-alpha measurement’ [65]). Since the pulsing was started at t = −∞, the system is already in a periodically stationary state at t = 0. It is supposed that in one pulse a random number of neutrons are injected with a probability distribution pq (n) and corresponding generating function r(z), i.e the same notations are used as before. The following quantities are defined: pd (n, t, T )

(10.108)

Ch10-I045064.tex

27/6/2007

9: 42

Page 284

284

Imre Pázsit & Lénárd Pál

is the probability that one single neutron, emitted into the system at t ≤ T will lead to n detections during the interval [0, T ), and ∞
zn pd (n, t, T ) (10.109) gd (z, t, T ) = n=0

is its generating function. Likewise, Pd (n, T ) is the probability that the pulse train together will lead to n detections during [0, T ), and Gd (z, T ) =

∞


zn Pd (n, T )

(10.110)

(10.111)

n=0

is the corresponding generating function. For the calculations, the first two factorial moments of gd (z, t, T ) will be needed. These have already been determined in Chapter 4. Due to the independence of the progeny production of the neutrons partly for all neutrons injected within one pulse, and partly for the different pulses, for the generating function Gd (z, T ) one can write down directly the following equation: Gd (z, T ) =

∞ 



[T /T0 ]

r[gd (z, −nT0 , T )]

n=0

r[gd (z, kT0 , T )],

(10.112)

k=1

where [T /T0 ] stands for the largest integer not exceeding T /T0 . Here the fact was accounted for that the measurement starts at t = +0, so the pulse arriving at t = 0 does not fall into the measurement period. The above expression could of course have been written in one single product with n running from −∞ to [T /T0 ]. However, the factorisation is useful for the calculations, because the one-particle-induced detection generating function gd (z, t, T ), and hence also its factorial moments, have different functional forms depending on whether t is negative or positive (i.e. if the source particles from the pulse arrived inside or outside the measuring period). To keep track of this difference, we shall use the subscripts ‘1’ and ‘2’ for the generating function and its factorial moments for the cases t ≤ 0 and t ≥ 0, respectively, i.e.  gd1 (1, t, T ), t ≤ 0, gd (1, t, T ) = gd2 (1, t, T ), t ≥ 0. For the calculation of the first moment M1 (T ) = Gd (1, T ) one obtains from (10.112)

⎡

Gd (1, T ) = r1 ⎣

∞
n=0

 gd1 (1, −nT0 , T ) +




[T /T0 ]

(10.113) ⎤  gd2 (1, kT0 , T )⎦ ,

(10.114)

k=1

where r1 is the expectation of the number of neutrons emitted in a pulse. To evaluate (10.114), one needs the expectation gd (1, t, T ). This was given in (4.24), and by a suitable redefinition of the arguments, it is obtained from  λd αt −αT ) = g  (1, t, T ), t ≤ 0, d1 α e (1 − e  gd (1, t, T ) = (10.115) λd −α(T −t) ) = g  (1, t, T ), t ≥ 0. d2 α (1 − e Substituting (10.115) into (10.114) and performing the summations leads to the result r1 λd [1 + [T /T0 ](1 − e −αT0 ) − e −α(T −[T /T0 ]T0 ) ] α(1 − e −αT0 ) r 1 λd = [1 + [T /T0 ](1 − e −αT0 ) − e −αu ], α(1 − e −αT0 )

M1 (T ) =

(10.116)

Ch10-I045064.tex

27/6/2007

9: 42

Page 285

285

Reactivity Measurements in Accelerator Driven Systems

where the notation u = T − [T /T0 ]T0

(10.117)

was introduced, in analogy with its definition and use in Chapter 3. Obviously, 0 ≤ u ≤ T0 is a running parameter, describing the variation of the measurement time between the kth and (k + 1)th period, such that T = kT0 + u, with k = [T /T0 ]. The second moment M2 (T ) of Gd (z, T ) is calculated from (10.112) with twofold derivation. A brief calculations yields

M2 (T ) =

∞


2  [r2 gd1 (1, −nT0 , T ) + r1 gd1 (1, −nT0 , T )]

n=0

+ r12

∞


 gd1 (1, −nT0 , T )

+




 gd1 (1, −n T0 , T )

n =n

n=0 [T /T0 ]

∞


2  [r2 gd2 (1, −nT0 , T ) + r1 gd2 (1, −nT0 , T )]

n=0




[T /T0 ]

+ r12

 gd2 (1, −nT0 , T )

∞


 gd2 (1, −n T0 , T ).

(10.118)

n =n

n=1

For the Feynman-alpha one actually needs the modified variance M2 (T ) − M12 (T ). Adding and subtracting the terms n = n and subtracting M12 (T ) leads to

M2 (T ) − M12 (T ) =

∞


2  [(r2 − r12 )gd1 (1, −nT0 , T ) + r1 gd1 (1, −nT0 , T )]

n=0




[T /T0 ]

+

2  [(r2 − r12 )gd2 (1, nT0 , T ) + r1 gd2 (1, −nT0 , T )].

(10.119)

n=1

To evaluate (10.119), the second factorial moment of the single-particle-induced detector count is needed. From (4.37), with a proper change of arguments, it is given as   (1, t, T ) + (t)gd2 (1, t, T ) = gd (1, t, T ) = (−t)gd1

 λ2d λf ν(ν − 1) U (−) (t, T ), t ≤ 0, U (+) (t, T ), t ≥ 0, α3

where U (−) (t, T ) = e αt [1 − 2αTe −αT − e −2αT − (1 − e −αT )2 (e αt − 1)] and U (+) (t, T ) = 1 − 2α(T − t)e −α(T −t) − e −2α(T −t) .

(10.120)

Ch10-I045064.tex

27/6/2007

9: 42

Page 286

286

Imre Pázsit & Lénárd Pál

Using (10.120) in (10.119) and performing the summations leads to the following result for the Feynman Y (T ) function: Y (T ) = =

M2 (T ) − M12 (T ) M1 (T )

 λ2d λf r1 ν(ν − 1) −2αue −αu + 2(1 − e −αT ) [T /T0 ] + 3 M1 (T )α 1 − e −αT0

 e −2αu − e −2αT + (1 − e −αT )2 2αT0 e −αT0 (e −αu − e −αT ) − 1 − e −2αT0 (1 − e −αT0 )2   λ2d (r2 − r12 ) 2(e −αu − e −αT ) e −2αu − e −2αT + (1 − e −αT )2 + [T /T0 ] − + . (10.121) M1 (T )α2 (1 − e −αT0 ) (1 − e −2αT0 ) −

This result can be written in a somewhat simpler form by utilising some identities between the parameters and also re-casting the first moment into M1 (T ) =

r1 λd ∗ M1 (T ) α

with M1∗ (T ) =

1 + [T /T0 ](1 − e −αT0 ) − e −αu . 1 − e −αT0

One then obtains Y (T ) =

λd λf ν(ν − 1) α2 M1∗ (T )



2αT0 e −αT0 (e −αu − e −αT ) −2αue −αu + 2(1 − e −αT ) − + (1 + δ∗ )[T /T0 ] −αT 0 1−e (1 − e −αT0 )2 −αu − e −αT )  e −2αu − e −2αT + (1 − e −αT )2 ∗ 2(e − (1 − δ∗ ) − δ . (10.122) 1 − e −2αT0 (1 − e −αT0 )

Here δ∗ is a ‘source enhancement factor’ δ∗ ≡

r1 (Dq − 1) r1 |ρ| |ρ| = δ − . ν Dν ν Dν

(10.123)

Its form is similar to the factor δ, defined in (10.15) and found in connection with multiple emission stationary and non-stationary sources. The difference is that instead of the factor Dq , the expression Dq − 1 appears in the numerator. The fact that in systems pulsed with instantaneous injection, i.e. with strictly periodic pulses of zero width, Dq − 1 replaces Dq in the second moments expressions was already observed in Chapter 3, in connection with the variance of the particle number. Although the present expressions refer to the variance of the detected neutrons, the difference remains the same. In particular, as it was also discussed in Chapter 3, this difference remains even if the width of the finite pulses is decreased to zero, while the source intensity is increased accordingly. This property of the variance of the neutron number with periodic instantaneous injection is sometimes expressed as if the periodicity of the instantaneous pulsing reduces the variance of the neutron number and the detector counts [60, 65]. An illustration of the realisation is given in Fig. 10.10. This figure was plotted with the same data for T0 and α as those in Fig. 10.5a. A comparison of the figures shows that qualitatively, the shape of the pulse is very similar to the one obtained from using the model of a narrow, but finite source. In particular, the fitting of the theoretical curves to the measured ones has approximately the same performance, irrespective of which model is used. This means that although the correct parameter α can be found also with the model of instantaneous pulses, the shape of the fitted curve shows the same deviations from the measured one as in the case of the model with finite pulse

Ch10-I045064.tex

27/6/2007

9: 42

Page 287

287

Reactivity Measurements in Accelerator Driven Systems

1.2

Y (T )

1 0.8 0.6

T0
0.02, α
266, δ
0.2

0.4 0.2 0

0.01

0.02

0.03

0.04

0.05

Gate width (T )

Figure 10.10 The same Feynman-alpha curve as in Fig. 10.5a, but with periodic instantaneous pulses with T0 = 0.02 s, α = 266 s−1 and δ∗ = 0.2.

width. One would expect that the present model of instantaneous pulses is more flexible, since it contains the extra parameter δ∗ , containing the first two factorial moments of the number of neutrons per injection (i.e. in one pulse of the neutron generator). However, in the fitted curve, the oscillating part still exceeds that of the measurements, indicating that the way the parameter δ∗ appears in the formula does not reduce the oscillating part with the optimum fitting. Presumably the same considerations apply to the difference as those mentioned in Section 10.2.4.

10.4.2 Feynman-alpha with stochastic pulsing Here it is assumed, similarly as in Section 10.2.1, that the arrival of the pulses and the start of the measurement are not synchronised. To extend the description of the previous section to this case, it is assumed that the arrival time of the last pulse that entered the system with t ≤ 0, i.e. outside the measurement period, is equal to t = −ξ, where ξ is a realisation of the random variable ξ which is uniformly distributed in [0, T0 ]. Following the same arguments, the generating function of the distribution of the detector counts is now given as Gd (z, T ) =

1 T0



T0 0

Gd (z, ξ, T )dξ =

1 T0



∞ T0  0

[(T +ξ)/T0 ]

r [gd (z, −nT0 − ξ, T )]

n=0



r[gd (z, k − ξT0 , T )]dξ.

k=1

(10.124) The factorial moments are to be determined by derivations from (10.124). The first moment is simply given as ⎡ ⎤  T0
[(T +ξ)/T0 ] ∞
r 1   ⎣ gd1 (1, −nT0 − ξ, T ) + gd2 (1, kT0 − ξ, T )⎦ dξ. (10.125) M1 (T ) = Gd (1, T ) = T0 0 n=0 k=1

This expression can be evaluated by two different ways, and actually both will be used. For the calculation of the first moment, the simpler way of evaluation goes by recognising that with a suitable change of variables, the sum of the two integrals can be written as two contiguous integrals in the form  0  T r1   g (1, t, T ) + gd2 (1, kt, T ) dt. (10.126) M1 (T ) = T0 −∞ d1 0  and g  in (10.126), the integrals can readily be evaluated and one Using the concrete form (10.115) of the gd1 d2 obtains r1 λd M1 (T ) = T. (10.127) α T0 This expresses the expected result that with random pulsing, the expectation of the detector counts is a smooth linear function of the measurement time, just as in the case of the random pulsing with finite width pulses.

Ch10-I045064.tex

27/6/2007

9: 42

Page 288

288

Imre Pázsit & Lénárd Pál

Calculation of the second factorial moment is more laborious, although still quite straightforward. Apart from the extra variable ξ and the integral, the structure of the expressions will be identical with those in (10.118), including the appearance of the terms n = n . As previously, these can be eliminated by adding and subtracting the terms n = n , leading to r1 M2 (T ) = T0 +

+

+



T0 0

r1 T0 r1 T0 r12 T0



Gd (1, ξ, T )dξ ∞ T0


0





∞ T0
0

 gd1 (1, −nT0 − ξ, T )dξ +

n=0 T0 [(T +ξ)/T
0]

0

⎡



(r2 − r12 ) = T0

T0

⎣

2 gd1 (1, −nT0 − ξ, T )dξ

n=0

(r2 − r12 ) T0



T0 [(T +ξ)/T
0] 0

n=0

 gd2 (1, −nT0 − ξ, T )dξ

n=0 ∞


0

2 gd2 (1, −nT 0 − ξ, T )dξ

 gd2 (1, −nT0 − ξ, T ) +

n=1

[(T +ξ)/T0 ]




⎤2  gd2 (1, kT0 − ξ, T )⎦ dξ.

(10.128)

k=1

In the first four terms on the right-hand side, the sum of the integrals can again be converted into single contiguous integrals, but this is not possible in the last term which contains products of sums. Recognising that the last term in the square brackets is equal to Gd2 (1, ξ, T ), the modified variance can be written as M2 (T ) − M12 (T )

(r2 − r12 ) = T0



0

−∞

r1 + T0



0

−∞

 gd1 (1, t, T )dt

 T r1  gd2 (1, t, T )dt T0 T 0 0 0   T0 2  T0 1 1 Gd2 (1, ξ, T )dξ − Gd (1, ξ, T )dξ . + T0 0 T0 0

+

(r2 − r12 )



2 gd1 (1, t, T )dt

T

2 gd2 (1, t, T )dt +

(10.129)

The very last term of (10.129) is just the square of M1 (T ) which was already calculated and available from (10.127). The reason for writing it in the above form is to show that unlike for the case of stationary sources or for the deterministic pulsing, but in agreement with the case of stochastic finite width pulses, equation (10.53), the last two terms do not cancel each other. This fact leads to some complications, but without principal difficulties. To begin with, the integrals in the first four terms of the right-hand side of (10.129) can be readily performed by using (10.115) and (10.120) for the gd and the gd , respectively, leading to the compact result 

  1 − e −αT r1 λ2d λf ν(ν − 1) (r2 − r12 ) λ2d . T 1 − + T0 α3 αT T 0 α2

(10.130)

Calculation of the second last term, the integral of Gd2 (1, ξ, T ), poses some more difficulties. By recalling (10.124) and (10.125), it can be written as

1 T0



T0 0

r1 Gd2 (1, ξ, T )dξ = T0



T0 0

⎡ ⎤2 [(T +ξ)/T0 ] ∞

  ⎣ gd1 (1, −nT0 − ξ, T ) + gd2 (1, kT0 − ξ, T )⎦ dξ. n=0

k=1

(10.131)

Ch10-I045064.tex

27/6/2007

9: 42

Page 289

289

Reactivity Measurements in Accelerator Driven Systems

This expression can be evaluated by first calculating the sums and then squaring and integrating. Performing the sums leads to Gd (1, ξ, T ) =

∞
n=0

 gd1 (1, −nT0 − ξ, T ) +

[(T +ξ)/T0 ]




 gd2 (1, kT0 − ξ, T )

k=1



r1 λd = [(T + ξ)/T0 ] + α

 − e −α (T ,ξ) , 1 − e −αT0

e −αξ

(10.132)

where (T , ξ) = T + ξ − [(T + ξ)/T0 ]T0 .

(10.133)

The integral of the square of (10.132) can be performed by noticing that, by writing again T = kT0 + u, 0 ≤ u ≤ T0 , or in other words, u = T − [T /T0 ]T0 , one has [(T + ξ)/T0 ] = k + (ξ + u − T0 ) and in a similar way

 (T , ξ) =

u + ξ, u + ξ − T0 ,

ξ ≤ T0 − u, ξ > T0 − u.

(10.134)

(10.135)

The integration of the six terms of the square of (10.132) is rather tedious and lengthy, but the reward is a stunningly simple expression. Summing up all terms, subtracting M12 (T ) and dividing with M1 (T ) yields the final result   % λd λf ν(ν − 1) $ 1 − e −αT ∗ 1 + δ 1 − Y (T ) = α2 αT   r1 λd u(T0 − u) e α(u−T0 ) + e −αu − e −αT0 − 1  + . (10.136) + α T0 T αT 1 − e −αT0 Similarly to the results of the stochastic pulsing with finite pulses, this expression consists of two terms; a smooth term corresponding to the traditional Feynman-alpha solution with stationary sources, although with an amplitude enhanced with the factor δ∗ , defined in (10.123), and a non-negative oscillating term with an oscillation period equal to T0 . In both cases, the smooth, traditional term constitutes a lower envelop of the oscillatory part. This expression is in a close analogy with the one obtained for the case of stochastic finite width pulses with compound Poisson statistics within the pulse, expression (10.106), which also contains a traditional term enhanced with the source multiplicity, and an oscillating part. Another, although more remote, similarity is that in the stochastic pulsed case with finite width pulses, the oscillating term contains the source intensity, which is not present in the case of stationary sources or that of deterministic pulsing. For instantaneous pulses no source intensity can be defined, but in (10.136) the expression r1 /T0 occurs in the first term of the oscillating component, which plays a similar role. The difference between the two formulae is that the source enhancement factors are different. While physically this is a rather significant difference, concerning the use of the formulae for the determination of the reactivity, the two solutions seem to be equally suitable. Each of them contains a separate parameter for the smooth and the oscillating parts, that can be determined from the fitting procedure. The formula for the instantaneous injection has the clear advantage that its form is much simpler than that of the finite pulse case, (10.106), the latter containing the searched parameter α not only in the exponent of the smooth part, but also in each coefficient an and bn of the oscillating part. Because of its much simpler form, the instantaneous pulsing formula seems to be superior in applications to reactivity determination in pulsed subcritical systems. One illustration of the stochastic instantaneous pulsed Feynman-alpha curve is shown in Fig. 10.11. Again, a very good similarity is seen with the corresponding curve in Fig. 10.6d calculated with the same T0 and α parameters.

Ch10-I045064.tex

27/6/2007

9: 42

Page 290

290

Imre Pázsit & Lénárd Pál

1.2

Y (T )

1 0.8 0.6 0.4

α
598, T0
0.02

0.2 0 0

0.01

0.02 0.03 Gate width (T )

0.04

0.05

Figure 10.11 Feynman-Y curve for stochastic instantaneous pulsing, for the same case as in Fig. 10.6d. T0 = 0.02 s, α = 598 s−1 . The broken line shows the traditional Feynman-Y curve, corresponding to the first term of (10.136).

10.4.3 Rossi-alpha with stochastic pulsing Here the concepts applied remain the same, but the formalism will change somewhat. It is assumed that the measurement is triggered by a detection at time t = 0, and is completed by another detector count at t = τ. The detections take place in infinitesimal intervals dt and dτ around t = 0 and t = τ, respectively, hence the detection intensity of the joint detection is λ2d times the joint expected value of the neutron number at times t = 0 and t = τ. As in the previous case, at time t = 0 the system is already in a periodic stationary state, since the pulsing started at t = −∞. For clarity of the derivation, the intensity of the joint detections at t = 0 and t = τ will be denoted as P2 (0, τ) = λ2d N(0)N(τ), such that P2 (0, τ)dt dτ = λ2d N(0)N(τ)dt dτ is the probability, in first order of dt dτ, of having one detection within (0, dt) and another one in (τ, dτ). Again it is practical to consider the events before and after t = 0 separately. It will be assumed that the last pulse for t ≤ 0 arrived at t = −ξ, where ξ is a random variable distributed uniformly in [0, T0 ]. With these preliminaries, and accounting for the additive character of the expected values, P2 (0, τ) can be written as follows:  ∞ λ2d T0
nQ(n, ξ) X (n, τ, ξ)dξ, (10.137) P2 (0, τ) = T0 0 n=0 where X (n, τ, ξ) =

∞


m[P(m, τ|n − 1, 0) + R(m, τ|0, 0, ξ)].

m=0

Here Q(n, ξ) is the probability that at time t = 0 there will be n neutrons in the system on the condition that the last pulse for t ≤ 0 arrived at t = −ξ; P(m, τ|n, 0) is the probability that in a source-free system there will be m neutrons at time t = τ given that at time t = 0 there were n neutrons present (the argument n − 1 in (10.137) accounts for the detection of one neutron at t = 0); and finally R(m, τ|0, 0, ξ) is the probability that in a pulsed system, there will be m neutrons at time t = τ given that the last pulse for t ≤ 0 arrived at t = −ξ, and that there were no neutrons (n = 0) at time t = 0 in the system. This equation can be converted into one for the corresponding generating functions on the right-hand side, by introducing
FQ (z, ξ) = Q(n, ξ)zn (10.138) n

and FR (z, τ, ξ) =


m

R(m, τ|0, 0, ξ)zm ,

(10.139)

Ch10-I045064.tex

27/6/2007

9: 42

Page 291

291

Reactivity Measurements in Accelerator Driven Systems

whereas the generating function of P(m, τ|n − 1, 0) is equal to g(z, τ)n−1 where g(z, τ) is the basic generating function of the probability of finding n particles in the system at time τ given that there was one particle in the system at t = 0. This generating function was defined in Chapter 1 and its first and second moments are given in (1.57) and (1.59) as g  (1, τ) = e −ατ

(10.140)

and λf ν(ν − 1) −ατ (10.141) e (1 − e −ατ ). α Here account was taken of the fact that the α in the above is equal to −α of Part I, and that Qq2 = λf ν(ν − 1). With the above definitions, (10.137) can be re-written as g  (1, τ) =

P2 (0, τ) =

λ2d T0



T0 0

  {FQ (1, ξ)g  (1, τ) + FQ (1, ξ)FR (1, τ, ξ)}dξ.

(10.142)

Equation (10.142) can be evaluated by noticing that the generating functions FQ and FR can be constructed, in a similar way as in the previous section for the Feynman-alpha method, from g(z, τ) as follows: FQ (z, ξ) =

∞ 

r[g(z, nT0 + ξ)]

(10.143)

n=0

and 

[(τ+ξ)/T0 ]

FR (z, τ, ξ) =

r[g(z, τ + ξ − kT0 )].

(10.144)

k=1

The calculation goes exactly along the same lines as in the previous section. Applying the same tricks as before, i.e. adding and subtracting the terms n = n and converting the single sums of integrals from 0 to T0 into single integrals from 0 to infinity, one arrives at ⎧  ⎨ λ2 (r − r 2 )  ∞ λ2d r1 ∞  1 2 d 2 g (1, t)dt + g (1, t)dt P2 (0, τ) = ⎩ T0 T0 0 0 2 ⎫  ∞ ⎬ λ2d r12 T0
 + g (1, nT 0 + ξ) dξ e −ατ ⎭ T0 0 λ2 r 2 + d1 T0

⎡



T0 0

⎣

n=0

∞
n=0

g  (1, nT 0 + ξ)dξ




[(τ+ξ)/T0 ]

⎤ g  (1, τ + ξ − kT0 )⎦ dξ.

(10.145)

k=1

The integrals in the first two terms and the summations in the third and the fourth are carried out easily. One finds that the last term, containing the product of two sums, can be written as the difference of two terms, out of which one cancels out with the third term of (10.145). The remaining term contains the construction e −αξ e −α(τ,ξ) , (1 − e −αT0 ) where the function (τ, ξ) is the same as in (10.133) and (10.135), where now the running variable 0 ≤ u ≤ T0 is defined by u = τ − [τ/T0 ]T0 . The integral of this term with respect to ξ can be evaluated with the technique

Ch10-I045064.tex

27/6/2007

9: 42

Page 292

292

Imre Pázsit & Lénárd Pál

1.2 1

R (τ)

0.8 α
266, T0
0.02

0.6 0.4 0.2 0

0.02

0.04

0.06

0.08

Gate delay (τ)

Figure 10.12 Rossi-alpha curve for stochastic instantaneous pulsing. T0 = 0.02 s, α = 266 s−1 . The broken line shows the traditional Rossi-alpha curve, corresponding to the first two terms of (10.148).

described in the previous section. Summing up all terms leads to the final result λ2d P2 (0, τ) = 2αT0



 $ −αu % r1 λf ν(ν − 1) r12 2 −ατ −αT0 αu . e + +e e + r2 − r1 e α 1 − e −αT0

(10.146)

This expression contains the expected decaying exponential, plus an oscillating part with non-zero mean, which shows periodical stationarity. The expression can be converted into the more familiar form of the Rossi-alpha formula, expressed in the form R(τ) of (9.163), by taking into account that the detection rate Z reads in this case, from (10.127), as r1 λd . (10.147) Z= αT0 Moreover, it gives insight to split up the last, oscillating term into a constant value expressing the time average over the oscillation period, plus a term oscillating with a zero temporal mean. This leads to R(τ) =

λd λf ν(ν − 1) r 1 λd P2 (0, τ) = [1 + δ∗ ]e −ατ + Z 2α αT0  r1 λd 2(1 − e −αT0 ) −αu −αT0 αu + e +e e − , 2(1 − e −αT0 ) αT0

(10.148)

where the enhancement factor δ∗ was defined in (10.123). The result (10.148) contains a traditional smooth Rossi-alpha formula, i.e. a decaying exponential and a constant part corresponding to the so-called correlated and uncorrelated counts, plus an oscillating part. Comparison with (10.147) shows that the constant part is equal to the detection intensity in one detector, i.e. it agrees with the corresponding term in the traditional Rossi-alpha formula. The oscillating part has a discontinuous derivative, unlike in the case of the stochastically pulsed Feynman-alpha curve. Qualitatively the solution is similar to that obtained by the method of finite but narrow pulses, even if some of the factors appearing have different values. One illustration of (10.148) is given in Fig. 10.12. As both the formula and the figure show, in contrast to the pulsed Feynman-alpha method, but in agreement with the findings of the treatment with finite width pulses, the oscillations do not decay with time, rather they are periodic. Qualitatively a very good agreement can be seen in a comparison with the results shown in Fig. 10.7, referring to the case of finite but narrow pulses. This lends the possibility of eliminating the oscillating part from an experiment, as was demonstrated earlier, from the tail of the curve at large τ values. For an analysis of the difference between the results of finite width pulses and instantaneous injection, one has to compare (10.148) with the case of a pulsed compound Poisson source, (10.107). The comparison shows that the basic difference between the inhomogeneous compound Poisson source and the instantaneous injection is that in the former the source enhancement factor δ appears, whereas in the latter the factor δ∗ .

Ch10-I045064.tex

27/6/2007

9: 42

Page 293

Reactivity Measurements in Accelerator Driven Systems

293

This is again the same difference as the one observed earlier in the comparisons of finite narrow pulses and instantaneous pulses. The conclusions regarding the applicability of the two formulae are also the same as before. Whether the finite width pulse or the instantaneous injection description is a better model of the actual physical situation cannot be decided from purely mathematical principles. Regarding the applications of the results for the evaluation of measurements, the difference is insignificant. In that process, both δ and δ∗ are fitting parameters, hence there is no difference which formula is being used. From the practical point of view, since (10.148) contains a closed form expression, it might be more advantageous to use than (10.107) when evaluating measurements with narrow pulses.