SPATIALLY DEPENDENT INTEGRAL
NEUTRON TRANSPORT THEORY FOR
HETEROGENEOUS MEDIA
USING HOMOGENEOUS GREEN'S FUNCTIONS
By
JOHN PHILLIPS CHURCH
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE
DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
April, 1963
To my former teacher
Dr. Dudley E. South
ACKNOWLEDGMENTS
The author wishes to express his appreciation
to his chairman, Dr. G. R. Dalton, for his guidance and
assistance throughout the course of this investigation.
In addition, the author is indebted to the
staff of the Computing Center of the University of Florida
for their assistance, and to the Computing Center and the
Department of Nuclear Engineering for their combined
financial support of the computational work for this
dissertation.
The author has been ably assisted in the prepa
ration of this dissertation by his wife, Nancy, who typed
the numerous drafts and the final manuscript; its final
form owes much to her concern. It is difficult to express
in words the author's appreciation of her faith in him.
She has made the path to his goal less difficult.
TABLE OF CONTENTS
ACKNOWLEDGMENTS . . . . . . .
LIST OF ILLUSTRATIONS . . . . . .
LIST OF TABLES . . . . . . . .
Y TI D~ T I
.N1 TR .W I A
CHAPTER I
CHAPTER II
I1.1
11.2
CHAPTER III
II1.1
111.2
111.3
CHAPTER IV
IV.1
IV.2
CHAPTER V
V.1
V.2
Page
S. . iii
* . vi
S. . viii
1
. . . . . . . # .
DEVELOPMENT OF THE HOMOGENEOUS
MEDIUM GREEN'S FUNCTION METHOD . .
ANGULAR INTEGRATION . . . .
Expansion of Angular Dependence in
Spherical Harmonics and Zeroth
Harmonic Approximation ..
Advantages and Disadvantages of
Homogeneous Green's Function Method
SPEED DEPENDENCE . . . . .
General Speed Dependent Equation .
The Source Term . .. . .
OneSpeed Equation . . .
DIFFERENCE DENSITY METHOD . .
Development of Difference Density
Equations . . . . . . .
Equivalence of Equations for N and D
PLANE SLAB GEOMETRY . . . .
OneDimensional Problem . . .
Symmetry Considerations . . .
* J.
S 4
S 15
S 15
S 26
28
S 28
S 30
S 34
S 36
*
.
.
CHAPTER VI
CHAPTER VII
VII. 1
VII.2
VII.3
CHAPTER VIII
VIII.1
VII.2
CHAPTER IX
IX.1
IX. 2
IX.3
IX.4
NOMENCLATURE
APPENDIX A
APPENDIX B
REFERENCES
COMPARISON OF N AND D METHODS
FOR SLAB GEOMETRY . . . . .
THE FIRSTFLIGHT GREEN'S FUNCTION
FOR A HOMOGENEOUS MEDIUM . . .
Differential Equation and Boundary
Conditions . . . . .
Analytic Form of G(v;x x') . ...
Monte Carlo Generation of G(v;xlx')
SOLUTION OF THE TRANSPORT EQUATION .
Spatial and Speed Dependence . .
Iteration Techniques . .. . .
RESULTS AND CONCLUSIONS . . .
Comparison of the Homogeneous Green's
Function Method with High Order P
and Sn Methods . . . . .
Monte Carlo Generation of the
Green's Function . . . . .
Iteration Convergence . . .
Generation of the Green's Function
for the Multispeed Problem . . .
. . . . . . . . *
DERIVATION OF Gf(v;xIx') . .
SAMPLE PROBLEMS . . . .
4 0 * 0 0 0 . . . .
S 69
* 77
S 90
. 102
. 102
. 115
. 119
119
125
132
135
137
142
145
169
* 0
* 0
LIST OF ILLUSTRATIONS
Figure Title Page
1 The Spatial Dependence of the h Functions
of Equations (2), (7), and (8) . . . 13
2 Plane Geometry Transformation . . . 50
3 Symmetry of Neutron Flow'at the Boundaries
of Unit Cells in An Infinitely Repeating
Lattice . . . . . . .. . . 75
4 Slab Unit Cell with Mirror
Boundaries at x+b . . . . . .. 85
5 Mirror Symmetry of FirstFlight Green's
Function About Plane at x0 . . . . 97
6 Average Values of the WholeCell and
HalfCell Green's Functions for
Problem 2, Appendix B . . . . . 128
7 Average Value of the WholeCell and Half
Cell Green's Functions for Problem 4,
Appendix B . . . . . . . . 129
8 Effect of Convergence Acceleration
Techniques Applied to Problem 2 of
Appendix B. Solution Iterated Until
Residuals Were Equal to 105 of the
Neutron Density at Each Spatial Point . 133
9 Spatial Dependence of Scalar Neutron
Density for GoldGraphite Cell of Problem 1. 148
10 Spatial Dependence of Scalar Neutron
Density for FuelWater Cell of Problem 2 .155
11 Advantage Factor for FuelWater Cell of
Problem 2 by Various Approximation Methods 156
12 Advantage Factor Calculated by HGI
Method vs. Number of Iterations for
FuelWater Cell of Problem 2 . . . . 158
13 Advantage Factors for Problem 3 by
Various Calculations . . . . . . 161
14 Advantage Factor for Problem 4 by
Various Calculations . . . . . .. 165
15 Spatial Dependence of Scalar Neutron
Density for Cases 1 and 2 of Problem 4 . 167
16 Spatial Dependence of Scalar Neutron
Density for Cases 3 and 4 of Problem 4 . 168
vii
Title
Page
Figure
LIST OF TABLES
Table
Number Title Page
1 Relative Effort for Methods of Solution
of EndPoint and Midpoint Equations for
a Total of M Spatial Intervals . . . 109
2 Breakdown of Computation Time for
Problem 2, Appendix B . .. . .. 134
3 Parameters for Problem 1 . . . . 147
4 Parameters for Problem 2 . . . . 153
5 Parameters for Problem 3 .. . . 159
6 Parameters for Problem 4 . . . . 163
viii
INTRODUCTION
Probably the most fundamental information needed
to describe static behavior of a nuclear reactor is the
speed and spatially dependent neutron collision rate,
that is, Z(r,v)N(r,v)v. Such information is needed in
order to predict longterm kinetic behavior of the reactor,
burnup, power level, etc. Thus.one needs to know the
neutron density, N(r,v), over the entire speed spectrum
and space of the reactor.
Because of chemical binding and crystalline
effects, the scattering kernels are considerably different
for neutron speeds near thermal than for the higher speeds
implied when one speaks of epithermal and fast speed
ranges. In general, the absorption crosssections of
reactor fuels also have a different functional dependence
on neutron speed at thermal speeds than at higher speeds.
Finally, because of the tremendous range of speeds in
reactor spectra, and because a critical reactor with no
independent sources is an eigenvalue problem if the entire
speed spectrum is considered, the problem of predicting
the neutron density is usually divided into several speed
ranges.
The present paper deals exclusively with the
prediction of the spatial and speed dependence of the
steadystate thermal neutron density in heterogeneous
media, although the basic method to be presented is
equally valid for any speed range.
There have, of course, been numerous approaches
to various aspects of the problem which can be very
loosely divided into two classes according to whether
the integrodifferential or the integral form of the
transport equation was solved. The methods of solving
the integrodifferential equation differ to a large extent
only in the choice of functions used to expand the angular
dependence of the terms in the equation. The expansion
functions include spherical harmonics (13), Tchebycheff
polynomials (47), and first order polynomials, that is,
trapezoids (8, 9).
The methods of solving the integral equation
have been a bit more varied, due perhaps to a more
intuitive derivation of the equation to be solved. As
Osborn (20) has pointed out, some approaches, particularly
as applied to thin absorbers (10, 11), have separated the
problem into two partsfirst, a calculation of the
absorption rate in the absorbing media in terms of the
density at the surface of the absorber and second, a
density depression problem of evaluating the density
at the surface of the absorbing media in terms of the
density in an allmoderating media. Other integral
approaches (1214), based on transport theory, depend
on the fact that the absorption rate depends only upon
the zeroth angular moment of the neutron density, the
solution for which is independent of higher angular
moments if scattering is isotropic.
The present paper describes the latter type
of approach. In brief, the method to be presented uses
a kernel for a homogeneous finite medium to formulate an
integral transport equation which predicts the speed and
spatially dependent neutron density in a heterogeneous,
finite unit cell.
CHAPTER I
DEVELOPMENT OF THE HOMOGENEOUS
MEDIUM GREEN'S FUNCTION METHOD
Assume that the neutron density in a general
threedimensional medium satisfies the time independent
linear Boltzmann equation, and let N(r,vn)drdvd. denote
the steadystate angular neutron distribution, that is,
the number of neutrons in the volume element dr about r
having speeds in dv about v, and whose directions of
motion lie within the solid angle dn about n; N(r,vn)
will subsequently be called the angular density.
Define V(r,v) to be the probability per unit
time that a neutron of speed v will suffer a scattering
collision in dr about r with a nucleus which is in
motion with respect to the laboratory system. Further,
define y(r,v) to be the probability per unit time that
a neutron of speed v will be absorbed in dr about r by a
nucleus which is in motion with respect to the laboratory
system. Finally, define P(r,v,v4*nn) dv dgto be the
probability per unit time that a neutron of velocity v'r.'
will suffer a scattering collision in 4r about r with a
nucleus which is in motion with respect to the laboratory
system such that the neutron emerges from the collision
with speed in dv about v and direction of motion in dn
about n. Then let S(r,vr)drdvdn denote the source term,
that is, the number of neutrons per unit time being added
to the volume element dr about r, having speeds in dv
about v, and whose directions of motion lie within the
solid angle dn about n.
The transport equation which takes into
account the relative velocity between the neutrons
and collision nuclei which are in thermal motion can
then be written as
vn 7N(r,vn + [V(r,v)+ A(r,v)] N(r,vr
S(r,vrn) + fdv*d' P(r,vN(v,. N(r,v'rD. (1)
Equation (1) is an integrodifferential equation; to
make the problem of its solution determinate it is
necessary to specify boundary conditions. The boundary
conditions follow from the physical interpretation of
N(r,vn). For example, at an interface between two media,
any packet of neutrons characterized by the vectors
r + Rn and vn will contain exactly the same number of
neutrons when it enters one medium as when it left the
other; in other words, N(r + Rn,vn) is a continuous
function of R for r + Rn at the interface. For a
complete discussion of this and other boundary conditions
applied in neutron transport theory, the reader should
consult Davison (2).
Consider the space of interest to be divided
into R regions, the medium of each region being assumed
homogeneous (the derivation of Eq. (1) assumes that the
medium is isotropic) so that the spatially dependent
probabilities V, ), and P in Eq. (1) are constant with
respect to r throughout each individual region.
Define a unit function, h, such that
h(ps) 1 if p > s
= 0 ,if p s. (2)
For each interface between regions define an index k and
a position vector T such that 1k_ e JA. Then
h(rkr) h(Eklr) = 1, if Ekr1 r Ek
(3)
= 0, if r e k, or if r = rk
where r is on the line between rk and :kl. One can
then use Eq. (3) to separate the r dependence of the
probabilities V, Y, and P in Eq. (1) in a special way.
In general, for any operator 0, one can write that
R
S(r,) kl [h(Ekr)h(klr)] (k,) (4)
Use Eq. (4) to rewrite Eq. (1) in the form
R
yvt* N(r,vn) + [ [h(rkEr)h(rk_,r)]
S[v(k,v)+ (k,v)] N(r,vn)
R
8(r,vr) + 1 [h(rkr)h(!k1)]
dv'/ df' P(k,vIv,n',O)N(r,v'n. ') (5)
Choose a particular region, say k', and remove
the term for the k' region from the summation on the LHS
of Eq. (5), meanwhile placing all other terms of that
summation on the RHS. Thus Eq. (5) becomes
vrL N(r,v*) + [h(Ek,r)h(rk'_r1)
[V(k',v)+ /(k,v)] N(r,vn)
R
S(r,v) [(r)h( r]
k+k'
[V(k,v)+ y(k,v)] N(r,vn)
R
+ 1 [h(Ekr)h(k_1r)]
ki
dv' dn P(k,viv,n'~N(r,v'n') (6)
Now using the identity
h(ps) I h(sp) (7)
one can write
h(!k,r) h(Ek, r) [1 h(rk,)] h(k,_lr)
1 [h(rr,.)+h(r,_r)] (8)
Substituting Eq. (8) into Eq. (6) and putting the
[h(rr,) + h(Ek'~,i)] term on the RHS gives
va*VN(r,vm) + [V(k',v)+'(k',v)] N(r,vn)
SS(rvn) + [h(r,)+h(r .,r)
[V(k',v)+ '(k',v)] N(r,vn)
[ [h(rkr)h( ~."r) [V(k,v)+ '(k,v)] N(r,vn)
k+k'
R
+ E [h(rkir)_h(!k1L]
(9)
dv' dn' P(kv',v,nt.) N(r,v'n')
Define a function gk' (r,vLr' ,vtn.') by requir
ing that it satisfy the equation obtained by replacing
the RHS of Eq. (9) by the Dirac delta function product
6 (rr') d (vv') & (nn'), subject to the same boundary
conditions as Eq. (9). Thus
TVgk T (rn r'Iv' 1)
+ [v(k',v)+ W/(k',v)] g' (r,vnlr' ,v'i')
g(rr')&(vv')s(nn') (10)
The function gk (r,vnCr',v'.n') has a physical
significance; gk (r,vnr r',v'n') drdvdn is the number of
neutrons in the volume element dr about r with speeds
in dv about v and with directions of motion lying within
the solid angle dn about n due to a unit point source
at r' emitting one neutron per unit time in direction
n' with speed v' in a firstflight medium which has
the same properties as the medium in region k'.
It is important to understand fully the term
"firstflight medium" as used here. To be more explicit,
a firstflight medium is a medium in which all collisions
with nuclei result in removal of the neutron from the
population. Thus the removal probability for the Green's
function problem is equal to the total collision proba
bility in medium k' of the real problem, and the neutrons
denoted by gkt(rvn l r*,v') are firstflight neutrons
in the sense that they have suffered no collisions
enroute from the source point to the field point. The
term firstflight neutrons will be used throughout the
following work.
The subscript k' will be dropped from
gk (r,vr'iE*,v't) in the following work, it being
understood that the function depends on the properties
of the particular region, k', chosen.
The function g(rvnjr',v'n') can be used as
a Green's function to convert Eq. (9) to the following
integral equation (see Chapter 7 of reference (16)).
N(r,vn) dr' dv' dn' g(r,vn r v'n) ,v'f)
rt v' n'
+ [h(r'Ek)+hc(Ek.,l') [V(k',v')+ 2'(k',v')]N(r',v'n')
[h(!kr' )h(rklr')] [V(k' vt)+ Y'(k',v')] N(r' ,v')
k+k'
R
+ [h(Ekr)h(kr1]
dv"d P(kvyrv',n")N(r'vlg) .
Since Eq. (10) describes the neutron density
in a medium in which all collisions result in removal
of the neutron from the population, the only way a neutron
which is emitted at r' may arrive at r is by direct
flights, that is, with no collisions enroute. It will
be assumed and later justified that events at the bounda
ries do not change the speed of the neutron. Hence if a
neutron of speed v' is emitted at r' and arrives at r,
its speed must still be v'.
From this discussion it can be concluded that
the speed dependence of g(r,vfir',v'n') includes a delta
function in speed so that
g(r,vn r' ,v'') g(v;r,~ r' ,n.' )(vv'). (12)
The product solution, Eq. (12), is a result of the fact
that since there can be no scattering in the medium,
then there is no interdependence of nr*n' and vv'. Note
that g(v;r, lr',') ) still has an implicit dependence on
v due to the v dependence of V(k',v) and ((k',v) in
Eq. (10).
Substituting Eq. (12) into Eq. (11) and carry
ing out the integration over v' yields
N(rvn) dr d g(v;r,nl r,f' )S(r,vn')
+[h(r'rk )+h(r, lrt)] [V(k' v)+ Z(k',v) N(r',vn)
Sh(Ekr)h(r'i )][V(k,v)+ (k,v)] N(r',vn')
k+k'
R
+ [h(rkr')h(k1r')] (13)
k1l k  
dv" dn" P(k,v"V ,n"nr )N(r',v"n ) .
The integration over r' may now be considered.
The meaning of the functions [h(r'rk,)+h(Ek1r' ) and
[h(rk,r')h(rk_'r') may be easily deduced from Eqs.
(2), (7) and (8) and is illustrated in Figure 1. The
function [h(r'!k,)+h( k,lr')] has unit value every
where except in region k' where it has zero value, while
[h(kr')h(Eklr')] has zero value everywhere except
in region k where it has unit value. It can easily be
seen that the r' integration of the second term on the
RHS of Eq. (13) is an integration over all regions except
1
0 
1
h(r'k +r' )
0 
1
h(r'r) + h(rkr')
0 
1
h(r~ktre )h(rk' _1r')
0 
Figure 1 The Spatial Dependence of the h Functions
of Equations (2), (7), and (8)
the k'th one while the r' integration of the remaining
stepfunction terms is a sum over k of r' integration
over each individual region k. Thus one can write
Eq. (13) as
N(r,,vn) dr' dn' g(v;r,, r',n')S(r',vn')
rt n
R
+ ki [v(k',v)+ (k',v)] V(k,v)+ k,v)]
dr' dn' g(v;r,nrr',nt)N(r',vn')
r'in k r'
R
+ Z dr' dnY g(v;r,n jr,nt)
k1 ~ 
r'in k n' (14)
dv" d"P(k,v" v,n'L.' )N(r',v"r") .
Note that the kk' term in the first summation
is automatically zero, so that no special notation is
needed on the summation symbol.
CHAPTER II
ANGULAR INTEGRATION
11.1 Expansion of Angular Dependence in Spherical
Harmonics and Zeroth Harmonic Approximation
The angular integration of Eq. (14) may be
performed as follows. Expand all the angular dependent
terms of Eq. (14) in the complete orthonormal set of
spherical harmonics as defined in Chapter 7 of reference
(15). Thus
oo n
N(r,vn) Nm(r v)y(n) (15)
n0O mn
co n
S(r,vn) E Sm(r,v)Ym(n) (16)
nO mn n n
oo a
P(k,v'"Uv,n'!n) Pa(k,vv)yb()b*(). (17)
a0 ba a a
The expansion for the Green's function is more
complicated. Because of possible events at the bounda
ries enclosing the space of interest, g(v;r,nI r' ,')
depends explicitly on both n and nr' and not just on r*n'
as P(k,vn.v,n'n ) does. For the moment only the n.' de
pendence of g(v;r,n r',nh) is represented by an expansion
in the complete orthonormal set of complex conjugate
spherical harmonics. Thus
g(v;r,n r'l,n') gs;rn r')Y (). (18)
s0 ts
Substitute Eqs. (15) through (18) into Eq. (14) to obtain
N(r,v ndr' dn' gt (v;rnJf)Yt* (A) : m(r',v)Ym(n!)
Ss,t nm
r' n'
+ R v(k',v)+ (kt,v)0 [V(kv)+ Y(k,v)
Sfdr' dn E gg(v;r, lr')Y t*(n) N(r (n'
S s,t s n,
r'in k n'
+ rl drt dn' g (v;r,nr')Y (')
r'in k n1'
(19)
ab a a n 
Equation (19) may be simplified considerably
by interchanging the order of integration and summations
and by using the orthogonality relation of the spherical
harmonics, namely,
n)Yb n, a Smb (20)
where Sn,a is the Kronecker delta function.
After carrying out the angular integration,
Eq. (19) becomes
N(r,vn) dr' E ge(v;r,nlr')Sm(r',v)
n,m
rf
+ [V(k ,v)+ ](k ,v) [V(k,v)+ ?'(k,v)1
Sdr g gm(v;r,nr')Nm(r',v)
n,m
r'in k
(21)
+ R dr' dv" gn(v;rn l')Pn(k,vZv)o(r',v")
kl 1 f n,m
r'in k
Equation (21) is an expression for the total angular
neutron density and is exact; no approximation has been
made so far.
Now, the relative isotropy or anisotropy of
the n' dependence of g(v;r r',n') depends on the
boundary conditions applied to Eq. (10) and also on the
properties of region k'. Recalling the physical meaning
of g(v;r,nlr ,n.'), one observes' that n' denotes the
angular dependence of the unit point source in a medium
in which all collisions result in removal of the neutron
from the population. If the space of interest is of
infinite extent, then g(v;r,n r',n') has a double delta
function behavior in direction, namely S(n*') (r')] ,
which merely describes the fact that since scatters
are not considered in the calculation of the Green's
function the direction of motion at r not only must be
the same as the direction of motion when emitted at r',
but also must be pointed along the vector connecting the
field and source points. In this case, the complete set
of coefficients g (v;r,,nir) would be required in Eq. (18)
to satisfactorily describe g(v;r,n r',n').
The problem is entirely different, however, if
one imposes perfectly reflecting boundaries around the
space of interest including the source as would be done
if one were calculating a unit cell in a repeating array.
Collisions or events at the boundaries will reflect
neutrons back into the space of interest with some
direction nL. If the neutron then arrives at r it will
have direction n, such that n + n' and n (rr'),
although there llbe a definite relation between
although there will be a definite relation between
, n'_, and (rr*). Further, if the medium characterized
by [V(k',v)+ ,(k',v)] has a mean free path that is long
relative to the distance between boundaries, then the
neutron may be reflected at more than one boundary before
arriving in the element of volume dr about r.
Now the integral of the angular density due to
a unit point source per solid angle, integrated over all
directions of emission n', is simply the angular density
due to an isotropic point source of one neutron per unit
solid angle or a total of 41T neutrons and is denoted by
g(v;r,nar'). Thus
g(v;r,nlr') g(v;r,nlr',n')dn'. (22)
From the discussion above regarding reflections at the
boundaries, it is obvious that in a finite space enclosed
by reflecting boundaries an isotropic point source will
produce a very nearly isotropic angular density
g(v;r, i r').
In order to take advantage of the relative
isotropy of g(v;r,n r'), one proceeds as follows. If
one asks only for information about the number density
N(r,v), the dependence of Eq. (21) may be integrated
out to obtain
N(r,v) dr' Z g(v;rlr )Sm rv)
nm
r'
k V (k ,v)+ k',v) [V(k,v)+ /(k,v)
dr' g(v;r r)lE (r',v)
n,m
r'in k
(23)
R m
+ dr' dv" n g(Vr'r)P (kv"v)Nm(r',v")
k1 n,mn n 
r'in k
where
N(r,v) /N(r,v)dr (24)
n
Expanding the angular dependence of
g(v;rl r') in spherical harmonics gives
g(v;rn ') g(v;r r')Ym(n) (25)
n,m
The coefficients gn(v;rlr') in Eq. (23) are recognized
as being equal to the coefficients in the angular
expansion of g(v;r, nr'), a function which has been
discussed above and which has been concluded to be
nearly isotropic for a finite space enclosed by perfectly
reflecting boundaries and containing a medium having a
mean free path long relative to the distance between
boundaries.
Thus the coefficients gn(v;r r') form a rapidly
decreasing sequence for increasing n. Further, since
Eq. (1) includes a scattering term, the angular density
N(r,vn) should be rather isotropic, even for rather
strong absorbers. A unique advantage of the integral
equation method over the integrodifferential equation
method thus becomes apparent; namely, the zeroth harmonic
approximation disregards only products of small terms,
that is, one assumes
1 1
1 1
0 0 0
g0 N0 g1 N *
g1 N1
1 1
Due to the orthogonality of the spherical harmonics,
there are no cross products to consider such as g0 No
or go N1 which might be of the same order of magnitude
0 0
as g0 NO'
The situation with regard to the scattering
term in Eq. (23) is even more favorable. One need not
assume completely isotropic scattering, but merely make
the less stringent assumption that the product of the
zero order coefficients of the Green's function, the
angular density, and the scattering kernel, is much
greater than the product of higher order coefficients,
that is,
g1 PI Ni1
g0 P NO 9 10 P NO
gl p N1
1 11
On the basis of the above discussion the
zeroth harmonic approximation is made in Eq. (23) and
only the n0 terms are kept. Thus Eq. (23) becomes
N(r,v) dr' g0(v;r r')So(r',v)
r0
+ [v,(k.v)+/(kv) [V~f,v)+(kv
dr' g00V;r Ir)N (r',v)
r'in k
R O l o t
+ Z dr' dv" go(v;r r')P.(k,v.v)NO(r,v"). (26)
+1 kl d 
r'in k
The zero order coefficients may be put in a
more recognizable form as follows. Substitute Eq. (15)
into Eq. (24) to obtain
N(r,v) Na(r,v) Ya(n)dn
n,m
Nm(r,v) Ym(c) /f YO*()
n,m
/41 N(r,v)
00*
where use has been made of the fact that Y* (n)
 YO n) 1/1 Thus
N0(rv) N(r,v) (27)
Similarly, defining
S(r,v) s (r,vn)dc (28)
rL
and using Eq. (16), one obtains
o 1
S(r,v) S(r,v) (29)
0 
Define
g(v;rlr) g(v;r,lr')dn (30)
to be the number density of firstflight neutrons at r
having speed v due to a point source at r' emitting 4TT
neutrons isotropically with speed v in a homogeneous
finite medium. Substitution of the angular expansion
Eq. (25) into Eq. (30) then gives the result that
g0(v;r ') 1 g(v;r r') (31)
For the scattering kernel one has from Eq. (17)
that
P(k,v.v, )dn . Z Pa(k,vi.v)Y b Y b d)
b 0
abb* /yb ,d
a,b
P(kv".v) ,/4f YO* ()
PO(k,v!',v) (32)
which is just the probability per unit time per unit
speed v that a neutron having speed v" will suffer a
scattering collision in region k and emerge from the
collision with speed v. The zero subscript will be
dropped and the function will subsequently be denoted
simply by P(k,v!3.v).
Finally, substituting Eqs. (27), (29), (31),
and (32) into Eq. (26) gives
N(r,v) dr' g(v;r r')S(r',v)
r'
R
+ 4 1 [(V(ktv)+ (k',v)] [V(k,v)+ Y(kv)]
Sf dr/ g(v;rlr)N(r',v)
r'in k
I R
k+ kl f dr' g(v;r r') dv"P(k,v'4v)N(r ,v").
r'in k (33)
11.2 Advantages and Disadvantages of Homogeneous
Green's Function Method
The preceding method of converting the integro
differential transport equation to an integral equation
with a kernel defined by Eq. (10) has at least two
advantages over the more standard methods of solution.
First, as discussed above, the zeroth harmonic approxi
mation involves discarding only products of small terms
as in Eqs. (22) and (23). For example, rather than neglect
ing the higher moments of the neutron density, as must
be done when the integrodifferential form of the transport
equation is used, one need only neglect the products of
the higher moments of the Green's function times the higher
moments of the neutron density.
A second advantage and a unique feature of the
Green's function method described in this work is that the
Green's function defined by Eq. (10) depends only on the
dimensions and geometry of the outer boundaries of the real
problem for which the neutron density is to be calculated,
the total cross section of the medium in the k'th region
of that problem, and the boundary conditions imposed on
Eq. (10). Once the Green's function has been determined,
it may be used to determine the neutron density in any
heterogeneous medium, the only restrictions being that
the outer boundaries of the neutron density problem must
coincide with those of the Green's function problem,
the boundary conditions must be the same for both g
and N, and the total cross section of the k'th region
of the neutron density problem must be the same as
that for the homogeneous medium of the Green's function
problem.
A disadvantage of the method is implied above;
namely, if the outer boundaries and/or boundary conditions
of the neutron density problem do not agree with those for
which the Green's function has been determined, then a
new Green's function must be determined. Another dis
advantage is that, as yet, the method has not been
extended to higher order approximations such as keeping
gm N P, g Pn N terms for n 0, although such an
extension to higher moments would, in general, follow
the work of Wilf (29), and would not alter the above
mentioned advantages.
CHAPTER III
SPEED DEPENDENCE
III.1 General Speed Dependent Equation
Equation (33) is a general speeddependent
equation and is valid over the entire speed spectrum,
the limits of the speed integration being from v" 0
to v" oo. As indicated in the introduction, the
present investigation will be confined to speeds in
the thermal region. In practice, there exist no true
sources for the thermal energy region other than slowing
down sources, and hence S(r',v) 0 in Eq. (33). The
speed integration in Eq. (33) are then separated into
two parts, one for 0 < v i v* and one for v =>v*. The
speed v* represents the upper limit of the speed region
of interest, and is characterized by the following
assumptions (17):
1. For v = v*, the neutron density speed distri
bution is proportional to 1/v2 (that is, the
neutron flux is proportional to 1/E).
2. The probability of a neutron being scattered
from a speed below v* to a speed above v* is
zero.
3. The motion of the collision atoms is ignored
with respect to that of the neutrons with'
speeds above v*.
4. The Maxwellian distribution for thermal
neutrons or collision atoms has negligible
magnitude for speeds above v*.
5. For neutron energies above v*, the energy
transfer is computed assuming that the col
lision atoms are free from chemical bonds.
This implies that v* should correspond to a
neutron speed well above that required to
excite vibrations in the collision molecules,
so that the collision atoms would recoil
freely at the instant of impact.
6. For neutron energies below v*, the energy
transfer kernel P(k,vvv) will depend on the
particular scattering model chosen. A
discussion of particular energy transfer
kernels and their scattering models is given
in references (17, 31).
The S(r',v) term in Eq. (33) is then modified
to include all neutrons scattered into a speed interval
dv about v from speeds above v*. Thus Eq. (33) may be
rewritten as
N(r,v) drI g(v;r r')S(r',v)
r'in k
+ 4 ( [V(kv)+ /'(k',v)] [V(kv)+ Y(kv)]
dr' g(v;r r)N(r',v)
r'in k
V*
+ k dr' g(v;r r') dv" P(k,v"!.v)N(r',v")
r'in k v"0 (34)
where
co
S(r',v) dv" P(k,v!.v)N(r',v") (35)
v"v*
III.2 The Source Term
Assumption (5) of III.1 implies that ordinary
slowing down theory may be used to derive the scattering
kernel for v v* (18, 19). Define
m1
 (36)
m+l
Then ov is the minimum speed that a neutron of initial
speed v can have after collision with a target nucleus
of free atom mass m; the energy transfer probability,
that is, the probability that a neutron of energy E'
will after scattering have energy in the interval dE
about E is (see (18))
T"(r,E..E)dE (dE for a2E' E E'
0 ,for E < a E'
S0 for E E' (37)
or, changing variables,
7l(r,vi.v)dv 7T(rE',E)dE dv
rvdv
so that
7rf(r,v.v) 7r(r,EL.E) dE
1 dE mv
SE' [la2(r)] dv lo2(r)] (mv'
2v
2v
["I, vv for wav' v S v'
S0 for v Oav'
0 for v >v' (38)
Now, the differential scattering cross section,
Lg(r,viav) dv, that is, the probability that a neutron
of speed v' will be scattered at r into dv about v is
equal to the product of the total scattering cross
section at r for neutrons of speed v' times the proba
bility that the speed of the neutron after collision
will be in dv about v. Thus
Zg(r,v4..v)dv s(r,v')7r(r,v'.v)dv. (39)
Then the probability per unit time that a neutron of
speed v* will suffer a scattering collision at r such
that it emerges from the collision with speed in dv
about v is simply
P(r,v'.v)dv v'1%(r,v'Lv)dv (40)
so that for v' > v* one has from Eqs. (38), (39) and
(40) that
P(r,v.v) 2vv's(r,v)
[1 (r ) v'2
2vZs(r,v')
for Ov' s v s v'
[la(r)] v'
0, for v ctv'
(41)
0, for v => v
Assumption (1) of III.1 implies
A(r)
N(r,v') A for v' > v* (42)
vt2
where A(r) is an arbitrary constant with respect to
speed and is normalized to give a stationary distribution
for N(r,v'). Substituting Eqs. (41) and (42) into
Eq. (35) yields
S(r',v) dv" 2vls(k,v") A(r')
S1a2(k) v" v"2
v"'v* [I I
2vs(k)A(r,) v/(k)
n [ldv" (43)
[1a2(k)] j v"3 d
where Zs has been assumed to be independent of v in the
speed range v* 5 v" v v/a(k).
Carrying out the integration in Eq. (43) gives
8(r',v) V A(r') v2 (k) if ca(k)v* v s v.
[1o(k)] V*2 v2
0 if v < a(k)v* (44)
where the constant A(r) determines the source strength;
A(r) will be assumed to equal a spatially independent
constant A, that is, one assumes a spatially flat
slowingdown source distribution. As Honeck (22) has
pointed out, A(r) could, if necessary, be approximated
by the shape of the resonance group neutron density as
predicted by diffusion theory and the use of suitable
resonance group cross sections given in Weinberg and
Wigner (19).
II.3 OneSpeed Equation
To obtain the onespeed equation, one must
assume separability of the spatial and speed dependence
of the neutron density. Then integrating Eq. (33)
over all speeds and using Maxwellian averaged cross
sections, one has
N(r) /dr'g(r r')S(r')
r'
+ 1 [v(k')+ $(k')] [V(k)+ Y(k)) dr' g(r r')N(r')
r'in k
+ 1 V(k) dr g(r r)N(r (45)
r'in k
where use has been made of the fact that
Co
P(k,vu.v)dv V(k,v") (46)
v0
which is just the total probability per unit time that
a neutron in region k and having speed v" will suffer a
scattering collision.
In the onespeed case motion of the collision
atoms may be ignored, so that in conventional notation
Eq. (45) may be written as
N(r) 1 dr' g(r r')S(r')
r'
+ 1 (k') (k)] dr' g(rr')N(r')
r'in k
r'in k
kj g(rfrl')N(r ) (47)
Finally, the last two terms in Eq. (47) may be combined
to obtain
N(r) dr' g(r r')S(r')
+ k t(k')a(k) dr' g(rr')N(r'). (48)
r'in k
CHAPTER IV
DIFFERENCE DENSITY METHOD
IV.1 Development of Difference Density Equations
An alternative to the preceding method of
describing the real neutron density, N, in a hetero
geneous medium is to define a difference density, D,
as the deviation from a known density, N', for a similar
or reference space. Although it may appear at the outset
that such a method is analytically identical to that
developed above for the real density N, nevertheless,
a number of authors have found such a procedure to be
advantageous (13) and, in fact, necessary (12) in order
to obtain results for specific infinite geometry problems.
It will be of considerable interest, therefore,
to develop equations for the difference density, D,
parallel to those above for the real density N. A
direct comparison of the two methods will then be made
in order to determine their relative advantages when
applied to the specific problem to be considered in this
paper.
Thus the difference density D(r,vn) is defined
such that
D(r,vn) N(r,vri) N'(r,vn) (49)
or
N(r,vn) = D(r,vn) + N"(r,vn) (50)
Using primed symbols to denote the reference
space, one can write the transport equation for the
assumed known angular density as
v ~IN'(r,v) + V'(r,v) + (rv) N'(r,Cvn)
.L (51)
S'(r,vn) + dv' dn P'(r,v.v,n4 n)N' (r,v'n').
Substituting Eq. (50) into Eq. (1), subtracting Eq. (51)
from the resulting equation, and grouping all of the
N'(r,vn) terms on the RHS gives
vnVD(r,vn) + [V(r,v)+ (r,v)] D(rvi)
= S(r,vn) S'(r,vn)
v(r_,v)+ (r,v)[,v) (v)+ /(r,v)] N'(r,vr)
+ dv dn' P(r,vv,n'.*0D(r,v'n') (52)
+ dv' d.n P(r,vv,nI)P' (r_,vv,_)] ,n N(rv,',).
As before, Eqs. (3), (4), and (8) are used to separate
the spatial dependence of the terms in Eq. (52). Thus
38
vnVD(r,vn) + [v(k',v)+ (k',v)] D(r,vn)
 S(r,vn) S'(r,vn)
+ [h(rrk,)+h(k'l'r)] [V(k',v)+ K(k',v)] D(r,vn)
rh(kr)h(rklr)] [V(k,v)+ K(k,v)] D(r,vn)
k+k'
 k 1 h(Ekr)(k1^r)] [V(k.v)+ Y(k,v)]
[V'(kv)+ 't(k,v)] N'(r,vr)
+ r h(rkr)h(kl,r)] fdv dn' P(k,vi*v ,nn)D(r,v'nW)
k=1 f
+ h(kr)h('k1r) dv' di' [P(k,vI+v n!*n)
(53)
P' (k,vv,L).)] N' (r,v'n').
The boundary conditions imposed on Eqs. (51) and (53)
are the same as those for Eqs. (1) and (9). Thus the
Green's function which satisfies Eq. (10) may be used
to convert Eq. (53) into an integral equation as before.
Then, using Eq. (12), one may perform the integration
over v' to obtain
D(r,vn.) dr' d' g(v;r,g r' n') S(r' ,v')S' (r' ,vn )
+ [h(r' ,)+h(k, rl') [v(k,v)+ ) k',v)] D(r',vn')
R
 [hr ktr')h(h1r')1 [V(k,v)+ Y(k,v)]D(r',vn' )
k+k'
R
[h(Ekr')h(iklr')] [V(k,v)+ Y(kv)
yv' (k,v)+ 1 (k,v)] N' (r' ,v')
+ kl[h(Ekr')h(rek_r') dv" dr' P(k,v" vn'n*')D(r',v'tn")
+ k E [h(Ekr' )h(Ekir') dvt" d [P(kv' mt.')
i (54)
P' (k,vvn'".)] N' (r' vn") .
The r' integration for Eq. (54) is identical
to that for Eq. (13) and the result is
D(rv_) dr' dn' g(v;r,n r' ,n)[S(r',v f)S'(r_',vn)]
+ V(k',v)+ (k',v)] [V(k,v)+ Yk,v)
dr' da g(v;rn ir'L l )D(r' ,v')
r'1in k 1'
R
+ dr' dn' g(v;r,nlr'l,.')
r'in k n
dv" a" P(k,v" vn'1,c n')D(r',v"n')
 [V(k,v)+ Y(k,v)] [V'(k,v)+ (k,v)
S dr' dn' g(v;rnl r',n')N'(r',vn!')
rrin k '
R
+ dr, dn' g(v;r,nlr',n') dv"fdrl" P(kv4I.vI.nIn')
r'in k T.'
(55)
P' (k,v".v,n ".')] N' (r' ,vI'") .
The angular integration presents no new
problems and may be performed as in Chapter II. Making
the zeroth harmonic approximation and asking only for
the scalar neutron difference density, one obtains
D(rv) d 1 dri g(v;rjr9)[ (r',v)S'(r',v)]
rt
+ ([V(k',v)+ (k',v)]
47r ki
[V(k,v)+ Y(k,v)]) dr' g(v;r r')D(r' ,v)
r'in k
R
k dr' g(v;rlr') dv" P(k,vU*v)D(r',v")
klr 
r'in k
R
ki
+ 1
471
4"/
Sdr' g(v;rr')N'(r',v)
r'in k
dr' g(v; rr')
r'in k
(56)
*dv" P(k,vt v)P' (k,v"v) N' (r' ,v").
Writing Eq. (56) specifically for the thermal
spectrum yields
D(rv) dr' g(v;rr')[S(r'v)S'(r',v)]
r'in k
R
[V(k',v)+ (k',v)]
[V(k,v)+ /(k,v)] dr' g(v;r r')D(r',v)
r'in k
S[vk,+ (k,v) )]
 [V'(k,v)+ Y'I(kv))
1
+1
I R
+
kw,
v f
V*
1 R
4+ F i dr' g(v;r r') dv" P(k,v"!v)D(r',v")
r'in k v"'0
1 R i (k,v)+ (kv)]
[V'(k,v)+ '(kv dr' g(v;r r')N'(r,v)
r'in k
R
+ k~ dr' g(v;rlr')
r'in k
(57)
v*
dv" [P(k,v,.u.v)P' (k,v".v)] N' (r' ,v")
VVO
where, in region k,
ro
S(r',v) S'(r',v) dv" P(k,v!.v)D(r',v")
vlv*
00 (58)
+ dv" [P(k,~,!v)P'(k,v".v)] N'(r',v")
vv*
and the scattering kernels for v > v* are determined
from ordinary slowing down theory as in 111.2.
Finally, the one speed version of Eq. (56) is
D(r) dr' g(r r')S(') (
r'
1 R
+  k
+ Rl
1 R
+1 R
47( N1
v[zt(k') It(k)] dr' g(r r')D(r')
r'in k
v Zg(k) dr' g(r r')D(r')
r'in k
v [zt(k) Z((k) drI g(r rf)N'(r')
r'in k
v[(k) (k) dr' g(r r')N'(r')
r'in k
or if various terms are combined,
D(r) /dr' g(r r')[S(r')S'(r')]
r'
(59)
+ v k') Zak)] dr' g(r r)(r
r'in k
S1 v[Za(k) (k)] dr' g(rr')N'(r') (60)
47r T l k f i
r'in k
IV.2 Equivalence of Equations for N and D
In order to establish the equivalence between
Eqs. (34) and (57) and between Eqs. (48) and (60),
the onespeed equations will be considered first.
Subtracting Eq. (60) from Eq. (48) and using
Eq. (50) gives
N'(r) 1 dr' g(r r')S'(r')
rt
+ k V[Zt(k*) a(k)] dr' g(r )N(r')
r'in k
+ 1 v[Za(k) Z(k)] f drv g(rlr')N'(rt)
k 47j
r'in k
_ fdr' g(r r')S'(r')
47r f
rt
+ RF v[ l(k) .(k)] drt g(rr)N(r'). (61)
r'in k
Comparing Eqs. (61) and (48), one observes that the
requirement that Eqs. (60) and (48) be equivalent is
simply that It(k') = (k'), that is, the total cross
section for the k'th region of the reference space be
identical to the total cross section of the k'th region
of the real space of interest. This simply says that
the Green's functions in Eqs. (48) and (60) must be
identical, that is, the Green's functions in Eqs, (48)
and (60) must both satisfy Eq. (10) and its boundary
conditions. But this is just the condition imposed in
the original derivation, and hence Zt(k') is equal to
Z (k'); equations (48) and (60) are therefore proved
to be equivalent formulations of the same problem.
Similarly, subtracting Eq. (57) from Eq. (34)
and using Eq. (50) gives
N'(r,v) = dr' g(v;rr')' (r ,)
*E klf ~ # ~
r'in k
1 R
47T kl
[V(k,v)+ a(k' ,v)]
Iv(k,v)+ Y(k,v)] dr' g(v;r r')Nt(r',v)
r'in k
V*
+ k dr' g(v;r r') dv" P(k,v"..v)N (r' ,v")
r'in k vIO
1
+
k [V(k,v)+ ~(kv)
[V'(k,v)+ f'(k,v)]
r'in g(v;rlr)N'(r',v)
r'in k
Sdr' g(v;rlr')
r'in k
* dv" [P(k,vt..v)P' (k,v"!.v) N' (r' ,v")
v" 0
1'
R
EI
r
I
r'in k
dr' g(v;r r')8'(rt,v)
+ 1 [V(k'iv)+ K(k',v)]
['(k,v)+ (k.v)) dr' g(v;rjr')N'(r',v)
r'in k
v*
S R
+ 1 dr' g(v;r r') dv" P'(k,,v'la)N'(r'v,v")
r'in k V1=0
Comparing Eqs. (34) and (62) one again observes
that, just as in the onespeed equation, the requirement
that Eqs. (34) and (57) be equivalent is simply that
[v(k',v)+ (k' ,v)] be equal to [V'(k',v)+ ('(k',v)],
that is, the medium of the k'th region of the reference
space must be the same medium that is in the k'th
region of the real space of interest. As has been
mentioned above, this condition was imposed in the
original derivation so that the equivalence of Eqs. (34)
and (57) is established. A discussion of the relative
usefulness of the two formulations will be deferred
until the general equations have been specialized to
slab geometry.
CHAPTER V
PLANE SLAB GEOMETRY
V.1 OneDimensional Problem
The equations of the previous chapters were
written for a general threedimensional geometry, and
the Green's function g(v;r r') was interpreted as the
number of firstflight neutrons per unit element of
volume at r having speeds v due to an isotropic point
source at r' emitting 47r neutrons per unit time having
speeds v in a homogeneous medium. Before specializing
Eqs. (34) and (57) to an infinite slab geometry it will
be convenient to define a new Green's function G(v;r r'),
such that
G(v;r g(v;rr') (63)
Thus G(v;r r') is the number density of firstflight
neutrons at r having speed v due to a unit isotropic
point source at r' emitting neutrons of speed v in a
homogeneous medium. Rewriting Eq. (34) in terms of
G(v;rlr') gives
r'in k
+ 1 fV(k',v)+ /(k: ,v) [V(kv)+ /(k,v)]
*/ dr' G(v;rir')N(r' ,v)
r'in k
v* (64)
+ 17 dr' G(v;r r') dv" P(k,vlvm)N(r',v")
k.f I f
r'in k v"0O
and similarly for Eq. (57).
The transformation to plane geometry may be
performed as follows. Referring to Figure 2, the element
of volume dr' may be expressed as
dr' 2 7Ta da dx'
Define
R r r'
R r r'I (65)
From Figure 2, and Eq. (65) one obtains
R ( xx' 2 + a2)1/2 .x_
In I
r fixed field point
r' variable source point
Y cos G
z /
//'r O
/ X
r x
/ L._
da
dx'
Figure 2 Plane Geometry Transformation
so that for constant Ixx'l
dR = a da Ixx'
a_2 dr .
R 2
Hence
dr' 2 7 R dR dx'
_Hxx'lX 2 d_
S2[Ixx'I dx' (66)
where the i integration is from 1l to y0, and the R
integration is from R jxx'I to R oo .
Observing that in plane geometry the r de
pendence of the functions N, N', D, S, and P is simply
a spatial dependence on x, one may write Eq. (64) for
plane geometry as
R
N(x,v) 1 dx' S(x',v) G(v;r r') 2 7r R dR
x'in k R xx'I
+ 1 [v(k v)+ '(k',v)] (k,v)+ (kv)]
Sdx' N(x',v) G(v;r r') 277R dR
x'in k Rxx'l
v* 00
+ 1 dxt dv" P(k,vl.v)N(x',v") G(v;r r') 27TR dR
xin k v"O RIxx'i (67)
Define
G(v;xx') G(v;rr ') 2rrR dR (68)
R xx'l
The function G(v;x x') is the density of firstflight
neutrons at x having speeds v due to a unit isotropic
plane sourceat x' emitting one neutron per unit time
with speed v in a homogeneous medium. This function is
well known for an infinite medium (see Chapter 2 of
Case, et al. (21)) and will be discussed extensively
for a finite unit cell in Chapter VII.
Using Eq. (68) to rewrite Eq. (67) gives
R
N(x,v) 4 /dx' G(v;x ,x)S(x,,v)
x'in k
+ k1 [V(k',v)+ /(k',v)][V(kv)+ Y(k,v)]
dx, G(v;x x')N(x',v)
x'in k
V*
+ dx' G(v;x x')
x'in k v"0
Similarly, Eq. (57) for plane geometry becomes
D(xv) / dx' G(v;x x')B(x'.v)S'(x'.v)]
x'in k
+ k1 [(k',v)+ (k',v)][V(kv)+ (k,v)]
/dx' G(v;x x')D(x',v)
x'in k
+ f dx' G(v;x x')
x'in k v"0
R
kx
3v" P(k,vL.v)D(x',v")
[V(k,v)+ Y(k,v)] [V (k,v)+ YO (k,v)]
dx' G(v;x x')N'(x',v)
x'in k
(70)
v*
+ dx' G(v;x x') dv" P(k,v'wv)P'(k,vtlv) N' (x',v").
x'in k v0
dv" P(k,v"L.v)N(x',v") (69)
V.2 Symmetry Considerations
For subsequent work the space of interest
will be assumed to be a plane geometry symmetrical
about the plane x0. In the summations over k as
previously discussed, kl denoted the leftmost region
in space. The regions will now be counted differently,
letting kl denote the region adjacent to, and to the
right of, the x0 plane. Thus a summation over k from
kl to kR is replaced by a summation over k from
k R/2 to k R/2. Let
R2 " (71)
denote the number of regions in the halfspace. For a
general function f(x) which is symmetric about x0,
that is, f(x) f(x), one can then write that
R
E f(k)J dx' G(v;x x')N(x',v)
kl
x'in k
R2
R f(k) dx' G(v;xlx')N(x',v)
kR f(k)
x'in k
+ ( f) dx' G(v;x x')N(x',v)
kR2 k f(k)
x'in k
0 b
= dx' + dx f(x')G(v;xlx')N(x',v) (72)
x'b x'O
where b and b denote the right and left boundaries of
the space of interest. The notation used is defined
such that
cm c2 7C
+f ... + + k (73)
kl kn kc1 ken
and
(74)
C2 cm c2 cm
dx' + ... + dx' f(x') dx'f(x') + ...+ dx'f(x').
x'c1 X' cn x'c 1 X'cn
In the first integral in Eq. (72) let x'y,
dx'=dy. Observing that f(x)f(x) and N(x,v)N(x,v),
one can write
0 0
dx'f(x')G(v;x x')N(x',v) dy f(y)G(v;x y)N(y,v)
x' b yb
b
/dy f(y)G(v;xly)N(y,v) (75)
yZ0
Then let yx' so that Eq. (75) becomes
0 b
Sdx'f(x')G(v;xx')N(x',v) dx'f(x')G(v'x x')N(x',v).
xb 0 (76)
Substituting Eq. (76) into Eq. (72) yields
R
Z f(k) dx'G(v;x x')N(x',v)
fk
x'in k
b
dx'f(x') [G(v;x x')+G(v;x x')]N(x',V)
x'0
(77)
f (k) dx' [G(v;x x')+G(v;x x')]N(x',v)
k1n
x'in k
As noted above, G(v;x x') is the scalar
density of firstflight neutrons at x due to an isotropic
unit plane source at x'. If one defines
H(v;x x') G(v;x x') + G(v;x x') (78)
to be the halfspace Green's function, that is, the
scalar density of firstflight neutrons at x due to
isotropic unit plane sources at both x' and x', then
equation (69) may be rewritten for the halfspace in
the form
N(x,v) dx' H(v;x x')S(x',v)
x'in k
+ kV (V(k',v)+ Y(k',v) fV(k,v)+ 7(k,v)1
f dx' H(v;x x')N(x',v)
x in k
v* (79)
+ dx H(v;xlx') dv" P(k,v'.v)N(x',v")
x'in k v"0
while Eq. (70) becomes
D(x,v) R 2dx' H(v;x x') S(x',v)S'(x',v)]
x'in k
+ 2 ([V(k'.v)+ (k',v)][V(kv)+ (kv)]
dx' H(v}x x')D(x',v)
x'in k
v*
+ T dx' H(v;xlx') dv" P(k,v!".v)D(x',v")
x'in k v"0
 ([V(k,v)+ /(k,v)][V'(k,v)+ Y'(k,v)]
dx' H(v;xIx')N'(x',v)
x'in k
(80)
v*
+ f dx' H(v;x x') dv" [P(k,v!.Lv)P'(k,v!..v) N'(x'f,v")
x'in k v"10
The onespeed equations may also be written
for the onedimensional halfspace. Thus Eq. (48)
becomes
N(x) dx' H(x x,)S(x')
x >0
+ :2 v[Zt(k') Ia(k)] fdx' H(x x')N(x') (81)
x'in k
while Eq. (60) becomes
D(x) .dx H(x xI) [(x')BS(xv)
x') O
59
+ E2 v[t(k') Za(k)] dx' H(xx')D(x')
x'in k
i2 v (k)a( j) dx' H(xJx')N'(x') (82)
x'in k
CHAPTER VI
COMPARISON OF N AND D METHODS
FOR SIAB GEOMETRY
For the direct comparison of the formulation
in terms of real density,N, with the formulation in
terms of difference density, D, only the onespeed
equations will be considered, namely Eqs. (81) and (82).
Similar arguments would apply to the multispeed equations.
Suppose that all R regions of the reference
space contain the same medium and that the reference
space itself constitutes one spatial period of an
infinitely repeating lattice, that is, the reference
space is a unit cell. Then the neutron density N'(x)
in the reference space is that for an infinite homo
geneous medium and may be assumed spatially flat.
Since the neutron density and medium are both
spatially uniform in the reference cell, one may assume
that the statement of neutron conservation holds not
only over the entire reference space, but also over any
arbitrarily small element of volume in that space. Thus
at any point
B'(x) S' v IN'. (83)
Substituting Eq. (83) into Eq. (82) and
cancelling terms where possible, one obtains
D(x) dx H(xlx')S(x')
f',0
XI'0
+ a vt(k') (k)] dx' H(xjx')D(x')
x'in k
v a(k)N' dx' H(x x') (84)
x'in k
Then using Eq. (50), one may combine terms to write
D(x) dx, H(x x')S(x')
X'O
+ k v 2(k') f dx' H(x x')D(x')
v Z.a(k) dx' H(x lx)N(x') (85)
k1 f
x'in k
Adding and subtracting the quantity
v t(k') fdx' H(x x')N(x')
k1in k
x'in k
to Eq. (85) gives
D(x) dx' H(x x)S(x')
x 'O
+ 1k2 v[Zt(kt) Za(k) f dx' H(x x')N(x')
k1
x'in k
R2
I v lt(k') dx H(x x')N'
kin
x'in k
dxt H(x x')S(x')
x*=0
+ 2 v[t(k') Za(k)] dx' H(x ')N(x')
x'in k
v I(k')N' dx' H(x x') (86)
x '0
But from Eq. (10)
v Zt(k')H(x x')dx' 1 (87)
x 50
so that using Eqs. (50) and (87), one can reduce Eq. (86)
directly to Eq. (81).
The direct reduction of Eq. (82) to Eq. (81)
having been established, the question remains whether
the use of Eq. (82) has any advantages over the use of
Eq. (81). It will be of particular interest to consider
whether or not the complexity of Eq. (82) can be reduced
by relating the real source S to the source 8' in the
reference space. A few such special cases are considered
below. In each case the assumption is made, as it was
above, that all R regions of the reference space contain
the same medium and that 8' and N' are spatially flat.
Case 1 S(x) 8' Over All Space
Equation (82) then becomes
D(x) 2 v I(k') a(k)] dx' H(xJx')D(x')
x'in k
va(k) i(k)]N f dx' H(xxt)
R 2 tk' a(k)] dxz H(x x')D(x')
x'in k
v Za(k)N' fdx H(x x') + vE(k'
x'in k
where Eqs. (83) and (87) have been used to obtain the
last term of Eq. (88).
Since S(x) = S' implies that S(x) is spatially
independent, Eq. (87) may be used to show that the first
term in Eq. (81) becomes
dx' H(x xt')(xt) = 8 dx' H(x x')
x'=O x '=O
S'
v St(k') (89)
Thus Eq. (88) can be reduced directly to Eq. (81) in
the same manner that Eq. (82) was reduced to Eq. (81),
so that there is no advantage of the difference density
formulation over that for the real density N. The
reason for this is that, although it is given that
S(x) S S', one must still calculate N' based on 8',
and the N' terms then act as source terms in the differ
ence density equations.
Case 2 S(x') 0 For Region k1, S(x') 8'
For Regions k 1
Equation (82) becomes
D(x) dx' H(xIx')S'(x')
x'in 1
+ 2
ki
 2
kl
v[ZCk) (k)]
dx' H(xIx')D(x')
x'in k
Sdx' H(x x')N'(x)
:'in k
To specialize further, let R2 2 and k' 2.
Eq. (90) becomes
D(x)  S' dx'
x'in 1
+ v[Zt(2)~a(1)
H(x x')
dx' H(x x')D(x')
x'in 1
+ v Z~(2) dx' H(x x')D(x')
x'in 2
v[a(1) a(1)] N' dx' H(xl x').
x'in 1
Using Eq. (83) one may write Eq. (91) as
D(x) v[Zt(2), a(l)] idx' H(xI xtD(x')
x'in 1
(90)
Then
(91)
+ v zS(2) dx' H(x Ix)D(x')
x'in 2
a(l) N' dxv H(xlx') (92)
x'in 1
or, again substituting Eq. (83) into the last term of
Eq. (92), one finally obtains
D(x) v[Zt(2) Za()] f dx' H(xx')D(x')
x'in 1
+ v ZS(2) dx' H(x x')D(x')
x'in 2
(1a(') 8'
7 a(2l) dx' H(x x') (93)
x'in 1
For this special case, Eq. (81) becomes
N(x) SI dx' H(xlx')
x'in 2
+ v[t(2)Za(1)] /dx' H(xlx)N(xt)
x'in 1
+ v Zs(2) dx' H(xlx')N(x') (94)
x'in 2
The obvious difference between Eqs. (93) and
(94) is the source term. In a typical application of
calculating the scalar neutron density in a space
containing a gold foil in a graphite medium, the ratio
Za(1)/Za(2) would become
(1) Au 5.79
1 _7 1.94 x 104
a a
a(2) ^C .000299
It can be seen that in a finite medium containing a
highly absorbing region, even though it be of small
spatial extent, the neutron density will not necessarily
be a small deviation from a known density in a reference
space with no absorber. In fact, the deviation in
magnitude of the densities will be quite large, even
though the spatial shape of the neutron density N may
not be vastly different from the spatial shape of the
known density N'.
It seems obvious now that for problems involv
ing finite geometries the analytic method of solving for
a difference density, D, offers no advantages over
68
the direct solution for the unknown density N; the
difference density method may, indeed, be limited in
usefulness to the solution of infinite geometry problems.
CHAPTER VII
THE FIRSTFLIGHT GREEN'S FUNCTION
FOR A HOMOGENEOUS MEDIUM
VII.1 Differential Equation and Boundary Conditions
The firstflight Green's function
6k' (r,v r',v'n') for a general threedimensional homo
geneous medium is defined such that it satisfies Eq. (10),
namely
va* Vgk' (r, vn r', v'n' )
+ [V(k',v)+ (k' ,v)] gk'(rvnr' ,v'n')
&(rr )S(vv )S(cnn'); (95)
further, gk (r,vnr',v'_1') is subject to the same
boundary conditions that are applied to the angular
density N(r,vn).
Considering the RHS of Eq. (9) as a fictitious
source Q(r,vn), one may rewrite Eq. (9) as
vn.7VN(r,v) + [V(k',v)+ ((k',v)]N(r,v)) Q(r,vn) (96)
Thus N(r,vn) represents the firstflight angular density
in a medium as discussed in Chapter I, the removal pro
bability being equal to the total probability of col
lision in region k', with a source of neutrons Q(r,vn).
If the homogeneous medium k' is of infinite
extent, then the solution to the onespeed version of
Eq. (95) is known to be (see Chapter 2, of Case, et al.
(21))
t(k') d
gi (rlr S 2 nr!' )'(n fr') (97)
where the superscript oo on g signifies that Eq. (97) is
valid only for an infinite medium or for a finite medium
enclosed by a nonreentrant surface with vacuum boundary
conditions, and rr' denotes the unit vector (rr')/jrr'
The extension to the speed dependent problem is obvious
from the discussion of Chapter I and Eq. (12). Thus
Zt(k',v) rr'
v
S0, (r ', ,v'n') e
(98)
o(vv')&(nn')An ;Cc') .
Integrating Eq. (98) over v' and n! gives
t(k'v)lrr'
gg?(v;r,nr') e t(k ,v ); (99)
v rr (99)
integrating over n then yields
k (v;r r') e Ztk v) (0oo
v rr 2 (100)
where use has been made of Eq. (30).
The simple form of Eq.' (100) is due to the
fact that the medium k' is homogeneous so that the
collision probability is not a function of the flight
path of the neutron. As pointed out by Case, et al.
(21), the case of a spatially variable collision proba
bility is a considerable complication because absolute
coordinates instead of only relative coordinates between
the field and source points enter the problem. It might
be well to point out again that a unique advantage of
the homogeneous Green's function method being described
is that it completely avoids the complication of a
spatially variable collision probability in the evalu
ation of the Green's function.
From Eqs. (63) and (100) one can write the
expression for G(v;r r') for an infinite medium. Thus
(k';v _rr'1
G (v;rr') e(101)
47 v rIr' 2
From Eq. (68) one can then evaluate G(v;x x') for an
infinite medium. Thus, using the transformations of
V.1, one obtains
,3 ,,v Ir,rIt(k'
Gco (v;x x') 27r R dR
47T v r1r' 2
RIxx'I
S e t(k',v)R 1 dR (102)
RI xx'
The exponential integrals or En functions (21, 24) are
defined such that
En(x) exu un du (103a)
1
i pn2 eX/P dp (103b)
00
Sxn1 eu un du (103c)
ux
so that Eq. (102) may be written
Goo (v;x x') L El [Zt(kv)Xxx] (104)
The Green's functions discussed above were,
in general, for an infinite medium. The problem
eventually to be considered, however, is that of a unit
cell; that is, the space of interest represents one
unit volume of an infinitely repeating set of unit volumes.
The application of the unit cell approach is
based on the fact that most heterogeneous nuclear reactors
have fuel and moderator arranged in a repeating lattice
structure. For calculational purposes a lattice is sub
divided into a number of identical unit cells (18, 19).
The spatial and speed dependent neutron density, usually
just the scalar density, is calculated for the unit cell
and is subject to particular boundary conditions. Quanti
ties such as the thermal utilization are then calculated
for the unit cell and are assumed to be valid throughout
the lattice. The reactor is then treated as a homo
geneous reactor having the same thermal utilization as
the unit cell, and the scalar neutron density is calcu
lated for the entire reactor.
In most reactors the lattice is very large
relative to the dimensions of the unit cell, so that
it is not unreasonable to treat the unit cell as if it
were in an infinite lattice. For a complete discussion
of this assumption, see Chapters 7 and 18 of Weinberg
and Wigner (19). If the assumption is then made that
the lattice is infinitely repeating, then at steady
state the flow of neutrons in dn about direction ni
across a point on the boundary between adjacent unit
cells is equal to the flow in d_ about direction i
where L2 is the mirror image of nl as illustrated in
Figure 3. This is physically equivalent to the
assumption that the cell boundaries are perfectly
reflecting. Hence in the analysis that follows,
neutrons which actually arrive at r in the unit cell
by crossing the cell boundaries on firstflights can
be considered to have originated within the unit cell
at r' and suffered reflections at the cell boundaries
before arriving at r.
S2 180    
T2 T 1
Figure 3 Symmetry of Neutron Flow at the
Boundaries of Unit Cells in An
problem in the following manner. Due to reflections at
the cell boundaries, firstflight neutrons from an
isotropic point source at r' can arrive at r by various
paths, so that the total density at r is the sum over
all possible paths i, of the densities of neutrons
arriving at r by path i and having traveled a total
path length Li.
Alternatively, the neutrons at r can be con
sidered to have arrived by firstflights from i image
sources at r located a straightline distance L from
soure i at located a strahtine distance L from
the field point r. The total angular density at r is
then the sum over the i image sources of the contri
butions to the density at r from each image source at
rt. (See Chapter 7 of Morse and Feshbach (16) for a
discussion of the method of images.)
Thus for a general homogeneous medium, one
has from Eqs. (63) and (99) that
t(k',v)L
G(v;r,n r') e (t( rv)i)
i 4 TTv2 i (105)
where Li is the total path length traveled by a neutron
arriving at r along a firstflight path from an image
source at ri. From the alternate point of view, Li is
the path length traveled by a neutron along path i as
it suffered reflections at the cell boundaries before
arriving at r.
Integration of Eq. (105) over n yields
" t(k',v)Li
G(v;r r') 4 7tvvJLi (106)
O vLi2 (106)
4 7r vLi1
As one would expect, for an infinite medium or for a
finite medium bounded bya nonreentrant surface with
vacuum boundary conditions Eq. (106) reduces to Eq. (101);
for a finite medium bounded by a reentrant surface with
vacuum boundary conditions, Eq. (106) reduces to the
infinite medium Green's function with spatially variable
collision probability as given in Chapter 2 of Case, et al.
(21), namely
R (r sR/R, v) ds
G(v;r r') e 2(107)
4 7 vR2
where Zt(s,v) 0 outside the bounding surface.
VII.2 Analytic Form of G(v;xlx')
The physical effect of introducing reflecting
boundaries is to increase the density of neutrons at a
point due to reflections at the boundaries which, in
effect, send neutrons back into the finite unit cell.
The total neutron density at r may be separated into
two parts; namely, that due to firstflights from the
source at r' within the unit cell and that due to
neutrons that may be considered either as arriving
from i image sources located at r' outside the unit
cell or as having suffered one or more reflections at
the boundaries after emission at r'.
The firstflight angular density at r due to
a unit isotropic point source at r' in a homogeneous
medium is simply the infinite medium Green's function
obtained from Eqs. (63) and (99), namely,
0,,nv) rr,e
oG (v;r,r,) e (a ^ ').
47T v rr  (108)
Let Gfc(v;r,n r') denote the angular density at
r of neutrons which do not arrive by firstflights
from the unit isotropic point source at r' but instead
arrive at r only after having suffered reflections at
the mirror boundaries of the cell, or alternatively,
arrive by firstflights from the i image sources
located at ri. Then let Li,n denote the total path
length traveled by a neutron along a path having a
direction that is characterized by the unit vector
ri' pointing from the image source at r' to the
j ji
field point r, the neutron arriving at the field point
r after having suffered n reflections at the mirror
boundaries of the cell. One can then write from Eq. (105)
that
Gfc(v;rnE I) z e lt(k',v)Li,n _
nl i 4 VLln 2 
n (109)
The double summation over n and i is, of course, equiva
lent to a single summation over all the image sources.
The total Green's function for the unit cell is then
G(v;r,n Pr') G (v;r,nlr') + Gfc(v'r,rl r') (110)
Thus it is seen that Gfc represents a finite medium
correction term to be added to GO.
In V.1 the plane geometry equations were
obtained by integrating the point source kernel over
all the sources on the yzplane at x'. For a unit cell
with plane reflecting boundaries, a neutron which
arrives at r after having suffered n reflections, and
thus having traveled a path length Lin, can be con
sidered to have originated at a fictitious or image
plane source located at a distance pLi,n from the field
plane. The cell boundary at which the neutron made its
last or nth reflection lies between the field plane and
the image source plane.
The transformation from the point source
Green's function to the plane source Green's function,
that is Eq. (68), must be modified using the above
method of images to locate image plane sources at
distances pLi,n from the field plane. Referring to
Figure 2, G(v;x,plx') is the number of neutrons with
direction cosines between p and + du crossing a ring
element of area 27 Tper unit time, the ring element
being located in the plane at xt. If f is the azimuthal
angle indicated in Figure 2, then
dr du dP (111)
Hence, using the transformations of V.1,
Goo (v;x, Ix')dx dp dx' ,
2 0
x d d G I
fma0 f2l1 (112)
Since the sources are assumed to be symmetri
cal about the xaxis, G(v;r,flr') is a function only of
r,r', and u. Thus the integration over c may be
carried out in Eq. (112):
1
GC (v;x,p x') 47r2/ G(v;rlr') X 2 (113)
P O =0
For slab geometry G(v;r,nlr') has a delta function
behavior with respect to p' so that Eq. (108) becomes
Gc (v;rn.r') 
e t(k',v) (xx)/p
e r r .
4x (114)
47yv(Ixx'l /p)2 27r
The p'integration may then be carried out in Eq. (113).
The result is
Gc(v;x,pIx') 
 (k',v) rr
2v
or using Eq. (2),
G C(v;x,p x') 
Z'(k'.V) I (XX'f)^l
2v II
Zt(k,v) (xx')/P
h(xx') for f > 0
(115)
h(x'x) for < 0
2v I
where the step function has been introduced to account
for the fact that in an infinite firstflight medium,
the flow of neutrons in the positive P direction is
zero unless x is greater than x', while the flow of
neutrons in the negative p direction is zero unless x
is less than x'. For a more complete discussion see
Chapters 2 and 4 of Case, et al. (21).
In the same manner, one obtains the correspond
ing expression for the finite medium correction term.
Thus
Zt(k',v) (xxi)/p
GfC(v;x,pl x) e h(xxi) for p 0
1 2v 
(116)
g(k v) I(xx) I
Ee , / h(xix) for p 0
i 2v i
where the single summation represents a sum over all i
image source planes located respectively at xi. The
first two image source planes are illustrated in
Figure 4.
Substituting Eqs. (115) and (116) into
Eq. (110), one can verify that each term in G(v;x,p x')
satisfies the onedimensional form of the sourcefree
transport equation corresponding to Eq. (95), namely
oG(v;x,p x')
VU xy + v Zt(k',v) G(v;x,Ix') 0 (117)
where v Zt(k',v) is equivalent to [V(k',v)+ (k' ,v) .
Further, G(v;x,p x') satisfies the reciprocity theorem
(21), that is,
G(v;x,p x') G(v;x' x,) (118)
Finally, each term of Eqs. (115) and (116) satisfies
the symmetry condition at the boundary x +b that the
flow of neutrons in dp about p is equal to the flow in
dp about ', that is,
G(v;x,u x') xl b G(v;x, ix') .xb (119)
x+b x+b
Since each term of the Green's function G(v;x,ujx')
satisfies the original boundary conditions, it can be
concluded that the function itself satisfies the boundary
conditions. One can verify that, in agreement with
sections 4 and 5 of Case, et al. (21), a source plane
of strength qs at x' is equivalent to a discontinuity
of the normal component of the angular current, where
for an isotropic plane source, qs 1/2. The angular
current is simply the number of neutrons crossing
unit area perpendicular to the direction of flow per
unit time and unit solid angle. In particular, for
the infinite medium terms one can write
vpG(v;xuIx') xx'+ vuG(v;xGiux') x x' q5s() "
In an infinite medium the second term on the LHS
vanishes for p >0 while from Eq. (115) the first
term is equal to 1/2, that is, the source plane value.
In the solution for the scalar neutron density
in Eqs. (69) and (77) one needs the corresponding
scalar Green's function G(v;x x') obtained by integrat
ing over all p the function G(v;x,p x'). Thus,
1
G(v;x x') G(v;x,p x')dp (120)
P1
Recall that G(v;x x') is the firstflight
density of neutrons at x due to a unit isotropic plane
source at x' in a homogeneous medium. Then, referring
to Figure 4, one can see that for image sources to the
right of x, that is, x => x, the neutrons reaching x
must have a negative direction cosine u, while for
image sources to the left of x, that is, xi <= x, the
neutrons reaching x must have a positive direction
cosine. Thus the contribution of the image sources is
0 1
Gfc(vI;x x) = Gfc(v;x,plx.) dp + Gfe(v;x,,Ix')di (121)
0 1 0
where the first and second integrals represent the
contributions of image sources placed respectively to
the right and left of x.
c7G
x2bx'
Image Source
Plane at
xb(b+x')
Figure 4
xb
xO
01800
'
1,1
xb x2bx'
Image Source
Plane at
xb+(bx*)
Slab Unit Cell with Mirror
Boundaries at x+b
Letting Gtc denote one term of Eq. (116),
one may write the first integral in Eq. (121) as
0 0
f fe(v e t(k'tv) I(xX)/IU
Gi (v;x, ')du d
iI2vl I2
e Zt(k',v)(xxi)/l d (122)
O0
the last integral having been obtained from the preced
ing one by a successive substitution of variables, first
letting yu, and subsequently letting iy.
In general then, one can sum over all image
sources and obtain the result that
G (v;x x') f e t(k',v) (xxi) d (123)
Again, for a firstflight medium the range
of integration over ) for the infinite medium contri
bution is from 0O to P1. Thus one can write
1
G(v;xlx') f e t(k',v) (xx')/ d
O0
1
+ t(ck',v) (xx') 4)
i 2v (124)
P0
or from Eq. (103),
G(v;x x') El ( t(k'v)xx'j]
(125)
+ E El[Zt(k',v) xxil] )
To express G(v;xIx') as a function of the
unit cell dimensions it will be instructive to return
to the physical concept of reflection of neutrons in
order to locate the image sources with respect to the
cell boundaries. Using the notation for path length
previously discussed with regard to Eq. (109), one can
write the expression for GfC(v;x x') in Eq. (125) as
Gc(v;x x') E E El[7t(k',v)ILi.nl] (126)
nl
where PLi,n is, as before, the distance separating the
field plane from an image source plane, the image
plane corresponding to a neutron having made n re
flections at the mirror boundaries.
In Eq. (126), the summation over i extends
over all possible paths by which a neutron having its
direction cosine equal to u can arrive at the field
point after suffering n and only n reflections at the
boundaries. In plane geometry, it is easily seen that
for a firstflight medium there corresponds to any given
p only two such paths for any n. Hence the summation
in Eq. (126) runs from i1 to i2. Figure 4 illustrates
the two possible paths, L1,1 and L2,1 corresponding to
S= 0 and u < 0 respectively, by which a neutron can
arrive at x after suffering only one reflection. The
unit cell in Figure 4 has its midplane at xO and the
reflecting boundaries at x b. The cell width,
that is, the distance between the two boundaries is 2b.
The expressions for pLi,n can be obtained by
simply adding the total path length a neutron travels
while making n reflections. The expressions for pi,n
depend on the orientation of the source plane and field
plane x relative to each other and to the x0 plane.
It has been shown in V.2, however, that one need only
consider field planes for positive x, because the
solution of the halfcell is symmetrical with respect
to x.
The expressions for aLin are derived in
Appendix A where it is shown that Eq. (125) may be
written, for x > O, as
G(v;xx') 1. E (k',v)Ixx']
+ 1 EI[Zt(k',v) 2nbx'+(1)nIxl
2v n1
(127)
+ E[Zt(,v)12nb+x'(1)nx]) for x = 0.
Note that Eq. (127) can be shown to be term by term
equivalent to a more cumbersome series representation
of G(v;x x') derived by Aswad (25).
Finally, the halfcell Green's function
defined by Eq. (78) becomes
H(v;xlx') = E[l t(kO,v)xx'I] + E1[t(k',v)x+x'I
1 co
+ I E 1I[Zt(kv) 12"b'x'l)n]
n=1
+ EL[Zt(k',v) 2nb.+x(1)nxi]
+ E,[Zt(k',v) 2nb+x'+(l)nx]
+ El[Zt(k',v)12nbx'(l)nxi] )(128)
for x =0 and x' O 0
VII.3 Monte Carlo Generation of G(v;xlx')
For complex finite geometries, the complete
analytic form of the firstflight Green's function may
be exceedingly difficult, if not impossible, to derive.
In such cases, it may be advantageous and even necessary
to represent the physical or mathematical system by a
sampling operation satisfying the same probability laws.
Such a process has come to be called a Monte Carlo
method (26).
Essentially, Monte Carlo methods are "paper
experiments." The experiment consists of performing
specific sampling operations many times, the sampling
operations satisfying the criteria mentioned above. As
in many other experimental procedures, the geometry
of the problem does not greatly affect the difficulty
of performing the Monte Carlo experiment. The con
struction of a function, say neutron density, by
Monte Carlo is conceptually little more difficult for
a square or hexagonal unit cell than it is for a slab
or cylindrical unit cell. This is in striking contrast
to the analytic description of the same function;
there the solution for simple geometries does not, in
general, give any indication of how to proceed to more
complex geometries. Thus there is strong motivation
to generate the plane geometry Green's function G(v;x x')
by Monte Carlo. The transition to cylindrical or even
more complex unit cell geometries would not be difficult;
the general procedure, that is, the Monte Carlo sampling
operations would be essentially the same.
The Monte Carlo generation of G(v;xix'), that
is, the first flight Green's function for a homogeneous
finite medium in an infinite plane geometry with
perfectly reflecting boundaries has been discussed
extensively by Aswad (25) and will be outlined below.
In plane geometry, a cluster of neutrons is
considered to be emitted with unit statistical weight
at the source plane x'. Its direction cosine, P, is
chosen from a random distribution, the distribution
being such as to represent an isotropic source, that is,
a source which has equal probability for emission of
neutrons in an element of solid angle dn about any
direction n; thus in plane geometry the direction
cosine P occurs with equal probability for values
between 1 and +1. The random selection of P is
done by Monte Carlo sampling.
The projection of the path of the neutron
cluster on the x axis is divided into closed intervals
[Xml, xmJ [xm xm+1 ... and the statistical
weight of the neutron cluster is computed at points
Xml, xm, xm+l, *. The weight of the neutron cluster
at some point Xm, denoted by W(xmjx'), is subject only
to exponential attenuation along the path of travel, the
mean free path between removal collisions being equal
to the reciprocal of the total collision probability
Zt(k',v). Thus
W(xmIx') e t ( ')/p (129)
If xm+l > xm x', and if W(xmlx') and W(xim+xl')
represent respectively the statistical weight of the
neutron cluster at x, and xm+1, then the quantity
[w(xmJx') W(xm1lx')] must represent the number of
removal collisions per unit time that occurred in the
closed interval [xm, xm+ .
The neutron cluster continues along its
original path until it arrives at a cell boundary
where it is considered to be reflected and sent back
across the cell with direction cosine p', where u*u,
that is, Ip' Iu." The number of collisions that
occur in the interval [xm, Xm+1 as the neutron cluster
passes through that interval is again recorded. The
process is repeated as the neutron cluster suffers
