ELEMENTARY THEORY OF THE BOILER*)

G. Soodak and E. Campbell

1. INTRODUCTION

Our aim is to consider the chain reaction in a boiler in which, as a result of fission, fast neutrons are born. Some of these neutrons, after being slowed down in the moderator, are captured by other fissioning nuclei, which leads to the appearance of new neutrons and to the development of the chain reaction. The latter occurs only when just a small fraction of the neutrons escapes beyond the boundaries of the boiler or is absorbed unproductively before fission takes place. In order for the boiler to operate at a constant neutron density, an exact balance must be achieved between the rate of neutron loss and the rate of their production.

This can be expressed by means of formulas as follows: let \(L\) be the number of neutrons leaving the boiler each second, and \(A\) the total number of neutrons absorbed in one second in the boiler. Some part of these latter, \(A_f\), produces fission events, in each of which \(\nu\) neutrons arise. One may write the balance equation:

\[ \text{leakage} + \text{absorption} = \text{production}, \]

i.e.

\[ L + A = A_f \nu . \tag{1,1} \]

We must investigate in detail what happens to the fast neutrons that arise in the process of fission, how they move in space, and how they are slowed down in matter. A fast neutron, passing through a moderator, for example graphite, will collide with carbon nuclei. The path of the neutron will consist of many short segments of different lengths.

*) H. Soodak and E. C. Campbell, Elementary Pile Theory. New York—London, 1950. Translated from English.

These segments will be oriented more or less randomly, since in a collision a neutron may either be deflected through a small angle or undergo a head-on impact with a large deflection and maximum loss of energy.

2. CROSS SECTIONS

It is customary to express the probability of some process (scattering, absorption, fission) occurring when nuclear particles pass through matter by specifying the effective area of the bombarded nucleus for this process. This quantity is called the cross section of the nucleus for this process and is usually denoted by the letter $\sigma$. If one imagines nuclei as small solid balls of radius $R$, then one may expect the effective area of the target to be of order $\pi R^2$. Since it is known from other data that nuclear radii are of order $10^{-12}\ \text{cm}$, it is not surprising that some cross sections for neutron scattering are of order $10^{-24}\ \text{cm}^2$. This quantity is called the “barn.”

The exact meaning of the concept of cross section is as follows: if a beam of neutrons with density $I_0$ (neutrons per $1\ \text{sec}$ through $1\ \text{cm}^2$) passes through a region in which there are $N$ nuclei per $1\ \text{cm}^3$, each with scattering cross section $\sigma_s$, then the number of scattering processes occurring in $1\ \text{sec}$ in $1\ \text{cm}^3$ will be $I_0 N\sigma_s$. This concept can be illustrated by the practically unrealizable case in which $I_0=1$ and $N=1$. If one neutron enters $1\ \text{cm}^3$ in $1\ \text{sec}$ and there is only one nucleus in this volume, then the probability that it will be scattered is equal to the ratio of the area of the nuclear target to the entire cross section of the “beam,” i.e. to $1\ \text{cm}^2$. Thus this probability is numerically equal to $\sigma_s$.

Let a beam with flux $1\ \text{neutron}/\text{sec}\cdot\text{cm}^2$ enter a plate made of moderator; we shall denote the probability of neutron scattering in a thickness $\Delta x$ by $\Sigma_s\Delta x$. The probability that it passes through $\Delta x$ without scattering will be $(1-\Sigma_s\Delta x)$, and the probability that it passes through $n$ such thicknesses will be $(1-\Sigma_s\Delta x)^n$. Let $x$ be the total thickness of moderator traversed by the neutron; then, as $n$ tends to infinity, $\Delta x=\dfrac{x}{n}$ tends to zero, and the probability that the neutron passes through the thickness $x$ without scattering will be equal to $e^{-\Sigma_s x}$, since

\[ \lim_{y\to 0}(1-y)^{\frac{a}{y}}=e^{-a}. \]

From the preceding arguments it is clear that $\Sigma_s$ is related to the cross section $\sigma_s$ by the relation $\Sigma_s=N\sigma_s$.

The fraction of neutrons that travel a distance $x$ without scattering will be equal to $e^{-\Sigma_s x}$. The mean free path length

for scattering \(\lambda_s\) is defined as the distance that, on average, a neutron travels before scattering, i.e.

\[ \lambda_s=\int_0^\infty x e^{-\Sigma_s x}\Sigma_s\,dx=\frac{1}{\Sigma_s}. \tag{2,1} \]

In this integral, \(e^{-\Sigma_s x}\) is the probability that the neutron will reach \(x\) without scattering, while \(\Sigma_s dx\) is the probability that it will be scattered in the following interval \(dx\); therefore the product of these quantities gives the probability of scattering between \(x\) and \(x+dx\); the integral then gives the mean free path for scattering \(\lambda_s\), which is equal to the reciprocal of \(\Sigma_s\)—the probability of scattering per \(1\ \mathrm{cm}\) of path.

The quantity \(\Sigma_s=N\sigma_s\) is called the “macroscopic scattering cross section” or the “scattering cross section in \(1\ \mathrm{cm}^3\).” In the same way one may define analogous quantities for other processes, for example for absorption, denoting them by \(\lambda_a\), \(\sigma_a\), and \(\Sigma_a\). The number of atoms in \(1\ \mathrm{cm}^3\), \(N\), is equal to the number of moles of substance in \(1\ \mathrm{cm}^3\), multiplied by the number of atoms in 1 mole. Thus,

\[ N=\frac{\text{density }(\mathrm{g}/\mathrm{cm}^3)}{\text{atomic weight }(\mathrm{g}/\mathrm{mole})}\times \text{Avogadro's number }(\text{atoms}/\mathrm{mole}). \tag{2,2} \]

In the case of graphite we obtain:

\[ N=\frac{1.6}{12}\cdot 6.02\cdot 10^{23}=8\cdot 10^{22}. \]

Since the cross section for scattering of thermal neutrons in graphite is \(4.8\) barns, the macroscopic cross section is

\[ \Sigma_s=N\sigma_s=8\cdot 10^{22}\cdot 4.8\cdot 10^{-24}=0.38\ \mathrm{cm}^{-1}. \]

The mean free path \(\lambda_s\) will be equal to \(1/0.38=2.7\ \mathrm{cm}\). The fact that graphite is usually used as a moderator is partly due to its extremely small absorption cross section, which is only one thousandth of the scattering cross section, namely \(\sigma_a=0.0048\) barn. The mean free path for absorption of thermal neutrons in graphite will therefore be about \(2700\ \mathrm{cm}\).

In the general case, both \(\sigma_s\) and \(\sigma_a\) depend on the energy of the neutron, and it must be specified explicitly when giving the cross section. The absorption cross section is a more sensitive function of neutron energy than the scattering cross section. For many substances, for example boron, \(\sigma_a\) is inversely proportional to the neutron velocity \(v\). Other substances exhibit the phenomenon of “resonance,” in which \(\sigma_a\) reaches very large values for neutron energies close to some definite value.

3. SLOWING DOWN OF NEUTRONS

Let us consider a fast neutron \((v \sim 10^9\ \text{cm/sec})\) in a large block of moderator, for example graphite. We must now investigate in detail the processes of elastic collisions, as a result of which the neutron passes from high energies to low ones. When the energy of the neutron becomes so small that it is comparable with the energy of the thermal motion of the graphite atoms, then subsequently the neutron can with equal probability either gain or lose energy. Therefore, in the thermal region it will keep its energy constant until it is, finally, captured.

A simple application of the laws of conservation of momentum and energy leads to a relation expressing the relative

	System \(L\) (laboratory)	System \(C\) (center of inertia)
Before collision	neutron \(m\) with speed \(v_0\); center-of-inertia speed \(V_c=\left(\dfrac{m}{M+m}\right)v_0\); nucleus \(M\)	neutron \(m\) with speed \(v_0 - V_c\); center of inertia; nucleus \(M\) with speed \(V_c\)
After collision	neutron \(m\) and nucleus \(M\) recoil at angle \(\theta\)	neutron velocity \(v - V_c\) at angle \(\varphi\); nucleus \(M\) with speed \(V_c\)

Fig. 1.

energy loss of the neutron in terms of the neutron deflection angle and the ratio of the mass of the nucleus to the mass of the neutron.

However, it is more convenient to consider the collision from the point of view of an observer who moves together with the center of inertia of both particles. In this system (denoted \(C\)) the total momentum (the vector sum) of the particles before the collision is equal to zero. From the law of conservation of momentum it follows that after the collision the total momentum is again equal to zero. This means that an observer moving together with the center of inertia sees these two particles fly apart in strictly opposite directions. If, moreover, the collision is elastic (kinetic energy is conserved), then the velocities of the particles remain the same as they were before the collision, since a change in velocity would mean a change in the total kinetic energy of these two particles. The final result, observed in system \(C\), consists only in a change of direc-

of velocities, but not their magnitudes. In the laboratory system (denoted by \(L\)), in which the nucleus was initially at rest, the magnitudes of the velocities change, while their directions do not change to the opposite ones. In order to determine the new velocity of the neutron in system \(L\), we must perform the inverse transformation from system \(C\) to system \(L\).

Let the neutron move to the right with velocity \(v_0\), energy \(E_0\), and mass \(m\). The nucleus has velocity equal to zero and mass \(M\). The velocity of the center of inertia will be:

\[ V_C=\left(\frac{m}{M+m}\right)v_0. \]

In system \(C\) the neutron moves to the right with velocity

\[ v_0 - V_C=\left(\frac{M}{m+M}\right)v_0, \]

and the nucleus moves to the left with velocity

\[ V_C=\left(\frac{m}{m+M}\right)v_0. \]

The total momentum will then be equal to zero, since

\[ m\frac{M}{M+m}v_0 - M\frac{m}{M+m}v_0=0. \]

After the collision the neutron flies off in system \(C\) at an angle \(\varphi\), and the nucleus at an angle \(180^\circ+\varphi\). In system \(L\) the neutron flies off at an angle \(\theta\) and has velocity \(v\), which is the vector sum of the velocity of the neutron in system \(C\) and the velocity of the center of inertia. Especially interesting are two cases:

Fig. 2.

for a grazing collision

\[ \varphi=0,\quad v=v_0 \quad \text{and} \quad E=E_0; \]

in the case of a head-on collision

\[ \varphi=180^\circ,\quad v=\left(\frac{M-m}{M+m}\right)v_0 \quad \text{and} \quad E=\left(\frac{M-m}{M+m}\right)^2 E_0. \]

It is obvious that in the latter case the neutron loses the maximum energy in the collision. For the case of carbon

\[ E=\left(\frac{12-1}{12+1}\right)^2 E_0=0.72E_0. \]

Thus, in the collision of a neutron with a carbon nucleus the neutron can lose up to \(28\%\) of its initial energy. This means that in such a collision a neutron possessing an energy of \(1\) MeV can lose up to \(0.28\) MeV, while a neutron with an energy of \(1\) eV can lose up to \(0.28\) eV. Since the maximal

if the relative loss is constant, then it is convenient to use a logarithmic energy scale for calculations. From Fig. 2 one can obtain, by means of the cosine theorem,

\[ v^2 = v_0^2\left(\frac{M}{M+m}\right)^2 + v_0^2\left(\frac{m}{M+m}\right)^2 + \]

\[ +\,2v_0^2\left(\frac{M}{M+m}\right)\left(\frac{m}{M+m}\right)\cos\varphi, \tag{3,1} \]

and the ratio of the energy of the neutron after the collision \(E\) to its initial energy \(E_0\) is equal to:

\[ \frac{E}{E_0} = \frac{\frac{1}{2}mv^2}{\frac{1}{2}mv_0^2} = \frac{v^2}{v_0^2} = \frac{M^2+m^2+2mM\cos\varphi}{(M+m)^2}. \tag{3,2} \]

Introducing the mass ratio \(A=\frac{M}{m}\) and the quantity \(r=\left(\frac{A-1}{A+1}\right)^2\), one can transform (3,2) to the form:

\[ \frac{E}{E_0}=\frac{1+r}{2}+\frac{1-r}{2}\cos\varphi. \tag{3,3} \]

The smallest value of \(E\) is obtained for \(\varphi=180^\circ\), when \(\cos\varphi=-1\), and \(E=rE_0\), whereas for \(\varphi=0\), \(\cos\varphi=1\) and \(E=E_0\). The angle of deflection in the system \(C\) is connected, according to Fig. 2, with the corresponding angle \(\theta\) in the system \(L\) by the equation

\[ \operatorname{ctg}\theta = \frac{\cos\varphi+\frac{1}{A}}{\sin\varphi}, \]

\[ \cos\theta = \frac{1+A\cos\varphi}{\sqrt{1+A^2+2A\cos\varphi}}. \tag{3,4} \]

It may be noted that if \(A\gg 1\), then \(\varphi=\theta\) and the systems \(C\) and \(L\) are almost identical.

In order to obtain the averaged properties of neutrons slowing down in a moderator, it is necessary to know exactly how the probability of scattering in the system \(C\) depends on the angle \(\varphi\). From theory and experiment it follows that for neutron energies with which we shall be concerned (energies less than \(10\ \mathrm{Mev}\)), and for small \(A\), scattering in the system \(C\) may, to a good approximation, be considered spherically symmetric.

This means that the differential cross section \(d\sigma_s\) for neutrons scattered within the solid angle \(d\Omega\) will be \((\sigma_s/4\pi)d\Omega\), where \(\sigma_s\) is constant. Since the element of solid angle contained between \(\varphi\) and \(\varphi+d\varphi\) is

\[ 2\pi\sin\varphi\,d\varphi=-2\pi d(\cos\varphi), \]

then all values of \(\cos \varphi\) have equal probability. It is not difficult to see that, since \(E/E_0\) is a linear function of \(\cos \varphi\), all values of \(E/E_0\) from 1 to \(r\) are equally probable. The probability \(P\,dE\) that, in one collision, the neutron will lose energy from the initial value \(E_0\) to a final value lying in the interval from \(E\) to \(E+dE\) is therefore equal to \(dE\) divided by the full interval of energy values, \(E_0-rE_0\), which it may have after the collision. Thus:

\[ P\,dE=\frac{dE}{E_0-rE_0}. \]

Fig. 3.

Now we are in a position to calculate the quantity called the “mean logarithmic energy loss in one collision” and usually denoted by \(\xi\). The convenience of this quantity \(\xi\) is connected with the fact that it does not depend on the energy of the neutron. By definition,

\[ \xi=\overline{\ln E_0-\ln E}=\overline{\ln \frac{E_0}{E}}, \]

i.e.

\[ \xi=\int_{rE_0}^{E_0}\ln\frac{E_0}{E}\,P\,dE = \int_{rE_0}^{E_0}\ln\frac{E_0}{E}\,\frac{dE}{E_0-rE_0}. \tag{3,5} \]

Let

\[ \frac{E}{E_0}=x. \]

Then

\[ \xi=\frac{1}{1-r}\int_1^r \ln x\,dx, \]

or

\[ \xi=1+\frac{r}{1-r}\ln r, \tag{3,6} \]

where

\[ r=\left(\frac{A-1}{A+1}\right)^2. \]

A convenient approximate expression, accurate to 1% for \(A>10\), has the form:

\[ \xi=\frac{2}{A+\frac{2}{3}}. \tag{3,7} \]

For \(A=1\) \((r=0)\) and for \(A=\infty\) \((r=1)\) this function is undefined. We can, however, define it for these two values by passage to the limit. This gives \(\xi=1\) for \(A=1\) and \(\xi=0\) for \(A=\infty\).

The first case corresponds to the use of hydrogen as a moderator. The value \(\xi=1\) for this case means that, on the average, the energy of a neutron in collisions with hydrogen nuclei is reduced by a factor \(e\) at each collision; this means that the energy after a collision is only \(37\%\) of the initial energy.

On the other hand, if \(A\) is very large, then \(\xi=0\). The neutron practically does not lose energy in an elastic collision with a heavy nucleus, for example with \(U^{238}\).

Although in the system \(C\) the scattering is spherically symmetric, this is in general not the case in the system \(L\). The deviation from spherical symmetry is measured by the mean value of \(\cos\theta\), where the average is taken over all possible collisions. Using equation (3.4), we obtain:

\[ \overline{\cos\theta} = \frac{1}{4\pi}\int \cos\theta\,d\Omega = \]

\[ = \frac{1}{2}\int_{0}^{\pi} \frac{1+A\cos\varphi}{\sqrt{1+A^{2}+2A\cos\varphi}} \sin\varphi\,d\varphi = \]

\[ = \frac{1}{2}\int_{-1}^{1} \frac{1+Ax}{\sqrt{1+A^{2}+2Ax}}\,dx, \]

or

\[ \overline{\cos\theta}=\frac{2}{3A}. \tag{3.8} \]

When \(A\) is very large, \(\overline{\cos\theta}\) is very small, and the angular distribution of deflection is isotropic. This agrees with the preceding result for \(A\gg 1\). In this case the systems \(L\) and \(C\) coincide. Neutrons, colliding with heavy nuclei, are scattered equally often forward (positive \(\overline{\cos\theta}\)) and backward (negative \(\overline{\cos\theta}\)). In the case of hydrogen \((A=1)\), \(\overline{\cos\theta}=2/3\), and scattering occurs predominantly forward in the system \(L\), although it remains spherically symmetric in the system \(C\).

If we know the value of \(\xi\) (equation 3.6) for a moderator, we can very easily calculate the mean number of collisions that a neutron must undergo in order to reduce its energy from \(2\) MeV to thermal \((1/30\ \text{eV})\). The number of collisions is equal to the total loss of the logarithm of the energy divided by the mean

loss \(\xi\) in one collision. This gives, for the average number of collisions,

\[ \frac{\ln 2\cdot 10^6-\ln \frac{1}{30}}{\xi} = \frac{\ln 6\cdot 10^7}{\xi} = \frac{18}{\xi}. \tag{3.9} \]

Some numerical examples are given in Table 1.

Table 1

Moderator	\(A\)	\(\xi\)	Number of collisions corresponding to an energy loss from \(2\) MeV to \(\frac{1}{30}\) eV, equal to \(18/\xi\)
Hydrogen	1	1	18
Deuterium	2	0.725	25
Carbon	12	0.158	114
Beryllium	9	0.209	86

4. SLOWING-DOWN DENSITY

Let us consider a block of moderator whose volume is so large that we may neglect the leakage of neutrons from it. Suppose that in each cubic centimeter, every second, \(q\) neutrons with energy \(E_0\) are born. We shall be interested in the distribution of neutrons with respect to energy. In doing so we shall assume that the spatial distribution is isotropic. Then, in the process of slowing down (provided that there is no absorption of neutrons in the moderator), in \(1\ \mathrm{cm}^3\) in \(1\ \mathrm{sec}\), \(q\) neutrons will be slowed down to any value of the energy \(E\). The quantity \(q\) is called the slowing-down density. If absorption exists, then fewer neutrons will leave the interval \(dE\) than entered this interval, and \(q\) will be a function of the energy \(E\). In the general case, \(q(E)\) is defined as the number of neutrons slowed down in \(1\ \mathrm{sec}\) in \(1\ \mathrm{cm}^3\) from energies greater than \(E\) to energies less than \(E\).

Under stationary conditions, the number of neutrons leaving the interval \(dE\) must be equal to the number of neutrons entering it. The number of neutrons leaving is equal to the number scattered, i.e. \(n(E)\,dE\cdot v\Sigma_s\), where \(n(E)\,dE\) is the number of neutrons in \(1\ \mathrm{cm}^3\) having energy in the interval between \(E\) and \(E+dE\), and \(v\Sigma_s\) is the probability of scattering referred to \(1\ \mathrm{sec}\). This number is equal to the number of neutrons entering this interval as a result of the scattering of neutrons

from higher energies. Consider those of them which enter from the interval \(dE'\), lying somewhere between \(E\) and \(E/r\). (Neutrons with energy greater than \(E/r\) cannot pass into \(dE\).)

The number of neutrons scattered into the interval \(dE\) from the region of higher energies is equal to the integral over \(dE'\) of the number of scattering events in \(dE'\) (the prime refers to the value of the energy \(E'\)), equal to

\[ n(E')v'\Sigma_s' dE', \]

multiplied by the probability that the energy loss will be just such as to transfer the neutron from \(dE'\) into \(dE\). This probability is equal to

\[ \frac{dE}{E' - E'r}. \]

Fig. 4.

The balance equation can be written in the form:

\[ n(E)v\Sigma_s\, dE = \int_E^{E/r} n(E')\, dE'\, v'\Sigma_s'\,\frac{dE}{E' - E'r}. \tag{4,1} \]

This equation is satisfied if one sets \(n(E)v\Sigma_s = c/E\), where \(c\) is a constant, which can be verified by direct substitution:

\[ \int_E^{E/r} \frac{c}{E'}\,\frac{dE}{E' - E'r} = \frac{c}{E} = n(E)v\Sigma_s . \tag{4,2} \]

The quantity \(c\) can be expressed in terms of the slowing-down density \(q(E)\).

Let us proceed to the calculation of \(q(E)\). The number of neutrons (per \(1\) sec in \(1\ \mathrm{cm}^3\)) slowed down below \(E\) and scattered from the interval \(dE'\) is equal to the number of scattering events in \(dE'\),

\[ n(E')v'\Sigma_s' dE', \]

multiplied by the fraction of neutrons from \(dE'\) which lose an energy greater than \(E' - E\) and enter the shaded region in Fig. 5. This fraction is equal to

\[ \frac{E - E'r}{E' - E'r}. \]

Integrating over \(E'\) from \(E\) to \(E/r\), we obtain:

\[ q(E) = \int_E^{E/r} n'v'\Sigma_s' dE'\,\frac{E - E'r}{E' - E'r}, \tag{4,3} \]

that by the substitution \(n(E')v'\Sigma'_s=c/E'\) reduces to

\[ q=c\left[1+\frac{r}{1-r}\ln r\right]. \tag{4,4} \]

The expression in brackets is \(\xi\)—the mean logarithmic energy loss per collision. As we expected, \(q\) does not depend on the energy, and the constant \(c=q/\xi\). The equation for the “flux” has the form:

\[ n(E)v=\frac{q}{\xi\Sigma_s}\frac{1}{E}. \tag{4,5} \]

Fig. 5.

The same result can be obtained by a simpler, though less rigorous, method, using the scheme proposed by Fermi. Let us imagine that the neutron is slowed down not in jumps, but continuously. Let the slowing-down density be equal to \(q\). What is the value of \(n(E)\,dE\), if each second \(q\) neutrons enter \(dE\) and leave \(dE\)? The answer depends on the average time a neutron requires to traverse this interval. If the neutron is delayed in the interval under consideration for \(t\) seconds, then the number of neutrons in the interval will be equal to \(qt\). The quantity \(t\) can be expressed as the product of the mean number of collisions undergone by the neutron in \(dE\) \(\left(\text{equal to } \frac{1}{\xi}\frac{dE}{E}\right)\), and the interval between two successive collisions (equal to \(\lambda_s/v\)). We then obtain

\[ n(E)\,dE=qt=q\frac{1}{\xi}\frac{dE}{E}\frac{\lambda_s}{v} \tag{4,6} \]

or, since \(\lambda_s=1/\Sigma_s\),

\[ n(E)v=\frac{q}{\xi\Sigma_s}\frac{1}{E}, \tag{4,7} \]

which coincides with (4,5). The quantity \(\xi\Sigma_s\) is called the slowing-down power of the moderator. Since \(\Sigma_s=N\sigma_s\) is the probability that a neutron will undergo a collision on a path of \(1\ \mathrm{cm}\), and \(\xi\) is the mean logarithmic energy loss per collision, the slowing-down power is equal to the logarithmic loss of energy (the decrease in \(\ln E\)) per \(1\ \mathrm{cm}\) of path. A good moderator, therefore, must have a relatively high value of \(\xi\Sigma_s\).

Let us note that the neutron “flux” \(n(E)v\) in the energy interval \(dE\) is inversely proportional to the energy.

5. SLOWING-DOWN IN THE PRESENCE OF ABSORPTION

In the preceding calculations it was assumed that absorption is absent, so that the slowing-down density \(q\) was constant. We shall now consider the case in which the moderator absorbs neutrons with an effective cross section \(\Sigma_a\), which, in general, is a function of energy. The slowing-down density \(q(E)\) decreases as \(E\) decreases. The decrease in \(q\) \(\left(dq=\dfrac{dq}{dE}\cdot dE\right)\) is equal to the number of neutrons in \(1\ \mathrm{cm}^3\) absorbed in \(1\ \mathrm{sec}\) in the interval \(dE\), i.e. \(n(E)v\Sigma_a dE\). Therefore we have

\[ n(E)v\Sigma_a=-\frac{dq}{dt}. \tag{5,1} \]

The balance equation (4,1) may be rewritten in the form

\[ n(E)v\Sigma_s\,dE=\frac{q}{\xi}\,\frac{dE}{E}. \]

This expression must now be modified so as to take account of absorption in the moderator. It expresses the balance between the number of neutrons (in \(1\ \mathrm{sec}\) in \(1\ \mathrm{cm}^3\)) that have left \(dE\) as a result of scattering, and the number of neutrons that have entered this interval as a result of scattering from higher energies. In our case the loss of neutrons in \(dE\) occurs both as a result of scattering and as a result of absorption. We may therefore write:

\[ n(E)v(\Sigma_s+\Sigma_a)\,dE=\frac{q(E)}{\xi}\,\frac{dE}{E}. \tag{5,2} \]

After multiplying both sides of (5,2) by \(\dfrac{\Sigma_a}{\Sigma_a+\Sigma_s}\) and using (5,1), we obtain a simple differential equation for \(q(E)\)

\[ \frac{dq}{dE}=\frac{\Sigma_a}{\Sigma_s+\Sigma_a}\,\frac{q(E)}{\xi E}, \tag{5,3} \]

which can be integrated from \(E\) to \(E_0\) to give:

\[ \ln \frac{q(E_0)}{q(E)} = -\frac{1}{\xi}\int_E^{E_0} \frac{\Sigma_a}{\Sigma_s+\Sigma_a}\,\frac{dE}{E}. \]

Solving for \(q(E)\), we obtain:

\[ \frac{q(E)}{q(E_0)} = e^{-\frac{1}{\xi}\int_E^{E_0} \frac{\Sigma_a}{\Sigma_s+\Sigma_a}\,\frac{dE}{E}}. \tag{5,4} \]

In the case of taking \(\Sigma_a\) to zero, equation (5.4) reduces to equation (4.7), which was obtained without taking account of the absence of absorption. If \(\Sigma_a\) and \(\Sigma_s\) are known as functions of the energy \(E\), then \(q(E)\) can be calculated. Since \(\Sigma_s\) is usually a slowly varying function of the energy, one may determine a suitable mean value \(\overline{\Sigma_s}\) and take it out from under the integral, obtaining:

\[ \frac{q(E)}{q(E_0)} = e^{ -\frac{1}{\xi \overline{\Sigma_s}} \int_E^{E_0} \frac{\Sigma_a}{1+\Sigma_a/\Sigma_s} \frac{dE}{E} } . \tag{5.5} \]

The meaning of this formula is easy to understand if one also determines the mean value \(\overline{\Sigma_a}\). However, since \(\Sigma_a\) is, generally speaking, a rapidly varying function of the energy (especially in the case of resonance absorption), this operation is not performed in actual calculations. The integral of \(dE/E\) gives \(\ln E_0/E\), and equation (5.5) reduces to

\[ \frac{q(E)}{q(E_0)} = e^{ -\frac{\overline{\Sigma_a}}{\xi \overline{\Sigma_s}} \ln \frac{E_0}{E} \cdot \frac{1}{1+\Sigma_a/\Sigma_s} } . \tag{5.6} \]

As was noted earlier, \(\xi\Sigma_s\)—the slowing-down power—is equal to the average decrease of \(\ln E\) per \(1\ \mathrm{cm}\) of path. The factor

\[ X=\frac{1}{\xi\Sigma_s}\ln\frac{E_0}{E}, \]

therefore, is the mean distance traversed by a neutron during slowing down from \(E_0\) to \(E\). As before, the fraction of neutrons which has slowed down from an energy greater than \(E\) to an energy less than \(E\) is equal to:

\[ \frac{q(E)}{q(E_0)}=e^{-\overline{\Sigma_a}X}, \tag{5.7} \]

if \(\Sigma_a/\Sigma_s \ll 1\). The factor \(\dfrac{1}{1+\Sigma_a/\Sigma_s}\) is a measure of the self-shielding of the absorber. For a strongly diluted absorber the number of moderator atoms is so much greater than the number of absorber atoms that this factor is very close to unity.

Up to now we have considered only a homogeneous mixture of moderator and absorber atoms. Additional possibilities are obtained if the absorber, instead of being distributed uniformly, is collected into “blocks.” The advantage of such an arrangement in a graphite–uranium reactor is connected with the fact that resonance neutrons are very rapidly absorbed by \(U^{238}\). By collecting the uranium into blocks, one can increase the fraction of neutrons which avoid such “resonance leakage,” and they may be absorbed by \(U^{235}\), causing fission. If a resonance neutron enters a uranium block, it is absorbed by a thin outer layer of \(U\), which acts as

filter that screens out such neutrons. Neutrons of higher energy are slowed down by uranium only slightly; it is therefore clear that the resonance absorption per uranium atom is considerably smaller in a system with blocks than in the case of a homogeneous mixture.

According to the types of reactors, they are divided into reactors with slow, fast, or intermediate neutrons, depending on whether absorption takes place when the neutrons are in the slowed, fast, or intermediate energy region.

In each individual case it is important to know what fraction of the neutrons is found in the various energy intervals. This depends on how the cross section varies with energy.

From equations (5, 7) it can be seen that an important characteristic determining in which energy region neutrons are absorbed is the ratio of the average macroscopic absorption cross section to the moderating power of the moderator. If this ratio is small (high moderating power), then the neutron need not traverse a long path in order for its energy to reach a small value, and along this short path it has little chance of being absorbed. In this case the majority of neutrons can reach thermal energies before they are absorbed. In the other limiting case, if there is no moderator at all, the neutrons will be absorbed (if they are absorbed at all) at high values of energy.

6. INTRODUCTION TO DIFFUSION THEORY

We shall now return to the problem of the distribution of neutrons in space, knowledge of which is necessary if we wish to calculate the leakage from a reactor or its critical size.

Suppose that we are dealing with a substance characterized by the macroscopic absorption cross section \(\Sigma_a\), the macroscopic scattering cross section \(\Sigma_s\), and a constant multiplication factor \(k\), where \(k\) is the number of neutrons produced per absorbed neutron. In this case the number of neutrons absorbed per second in one cubic centimeter of the substance will be \(nv\Sigma_a\), whereas the number of neutrons produced per second will be \(nv\Sigma_a k\). The quantity \(nv\) is usually called the “neutron flux.” Its meaning is easy to understand if one notes that \(nv\) is equal to the total length of path traversed per second by all neutrons in \(1\ \mathrm{cm}^3\). Since \(\Sigma_a\) is the probability of absorption of a neutron on a path of \(1\ \mathrm{cm}\), the product \(nv\Sigma_a\) is the number of neutrons absorbed per second in \(1\ \mathrm{cm}^3\).

If we temporarily forget that these neutrons were born fast, then the balance equation (1,1), written for \(1\ \mathrm{cm}^3\), will have the form:

\[ nv\Sigma_a k = nv\Sigma_a + L \]

or

\[ -L+(k-1)\Sigma_a nv=0. \tag{6,1} \]

The problem is to obtain an expression for the leakage \(L\) per unit volume as a function of the flux \(nv\) and its spatial variation.

We can calculate the number of neutrons that pass per second from above through an elementary area \(dS\), whose normal is directed along the \(z\)-axis, as shown in Fig. 6. The number of such neutrons emerging from the volume element \(dV\) is equal to the number of collisions per second in \(dV\),

\[ nv\Sigma_s\,dV, \]

multiplied by the probability that the scattered particle moves in such a direction as to cross \(dS\),

\[ \frac{|\cos\theta|\,dS}{4\pi r^2}, \]

and multiplied by the probability that the particle traverses the distance \(r\) without being scattered,

\[ e^{-r/\lambda_s}. \]

Fig. 6.

In order to find the total current \(I_+ dS\), we must also integrate over the entire volume lying above \(dS\). Then we obtain:

\[ I_+ dS=dS\int nv\Sigma_s e^{-r/\lambda_s}\frac{|\cos\theta|}{4\pi r^2}\,dV. \]

In polar coordinates \(dV=r^2\sin\theta\,d\theta\,dr\,d\varphi\), so that:

\[ I_+=\frac{1}{4\pi}\int_0^{2\pi}\int_0^\infty\int_0^{\pi/2} nv\Sigma_s e^{-r/\lambda_s}\frac{|\cos\theta|}{r^2}r^2\sin\theta\,d\theta\,dr\,d\varphi. \tag{6,2} \]

Now we must make certain assumptions regarding the dependence of \(nv\) on the coordinates \(x, y\), and \(z\), substitute the corresponding expression into integral (6,2), and carry out the integration.

Because of the factor \(e^{-r/\lambda_s}\) in the integrand, the main contribution to the integral comes from values of \(nv\) in a region several mean free paths away from the origin. Therefore one may expand \(nv(x,y,z)\) in a Maclaurin series up to

terms of second order:

\[ nv(x,y,z)=nv_0+x\left(\frac{\partial nv}{\partial x}\right)_0+y\left(\frac{\partial nv}{\partial y}\right)_0+z\left(\frac{\partial nv}{\partial z}\right)_0+ \]
\[ +\frac{1}{2!}\left\{x^2\left(\frac{\partial^2 nv}{\partial x^2}\right)_0+y^2\left(\frac{\partial^2 nv}{\partial y^2}\right)_0+z^2\left(\frac{\partial^2 nv}{\partial z^2}\right)_0+2xy\left(\frac{\partial^2 nv}{\partial x\partial y}\right)_0+\right. \]
\[ \left.+2xz\left(\frac{\partial^2 nv}{\partial x\partial z}\right)_0+2yz\left(\frac{\partial^2 nv}{\partial y\partial z}\right)_0\right\}. \tag{6,3} \]

The index 0 means that the derivatives must be taken at the origin of coordinates. Express \(x, y\) and \(z\) in polar coordinates \(r, \theta\) and \(\varphi\):

\[ x=r\sin\theta\cos\varphi, \]

\[ y=r\sin\theta\sin\varphi, \]

\[ z=r\cos\theta. \]

The integration can be carried out elementarily. We obtain:

\[ J_+=\frac{nv_0}{4}+\frac{\lambda_s}{6}\left(\frac{\partial nv}{\partial z}\right)_0+ \]
\[ +\frac{\lambda_s^2}{16}\left\{\left(\frac{\partial^2nv}{\partial x^2}\right)_0+\left(\frac{\partial^2nv}{\partial y^2}\right)_0+2\left(\frac{\partial^2nv}{\partial z^2}\right)_0\right\}. \tag{6,4} \]

Analogously, the current density \(J_-\), due to particles passing through \(dS\) from below (negative \(z\)), can be calculated. In this case the integrand is the same, and the only difference is that the integration with respect to \(\varphi\) is carried out over the limits from \(\frac{\pi}{2}\) to \(\pi\).

Carrying out the integration, we obtain:

\[ J_-=\frac{nv_0}{4}-\frac{\lambda_s}{6}\left(\frac{\partial nv}{\partial z}\right)_0+ \]
\[ +\frac{\lambda_s^2}{16}\left\{\left(\frac{\partial^2nv}{\partial x^2}\right)_0+\left(\frac{\partial^2nv}{\partial y^2}\right)_0+2\left(\frac{\partial^2nv}{\partial z^2}\right)_0\right\}. \tag{6,5} \]

The terms with \(x, y, xy, xz\) and \(yz\) give nothing, because the integral with respect to \(\varphi\) of these terms vanishes.

The total current density \(J\) in the direction \(+z\) will have the form:

\[ J=J_- - J_+=-\frac{\lambda_s}{3}\left(\frac{\partial nv}{\partial z}\right)_0. \tag{6,6} \]

This expression is valid up to terms of second order.

The condition for applicability of formula (6,6) is that terms of third and higher orders in expression (6,4) would make a contribution to \(J\) that would be small in comparison with the contribution from the first-order terms.

If the surface element \(dS\) is oriented so that the normal to it makes angles \(\alpha, \beta, \gamma\) with the axes \(x, y, z\), and is not directed along the \(z\)-axis, then the expression for the total current through \(dS\) can be rewritten in the form:

\[ JdS=-\frac{\lambda_s}{3}\,dS\left[\left(\frac{\partial n v}{\partial x}\right)_0\cos\alpha+ \left(\frac{\partial n v}{\partial y}\right)_0\cos\beta+ \left(\frac{\partial n v}{\partial z}\right)_0\cos\gamma\right]. \tag{6,7} \]

This expression can be written more briefly by using the differential vector operator \(\operatorname{grad}\) (sometimes denoted by \(\nabla\)). The vector \(\operatorname{grad}(nv)\) has components \(\partial nv/\partial x\), \(\partial nv/\partial y\), and \(\partial nv/\partial z\) in the directions \(x, y\), and \(z\), respectively. If \(d\mathbf{S}\) is a vector with magnitude equal to the area of the surface element and directed normal to it, then formula (6,7) may be written as

\[ \mathbf{J}\cdot d\mathbf{S}=-\frac{\lambda_s}{3}\,d\mathbf{S}\cdot \operatorname{grad}(nv), \tag{6,7a} \]

or

\[ \mathbf{J}=-\frac{\lambda_s}{3}\operatorname{grad}(nv). \tag{6,7b} \]

It should be noted that the expressions for \(J_+\) or \(J_-\) contain terms of the second order. The approximation in which only the first two terms are retained in expressions (6,4) is valid only when the change of \(\operatorname{grad}(nv)\) upon displacement in the medium over a segment \(\lambda_s\) is small in comparison with \(\operatorname{grad}(nv)\) itself.

This condition is usually not satisfied in a region with dimensions of the order of the mean free path near the boundary between two heterogeneous media or near a heavy absorber.

In deriving formula (6,6) we assumed that there is no correlation between the direction of motion of the neutron before and after a collision. In fact, we implicitly assumed that the scattering is spherically symmetric \((\overline{\cos\theta}=0)\) in the laboratory system. As was indicated above (in Section 3), this is true only for collisions with heavy nuclei. For nuclei of mass \(A\),

\[ \overline{\cos\theta}=\frac{2}{3A}, \]

and one may introduce a correction allowing for such predominance of forward scattering. For this purpose, instead of the mean scattering length \(\lambda_s\), one must use a new quantity, called the mean transport length \(\lambda_t\), defined by the formula

\[ \lambda_t=\frac{\lambda_s}{1-\overline{\cos\theta}}. \tag{6,8} \]

If the scattering is predominantly directed forward, then \(\overline{\cos\theta}\) is positive and the mean transport length is greater than the mean scattering length.

scattering length. This is easy to understand, since it means that on the average a neutron will travel, for a given number of collisions, greater distances if there is such a correlation between the directions of motion before and after a collision.

Similarly, one may define the transport cross section \(\sigma_t\) and the macroscopic transport cross section \(\Sigma_t\) by the formulas:

\[ \sigma_t=\sigma_s(1-\overline{\cos\theta}) \tag{6,8a} \]

and

\[ \Sigma_t=\Sigma_s(1-\overline{\cos\theta}). \tag{6,8b} \]

In the case of graphite \(\overline{\cos\theta}=\dfrac{2}{3A}=\dfrac{2}{3\cdot 12}=0.056\), and therefore the transport length in graphite is \(\lambda_t=2.70/(1-0.056)=2.86\ \text{cm}\).

Fig. 7.

We may also write an expression for the total current through an area \(dS\),

\[ J\,dS=-\frac{\lambda_t}{3}(nv)'dS, \tag{6,9} \]

where the derivative \((nv)'\) is taken along the direction perpendicular to \(dS\).

Suppose now that the neutron flux \(nv(x,y,z)\) is known. We want to calculate the number of neutrons leaving the volume element \(dV=dx\,dy\,dz\), located at the point \((x,y,z)\). Let us first consider the neutrons leaving through two faces \(dx\,dy\), perpendicular to the \(z\)-direction. We have:

\[ L_z\,dV=(J_{z+dz}-J_z)\,dx\,dy= \]

\[ =-\frac{\lambda_t}{3}\left\{\left(\frac{\partial nv}{\partial z}\right)_{z+dz} -\left(\frac{\partial nv}{\partial z}\right)_z\right\}dx\,dy= \]

\[ =-\frac{\lambda_t}{3}\frac{\partial^2 nv}{\partial z^2}\,dx\,dy\,dz. \]

In an analogous way one may obtain, for the number of neutrons leaving through faces perpendicular to the \(x\) and \(y\) directions,

\[ L_x dV=-\frac{\lambda_t}{3}\frac{\partial^2 nv}{\partial x^2}\,dV,\qquad L_y dV=-\frac{\lambda_t}{3}\frac{\partial^2 nv}{\partial y^2}\,dV, \]

or, for the total number of neutrons leaving a unit volume:

\[ L=L_x+L_y+L_z =-\frac{\lambda_t}{3}\left[ \frac{\partial^2 nv}{\partial x^2} +\frac{\partial^2 nv}{\partial y^2} +\frac{\partial^2 nv}{\partial z^2} \right]. \tag{6,10} \]

Usually this expression is written in the following abbreviated form:

\[ L=-\frac{\lambda_t}{3}\Delta nv, \tag{6,11} \]

where the differential operator \(\Delta \equiv \dfrac{\partial^2}{\partial x^2}+\dfrac{\partial^2}{\partial y^2}+\dfrac{\partial^2}{\partial z^2}\) is the Laplace operator in Cartesian coordinates.

From the sign of equation (6,11) one may conclude that there exists a total leakage of neutrons from an element of volume if \(\Delta nv\) is negative, since in this case \(L\) is positive (Fig. 8). Neutrons diffuse from a place with a high neutron density in the direction toward a place where the neutron density is lower. Equation

Fig. 8.

(6,11) can be obtained by taking the divergence of the vector equation (6,76), which gives, as before,

\[ L\,\operatorname{div} J = -\frac{\lambda_t}{3}\operatorname{div}\operatorname{grad} nv = -\frac{\lambda_t}{3}\Delta nv. \]

7. SOLUTION OF THE DIFFUSION EQUATION

Our balance equation (1,1) now takes the form:

\[ -\frac{\lambda_t}{3}\Delta(nv)+nv\Sigma_a=Q \]

\[ (\text{leakage}+\text{absorption}=\text{production}), \]

where \(Q\) is the number of neutrons produced in \(1\) sec in \(1\ \mathrm{cm}^3\), which may also be a function of the coordinates. Thus, we obtain:

\[ \frac{\lambda_t}{3}\Delta nv-nv\Sigma_a+Q=0. \tag{7,1} \]

In order to avoid complications, we shall assume that all neutrons have thermal energies and that all sources emit only thermal neutrons. In reality, thermal neutrons arise only as the result of the slowing down of fast neutrons. We shall first consider the case \((Q=0)\), when inside the consider-

able volume no neutrons arise. We may have, for example, a source of neutrons located outside a graphite block, while we are interested in the distribution of the neutron flux inside the graphite. Mathematically this problem is analogous to the problem of the distribution of temperature inside a body under prescribed boundary conditions. In the case of neutron diffusion, the boundary conditions determine which solutions or combinations of solutions of equation (7.1) correspond to the problem. For \(Q=0\), equation (7.1) may be written in the form

\[ \Delta nv - K^2 nv = 0, \tag{7.1a} \]

where we have put

\[ K^2=\frac{3\Sigma_a}{\lambda_t}. \]

For convenience, in Table II we shall give the solutions of the equation

\[ \Delta u-\alpha^2 u=0 \tag{7.16} \]

in coordinate systems corresponding to various forms of surface. The solution for \(\alpha^2>0\) must be referred to such fluxes in a medium as cannot sustain a chain reaction; the solutions for \(\alpha^2<0\), as is easily seen, correspond to media in which a chain reaction can develop.

Table II

Shape of boundary	Variables	Solutions for \(\alpha^2>0\)	Solutions for \(\alpha^2<0\) \((\beta=i\alpha)\)
Plane	\(x\) \(\displaystyle \Delta=\frac{\partial^2}{\partial x^2}\)	\(\displaystyle e^{\pm \alpha x}\) or \(\displaystyle \begin{Bmatrix}\operatorname{sh}\\ \operatorname{ch}\end{Bmatrix}\alpha x\)	\(\displaystyle e^{\pm i\beta x}\) or \(\displaystyle \begin{Bmatrix}\sin\\ \cos\end{Bmatrix}\beta x\)
Cylinder	\(r,z\) \(\displaystyle \Delta=\frac{\partial^2}{\partial r^2}+\frac{1}{r}\frac{\partial}{\partial r}+\frac{\partial^2}{\partial z^2}\)	\(\displaystyle e^{\pm \alpha_z z}\times \begin{Bmatrix} I_0(\alpha_r r)\\ K_0(\alpha_r r)\end{Bmatrix}\) where \(\displaystyle \alpha_r^2+\alpha_z^2=\alpha^2\)	\(\displaystyle e^{\pm i\beta_z z}\times \begin{Bmatrix} J_0(\beta_r r)\\ Y_0(\beta_r r)\end{Bmatrix}\) where \(\displaystyle \beta_r^2+\beta_z^2=\beta^2\)
Sphere	\(r\) \(\displaystyle \Delta=\frac{\partial^2}{\partial r^2}+\frac{2}{r}\frac{\partial}{\partial r}\)	\(\displaystyle \frac{1}{r}e^{\pm \alpha r}\) or \(\displaystyle \frac{1}{r}\times \begin{Bmatrix}\operatorname{sh}\alpha r\\ \operatorname{ch}\alpha r\end{Bmatrix}\)	\(\displaystyle \frac{1}{r}e^{\pm i\beta r}\) or \(\displaystyle \frac{1}{r}\times \begin{Bmatrix}\sin \beta r\\ \cos \beta r\end{Bmatrix}\)

The functions \(I_0\), \(K_0\), \(J_0\), and \(Y_0\) are Bessel functions of zero order, whose definition and tables are given in Watson’s book A Treatise on the Theory of Bessel Functions (Moscow, 1949).* The general solution

* See also Karman and Biot, Mathematical Methods in Engineering, vol. 2, where a summary of the properties of these functions is given, or Smirnov, A Course of Higher Mathematics, vol. 3, part 2. (Translator’s note.)

in any case is obtained by means of a linear combination of the solutions given in the table. For example, for the spherical case the general solution has the form:

\[ nv=A\,\frac{e^{-\alpha r}}{r}+B\,\frac{e^{-\alpha r}}{r}, \]

where \(A\) and \(B\) are arbitrary constants.

Case 1. For simplicity we shall first consider the case of a half-space filled with matter, on whose only boundary \((x=0)\) sources are uniformly distributed, giving \(Q_0\) neutrons per \(1\ \mathrm{sec}\) per \(1\ \mathrm{cm}^2\). Owing to symmetry, \(nv\) does not depend on \(y\) and \(z\), and \(\Delta(nv)\) reduces to \(\dfrac{d^2}{dx^2}(nv)\). The equation is satisfied if

\[ \frac{d^2(nv)}{dx^2}-K^2nv=0 \qquad (x>0), \tag{7,2} \]

where we have put \(Q=0\) and again introduced \(K^2=3\Sigma_a/\lambda_t\). The general solution of this equation will be

\[ nv=ae^{-Kx}+be^{-Kx}, \tag{7,3} \]

where the arbitrary constants \(a\) and \(b\) are determined by the boundary conditions. The first condition is that \(nv\) be finite everywhere. This gives \(b=0\), since otherwise \(nv\to\infty\) when \(x\to\infty\). The other boundary condition is that at \(x=0\) the flux \(J\) is equal to \(\dfrac{1}{2}Q_0\), since only one half of the neutrons produced go to the right. This determines the constant \(a\). We obtain, with the help of (6,6):

\[ J(0)=\frac{1}{2}Q_0=-\frac{\lambda_t}{3}\left[\frac{d}{dx}(nv)\right]_{x=0} =\frac{\lambda_t}{3}Ka, \]

whence

\[ a=\frac{3}{2\lambda_t K}\,Q_0. \]

The solution of the problem has the form:

\[ nv=\frac{3}{2\lambda_t K}\,Q_0e^{-Kx}. \tag{7,4} \]

Fig. 9.

It is easy to see that the neutron flux in the medium decreases exponentia-

ally with distance from the plane source. The flux, as they say, attenuates; in so doing it decreases by a factor of \(e\) over the distance \(L=1/K\).

The length \(L\) is called the “diffusion length.” Its definition is

\[ L^2=\frac{1}{K^2}=\frac{\lambda_t}{3\Sigma_a}=\frac{1}{3\Sigma_a\Sigma_t}=\frac{\lambda_t\lambda_a}{3}. \tag{7,5} \]

Up to a factor \(\sqrt{3}\), the diffusion length is the geometric mean of the mean transport length and the mean absorption length. For graphite,

\[ L=\sqrt{\frac{2.86\cdot2700}{3}}\sim 50\ \text{cm};\qquad K=\frac{1}{50}=0.02\ \text{cm}^{-1}. \]

Case 2. In an analogous way we can solve the problem of the distribution of neutrons in an infinite medium with one point source of intensity \(Q_0\) neutrons per second. In this case \(nv\) is a function of the distance \(r\) from the source. In spherical coordinates the Laplace operator is equal to

\[ \Delta=\frac{\partial^2}{\partial r^2}+\frac{2}{r}\frac{\partial}{\partial r}; \]

the remaining terms depend only on derivatives with respect to the azimuth and polar angle. If we assume that the function \(nv\) is spherically symmetric, then these terms contribute nothing. We must then solve the equation

Fig. 10.

\[ \frac{d^2}{dr^2}(nv)+\frac{2}{r}\frac{d}{dr}(nv)-K^2nv=0\qquad (r>0) \tag{7,6} \]

with the boundary conditions: (A) the total neutron current through a small sphere of radius \(a\) surrounding the source is, in the limit \(Q_0\), when the radius of the sphere tends to zero; (B) the function \(nv\) is finite everywhere.

Two solutions of this equation are:

\[ \frac{C}{r}e^{-Kr}\quad \text{and}\quad \frac{D}{r}e^{Kr}. \]

The second function does not satisfy condition (B) as \(r\to\infty\). We might also have tried to discard the first solution, on the grounds that it tends to infinity at \(r=0\). However, such a conclusion is incorrect, for the equation without sources (7,6) is not applicable at the point \(r=0\), where the source is located.

Boundary condition (A) can be written in the form:

\[ Q_0=\lim_{a\to 0}4\pi a^2J(a)=\lim_{a\to 0}4\pi a^2\left[-\frac{\lambda_t}{3}\frac{d}{dr}(nv)\right]_{r=a} \]

or, substituting

\[ nv=C\frac{e^{-Kr}}{r}, \]

we obtain:

$$ Q_0=\lim_{a\to 0}\left[-{\lambda_t\over 3}C\left(-{K e^{-Ka}\over a}-{e^{-Ka}\over a^2}\right)4\pi a^2\right]= $$

$$ =\lim_{a\to 0}\left[{4\pi C\lambda_t\over 3}(aK+1)e^{-Ka}\right] ={4\pi C\lambda_t\over 3}. $$

Thus the complete solution has the form

$$ nv={3Q_0\over 4\pi\lambda_t}\,{e^{-Kr}\over r}. \tag{7,7} $$

The physical meaning of the diffusion length \(L=1/K\) becomes clear if we use the solution (7,7) for a point source (of thermal neutrons) in an infinite medium. One can determine the mean square \(\overline{r^2}\) of the distance from the source traversed by a neutron (up to absorption). We obtain:

$$ \overline{r^2}= {\int r^2 nv\,dV\over \int nv\,dV} = {\displaystyle \int_0^\infty r^2 {e^{-{r\over L}}\over r}\,4\pi r^2dr \over \displaystyle \int_0^\infty {e^{-{r\over L}}\over r}\,4\pi r^2dr} = {6L^4\over L^2}=6L^2. $$

The integral is evaluated with the aid of the formula

$$ \int_0^\infty x^n e^{-{x\over L}}\,dx=n!L^{n+1} \qquad (n\text{ is any positive integer}). $$

Thus, the square of the diffusion length \(L\) is equal to \(1/6\) of the mean square of the distance traversed by a neutron from the point of the source, at which the neutron becomes thermal, to the point at which it is absorbed.

Case 3. Let us now consider the distribution of neutrons in a plane layer of material, infinite in the directions of the \(y\) and \(z\) axes, but having a finite thickness \(t\) in the direction of \(x\). The planes \(x=0\) and \(x=t\) are the boundaries of the layer. The equation that must be solved is the same equation (7,2), but with different boundary conditions. Let the total current density of neutrons in the direction of the positive \(x\)-axis at \(x=0\) be given and equal to \(I_0\). The boundary condition at \(x=0\) will then be \(J(0)=I_0\). At the boundary \(x=t\), which separates the substance from empty space, an approximate boundary condition will be the vanishing of the neutron flux at this boundary. A more exact condition can be obtained in the following way. All neutrons that have diffused beyond the boundary \(x=t\),

cannot return to the medium as a result of scattering. Empty space thus plays the role of an ideal absorber. The boundary condition (C) will be the condition that the current density \(J_+\) in the \(-x\) direction becomes zero at the boundary. Then

\[ J_+(t)=\left[\frac{nv}{4}+\frac{\lambda_t}{6}\frac{d}{dx}(nv)\right]_{x=t}=0, \tag{7,8} \]

whence we obtain:

\[ \left[\frac{nv}{\dfrac{d(nv)}{dx}}\right]_{x=t}=-\frac{2}{3}\lambda_t. \]

If we assume that inside the medium \(nv\) can be expressed, near \(x=t\) (for \(x<t\)), by a linear function which vanishes at \(x=t+d\) (at a distance \(d\) beyond the boundary of the medium), then

Fig. 11.

\[ nv(x)=nv(t)\left[1-\frac{x-t}{d}\right], \tag{7,9} \]

whence:

\[ \left[\frac{nv}{\dfrac{d(nv)}{dx}}\right]_{x=t}=-d \]

or

\[ d=\frac{2}{3}\lambda_t. \tag{7,10} \]

Thus, boundary condition (C) is equivalent to the following: the neutron flux in a medium with a plane boundary (with vacuum) changes in such a way that the linearly extrapolated distribution becomes zero at a distance \(d=\dfrac{2}{3}\lambda_t\) beyond the boundary. This distance is called the extrapolated boundary. A more rigorous treatment of the problem leads to the value

\[ d=0.71\lambda_t. \tag{7,11} \]

It should be noted that equation (7,9) in reality does not represent the flux outside the medium. One may also use in this problem the pseudoboundary condition \(nv=0\), if we introduce a fictitious boundary of the medium at \(x=t'=t+d\). Since the diffusion equation for the neutron flux (7,1) is invalid at distances of several mean free paths from the surface, the value of the flux calculated by means of this equation will be inaccurate. However, it turns out that calculations of critical

of the dimensions of a boiler without a shell, by means of the condition that the neutron flux vanish at a distance \(d=0.71\lambda_t\) from the boundary, proves to be sufficiently accurate. In what follows it is convenient to use the boundary condition \(nv=0\). When using this condition, however, it is necessary to remember that all dimensions must be corrected to take account of the extrapolated boundary \((7,11)\).

The solution for case 3 is now written down directly. It is equal to a linear combination of the functions \((7,3)\) with arbitrary constants \(a\) and \(b\), chosen so as to satisfy the boundary conditions. This gives:

\[ I_0=J(0)=-\frac{\lambda_t}{3}\left[\frac{d}{dx}(nv)\right] =-\frac{\lambda_t}{3}(-aK+bK), \]

\[ nv(t')=ae^{-Kt'}+be^{Kt'}=0, \]

where \(t'\) is the extrapolated thickness. We then obtain, solving the equation for \(b\),

\[ b=-ae^{-2Kt'}. \]

Substituting this value of \(b\), we obtain:

\[ I_0=\frac{\lambda_tKa}{3}\left(1+e^{-2Kt'}\right) \quad \text{or} \quad a=\frac{3I_0}{\lambda_tK}\frac{1}{1+e^{-2Kt'}}. \]

Equation \((7,3)\) takes the form:

\[ nv=ae^{-Kx}+be^{Kx}=a\left(e^{-Kx}-e^{-2Kt'}e^{Kx}\right)= \]

\[ =ae^{-Kt'}\left[e^{K(t'-x)}-e^{-K(t'-x)}\right] =2ae^{-Kt'}\operatorname{sh}K(t'-x), \]

where the hyperbolic sine has been introduced: \(\operatorname{sh}z\equiv \dfrac{e^z-e^{-z}}{2}\). Thus the complete solution is:

\[ nv=\frac{3I_0}{\lambda_tK\,\operatorname{ch}Kt'}\operatorname{sh}K(t'-x). \tag{7,12} \]

As \(t'\to\infty\), this solution tends to \((7,4)\) with \(I_0=\dfrac{Q}{2}\).

Case 4. The solution of the problem of neutron diffusion in a medium consisting of a series of plates of different material can be carried out by the method of case 3. The only difference lies in the boundary conditions imposed on the planes separating two different media. At these boundaries we must require continuity of the neutron flux \(nv(x)\) and of the current density \(J(x)\). The physical meaning of these conditions is obvious. If the neutron flux were not continuous, this would mean that neutrons are being born or destroyed at an infinitely thin boundary; on the other hand, if there is a discontinuity in the value

neutron current, it will be rapidly smoothed out by diffusion of neutrons from places with a large flux to places with a small flux. We can then write, as boundary conditions for each of the boundaries:

\[ \begin{array}{ll} (\mathrm{A}) & nv_{\mathrm{I}}=nv_{\mathrm{II}},\\[6pt] (\mathrm{B}) & -\dfrac{\lambda_{\mathrm{I}}}{3}(nv_{\mathrm{I}})'=-\dfrac{\lambda_{\mathrm{II}}}{3}(nv_{\mathrm{II}})'. \end{array} \tag{7,13} \]

Suppose that it is necessary to solve the problem of finding the neutron flux inside a medium consisting of \(n\) plane layers, under the condition that the total density of neutrons at \(x=0\) is equal to \(I_0\).

Fig. 12.

We shall describe only the essence of the method of solving the problem, without carrying out the calculations in detail. In each of the \(n\) layers (denoted \(1, 2\), etc.) we shall have a solution of the form (7,3):

\[ nv_i=a_i e^{-K_i x}+b_i e^{K_i x} \]

with two arbitrary constants. In all there will be \(2n\) arbitrary constants \((a_1 b_1,\ a_2 b_2,\ \text{etc.})\), which must be determined with the aid of the boundary conditions. We must therefore have, for their determination, \(2n\) independent equations. Since there are two boundary conditions (of types A and B) at each of the \(n-1\) boundaries and one boundary condition at each of the external boundaries, we have just the required number of linear equations. Thus the complete solution, giving the neutron flux at any point inside the medium, is reduced to the solution of \(2n\) linear equations with \(2n\) unknowns.

8. ALBEDO

Not all neutrons which fall, for example, on a graphite block remain in the graphite. Some of them are scattered back in the graphite and emerge through the same face of the block through which they entered. The ratio of the reverse current \(J_+\) to the direct current \(J_-\) is called the albedo \(A\), or the reflection coefficient of the medium. By definition,

\[ A=\frac{J_+}{J_-} = \frac{\dfrac{nv_0}{4}+\dfrac{\lambda_t}{6}(nv)'_0} {\dfrac{nv_0}{4}-\dfrac{\lambda_t}{6}(nv)'_0} = \frac{1+\dfrac{2}{3}\lambda_t\left(\dfrac{nv'}{nv}\right)_0} {1-\dfrac{2}{3}\lambda_t\left(\dfrac{nv'}{nv}\right)_0}, \tag{8,1} \]

where the expressions for \(J_+\) and \(J_-\) are obtained from formula (6,5) if only the first two terms are taken in it and \(\lambda_t\) is substituted for \(\lambda_s\). In this case the albedo will be a function of the free path and of the logarithmic derivative of the neutron flux.

In the one-dimensional case, in an infinitely thick layer, the neutron flux is

\[ nv \sim e^{-Kx}, \]

and consequently the logarithmic derivative is equal to

\[ \left(\frac{nv'}{nv}\right)_0=-K, \]

which gives

\[ A=\frac{1-\frac{2}{3}\lambda_t K}{1+\frac{2}{3}\lambda_t K}. \tag{8,2} \]

In the case of graphite, the albedo for an infinitely thick layer can be determined using the values \(\lambda_t=2.86\) and \(K=0.02\). Then

\[ A=\frac{1-0.0375}{1+0.0375}=0.93. \]

This means that \(93\%\) of the neutrons incident on a thick (infinite) block of graphite are reflected and only \(7\%\) are absorbed by the graphite.

For a layer of finite thickness \(t\), \(nv\) is determined by formula (7,12):

\[ nv=C\,\operatorname{sh} K(t-x);\qquad (nv)'=-KC\,\operatorname{ch} K(t-x). \]

Substitution into formula (8,1) gives:

\[ A=\frac{1-\frac{2}{3}\lambda_t K\,\operatorname{cth} Kt}{1+\frac{2}{3}\lambda_t K\,\operatorname{cth} Kt}. \tag{8,3} \]

For \(Kt\to\infty\), \(\operatorname{cth} Kt\to 1\), and expression (8,3) reduces to the expression for the case of infinite thickness (8,2).

It is evident from formula (8,3) that, in order for the albedo of the medium to be large, the quantity \(K\lambda_t\) must be small. In other words, this means that the mean transport length must be small in comparison with the diffusion length in the medium under consideration.

The albedo defined above depends both on the nuclear constants of the substances composing the pile and on its shape. This is most easily seen in the case of a reflector having spherical symmetry (a cavity). The neutron flux for \(r>R\), determined by formula (7,7), in this case will be:

\[ nv=C\frac{e^{-Kr}}{r}, \]

whence we obtain:

\[ \left(\frac{nu'}{nu}\right)_{r=R}=-K-\frac{1}{R}. \]

Substitution in formula (8.1) gives for the albedo the expression

\[ A=\frac{1-\frac{2}{3}\lambda_t\left(K+\frac{1}{R}\right)} {1+\frac{2}{3}\lambda_t\left(K+\frac{1}{R}\right)}. \tag{8.4} \]

The fact that the albedo for a spherical cavity is smaller than the albedo for a plane layer can be explained as follows. Neutrons which diffuse from the cavity into the medium have a smaller probability of being scattered back into the cavity the faster they diffuse in the medium, since the probability of backscattering (even in the absence of absorption), roughly speaking, is proportional to the solid angle under which the cavity is seen from the point at which the neutron is located.

Fig. 13.

The albedo for a medium is directly related to the mean number of times a given plane situated in the medium is crossed.

In order to diffuse through the boundary, a neutron must cross it an odd number of times. Let the probability that the neutron crosses the boundary 1, 3, 5, ... times be \(p(1)\), \(p(3)\), \(p(5)\), ... Then, if the albedo is equal to \(A\):

\[ p(1)=1-A, \]

\[ p(3)=A(1-A), \]

\[ p(5)=A^2(1-A), \]

\[ \ldots\ldots\ldots \]

or, in the general case:

\[ p(2i+1)=A^i(1-A) \quad \text{for any } i. \]

The sum of the probabilities is equal to:

\[ \sum_{n\ \text{odd}} p(n) =1-A+A(1-A)+A^2(1-A)+\ldots \]

\[ =(1-A)(1+A+A^2+A^3+\ldots) \]

\[ =\frac{1-A}{1-A}=1, \]

where the geometric progression is summed for \(A<1\).

The average number of intersections \(\bar n\) is equal to:

\[ \begin{aligned} \bar n &= \sum n p(n) = p(1) + 3p(3) + 5p(5) + \ldots \\ &= 1 - A + 3A(1 - A) + 5A^2(1 - A) + \ldots \\ &= (1 - A)(1 + 3A + 5A^2 + \ldots) = 1 + 2A + 2A^2 + 2A^3 + \ldots \\ &= 1 + \frac{2A}{1 - A} = \frac{1 + A}{1 - A}. \end{aligned} \]

For graphite (of infinite thickness)

\[ A = 0.93 \quad \text{and} \quad \bar n = \frac{1 + 0.93}{1 - 0.93} = 27 . \]

Thus, it should be expected that a neutron diffusing in a thick layer of graphite intersects the given plane on average 27 times.

9. SPATIAL DISTRIBUTION OF THE MODERATION DENSITY

In Sections 4 and 5, neutron moderation was considered only from the point of view of their loss of energy; the distribution of neutrons in space during the process of moderation was not considered at all. In reality, the moderation density \(q\) is a function not only of the energy \(E\), but also of the coordinates \(x, y, z\). This circumstance proves to be very important for chain reactors, since if a neutron after its birth diffuses too far while remaining fast, then it may not become thermal at all before leaving the reactor, and therefore may prove incapable of inducing a new fission.

The term describing the birth \(Q(E_0)\), where \(E_0\) is the initial energy of the neutron that arose during the fission event, is proportional to the flux of thermal neutrons \(nv_\tau\). We shall assume that for each absorbed neutron there arise, on average, \(k\) new neutrons. Since the number of neutrons absorbed per second in a cubic centimeter is equal to \(nv_\tau \Sigma_a\), the number of neutrons arising per second in a cubic centimeter is equal to \(nv_\tau \Sigma_a k\). If we had a uniform distribution of the thermal-neutron flux in space and no absorption during the process of moderation, then all these neutrons would reach thermal velocities. In that case we would obtain \(Q(E_\tau) = Q(E_0)\). However, this equation is not valid for a real reactor, as a result, first, of absorption in the process of moderation and, second, of leakage in the process of moderation. The first cause was discussed in Section 5; here we shall consider the second. It is not difficult to see that in a reactor with a chain reaction the flux of thermal neutrons is not

constant, but has a higher value at the center of the boiler and a lower one at its boundary. As a result, there arises a continuous flux of neutrons directed from the center to the boundary of the boiler, where they escape outward; therefore the multiplication of neutrons in a chain reactor must be such as to compensate both these losses and the losses due to nonproductive absorption in the boiler itself.

In addition, there will also be losses of fast neutrons, especially those born near the boundaries of the boiler. However, since the flux of thermal neutrons at the boundaries of the boiler is small, the number of fast neutrons arising in this region will likewise be small. In any case, the fluxes of both fast and slow neutrons will be small near the boundaries of the boiler. In the general case, if we examine the spatial distribution of neutrons of various energies arising from a point source, we find that the smaller the neutron energy, the more diffuse their distribution will be in space. Let us give an exact expression of this fact. It is not difficult to see that each point source of fast neutrons (arising in fission) leads to the appearance of a distributed source of thermal neutrons. In the stationary state, the number of neutrons in the energy interval from \(E\) to \(E+dE\), leaving \(1\ \mathrm{cm}^3\) per second, is equal to

Fig. 14.

\[ -\frac{\lambda_t}{3}\Delta n(E)v\,dE. \]

If absorption in the process of slowing down is neglected, then this flux will be equal to \(Q(E)dE\), i.e. to the excess of the neutrons entering the interval \(dE\) over the number of neutrons leaving this interval.

Therefore we have:

\[ Q(E)dE=q(E+dE)-q(E)=\frac{\partial q}{\partial E}dE, \]

where \(q(E)\) is again the slowing-down density. We can now write:

\[ -\frac{\lambda_t}{3}\Delta n(E)v\,dE=\frac{\partial q}{\partial E}\cdot dE \tag{9,1} \]

and, substituting \(n(E)v=q/\xi\Sigma_s E\) (formula (4,7)), we obtain:

\[ \frac{\lambda_t}{3}\frac{1}{\xi\Sigma_s E}\Delta q+\frac{\partial q}{\partial E}=0. \tag{9,2} \]

This equation can be written in a simplified form if we introduce the quantity \(\tau\), defined by the formula

\[ d\tau=\frac{\lambda_t}{3\xi\Sigma_s}\,dE. \tag{9,3} \]

Then equation (9.2) takes the form:

\[ \Delta q+\frac{\partial q}{\partial \tau}=0. \tag{9.4} \]

Let us note that \(q\) is a function of the coordinates and of the energy, which enters the equation through \(\tau\). When integrating equation (9.3), we choose the constant of integration in such a way that for thermal energies \(\tau=0\). Then \(\tau(E_T)=0\). This gives us the formula for determining \(\tau(E)\):

\[ \tau(E)=\int_{E_{th}}^{E}\frac{\lambda_t}{3\xi\Sigma_s}\,\frac{dE}{E} =\int_{E_{th}}^{E}\frac{1}{3\xi\Sigma_s\Sigma_t}\,\frac{dE}{E}. \tag{9.5} \]

The quantity \(\tau\), which plays the principal role in the theory of slowing down, is called the “Fermi age,” or simply the “age.” This name is connected with the fact that this quantity enters equation (9.4) analogously to the way in which time enters the equation of heat conduction, although, of course, it has nothing in common with time.

Equation (9.5) can be integrated if we introduce certain averaged values of \(\lambda_t\) and \(\Sigma_s\). Then:

\[ \tau(E)=\frac{\overline{\lambda_t}}{3}\left[\frac{\ln \dfrac{E}{E_T}}{\xi\overline{\Sigma_s}}\right]. \]

Fig. 15.

As we saw in deriving equation (5.7), the quantity in brackets is the mean distance \(X\) (measured along the neutron’s zigzag path) traversed by the neutron during slowing down from \(E\) to \(E_T\). Therefore \(\tau\) may be written in the form

\[ \tau(E)=\frac{\overline{\lambda_t}}{3}X. \]

It can be seen that \(\tau(E)\) has the dimension of the square of a length and is, in the main, proportional to the logarithm of the neutron energy.

If a noticeable fraction of neutrons is absorbed during slowing down, then equation (9.1) becomes inapplicable. In this case one should add one more term, \(-n(E)v\Sigma_a\,dE\), in order to take into account the absorption occurring in the energy interval \((E,E+dE)\). We then have:

\[ \frac{\lambda_t}{3}\Delta n(E)v-n(E)v\Sigma_a+\frac{\partial q'}{\partial E}=0, \tag{9.6} \]

where \(q'(E)\) is the slowing-down density in the presence of absorption. Substituting:

\[ q'(E)=n(E)v\xi\Sigma_s E \quad\text{and}\quad d\tau=-\frac{\lambda_t}{3\xi\Sigma_s}\,\frac{dE}{E}. \]

We obtain, as before,

\[ \Delta q' - \frac{3\Sigma_a}{\lambda_t} q' + \frac{\partial q'}{\partial \tau}=0. \]

If we multiply this equation by \(e^{+\int_{\tau}^{\tau_0}\frac{3\Sigma_a}{\lambda_t}\,d\tau}\) and set

\[ q(E) \equiv q'(E)e^{+\int_{\tau}^{\tau_0}\frac{3\Sigma_a}{\lambda_t}\,d\tau}, \tag{9,7} \]

then we obtain:

\[ \Delta q + (\partial q/\partial \tau)=0, \]

an equation coinciding with the analogous equation (9,4). The meaning of the substitution (9,7) can be understood if we find the solution for \(q'(E)\) and express the integral in terms of \(E\), and not in terms of \(\tau\). The final expression will have the form:

\[ q'(E)=q(E)e^{-\int_E^{E_0}\frac{\Sigma_a}{\xi\Sigma_s}\frac{dE}{E}}, \tag{9,8} \]

which, accurate to the factor \(\frac{1}{1+\Sigma_a/\Sigma_s}\), under the integral coincides with the more rigorous expression (5,4) obtained by Wigner. From this expression it is seen that the influence of absorption in the slowing-down process on the slowing-down density \(q\) reduces to multiplying this quantity by the factor

\[ p(E)=e^{-\int_E^{E_0}\frac{\Sigma_a}{\xi\Sigma_s}\frac{dE}{E}}, \tag{9,9} \]

which depends only on the energy and does not change the dependence of \(q\) on the coordinates. The quantity \(p(E_T)\) is the probability that the neutron will reach thermal energies without being absorbed (the probability of avoiding resonance capture).

To solve equation (9,4) for \(q(x,y,z,\tau)\), we shall try to seek a solution in the following form:

\[ q=e^{\alpha x+\beta y+\gamma z+\varepsilon \tau}. \]

Substituting this expression into the equation, we find that it is a solution under the condition that:

\[ \alpha^2+\beta^2+\gamma^2+\varepsilon=0. \]

The solution may be written in the form

\[ q=e^{\alpha x+\beta y+\gamma z}\cdot e^{-(\alpha^2+\beta^2+\gamma^2)\tau}, \tag{9,10} \]

where \(\alpha,\beta,\gamma\) are arbitrary. Since equation (9,4) is linear, any linear combination of such solutions is also a solution.

We can now, with the aid of the Fourier integral method, find such a linear combination that will satisfy the prescribed boundary conditions.

In particular, if we have a point source located at the point \((x, y, z) = (0, 0, 0)\) and emitting \(Q_0\) neutrons per second having one and the same energy \(E_0\) (corresponding to \(\tau = \tau_0\)), then the boundary conditions have the form:

\[ q(x, y, z, \tau_0) = Q_0 \delta(x, y, z), \]

where \(\delta(x, y, z)\) is the “delta function,” which is equal to zero everywhere except at the origin and for which \(\int \delta(x, y, z)\,dV = 1\), if the region of integration contains the origin. It can be shown that, under such boundary conditions, the solution of equation (9.4), obtained by superposing (more precisely, integrating with respect to \(\alpha, \beta, \gamma\)) solutions of the form (9.10), will have the form:

\[ q(x, y, z, \tau) = Q_0 \frac{ e^{-\frac{x^2+y^2+z^2}{4(\tau_0-\tau)}} }{ [4\pi(\tau_0-\tau)]^{3/2} } = \frac{Q_0}{[4\pi(\tau_0-\tau)]^{3/2}} e^{-\frac{r^2}{4(\tau_0-\tau)}}, \tag{9.11} \]

where \(r^2 = x^2 + y^2 + z^2\).

For each energy (corresponding to some value of \(\tau\)), \(q(r)\) is represented by a Gaussian error curve, whose width and height depend on \(\tau\) (or on the energy). The width (corresponding to a decrease of the ordinate by a factor of \(e\)) is equal to \(2\sqrt{\tau_0-\tau}\).

In Fig. 16 this function is shown schematically for two values of the energy. As we see, for large energies the \(q(r)\)-distribution has a high and narrow maximum; in the case of small energies this maximum is low and diffuse. This corresponds to our physical picture: neutrons that have undergone very few collisions and have lost only a very small fraction of their initial energy can be detected almost exclusively near the origin; neutrons, however, that have lost a large part of their energy are distributed over a wide region, since on average they have made a considerably larger number of collisions and have traveled a greater distance.

Fig. 16.

The physical meaning of the quantity \(\tau_0\) becomes clearer if we calculate the mean square distance from the source for neutrons,

entering the region of thermal energies (for which \(\tau_0 = 0\)). We obtain:

\[ \overline{r^2} = \frac{\int r^2 q_{th}(r)\,dV}{\int q_{th}(r)\,dV} = \frac{ \int\limits_0^\infty r^3 e^{-\frac{r^2}{4\tau_0}}\,4\pi r^2 dr }{ \int\limits_0^\infty e^{-\frac{r^2}{4\tau_0}}\,4\pi r^2 dr } = 6\tau_0, \tag{9,12} \]

where \(\tau_0\) is the age of fission neutrons. The age of fission neutrons can be measured experimentally in various substances with the aid of a detector made of indium foil coated with cadmium, which becomes radioactive (\(\beta\)-decay with a period of 54 min.) when neutrons with the resonance energy of indium (\(1.4\) ev) fall upon it. The dependence of the activity of the foil on its position in the system gives the neutron distribution, from which the age can be calculated by means of formula (9,12).

Equation (9,4) and the results based on it contain the assumption that the slowing down of neutrons occurs continuously, and not in the form of a discrete process. Although such a description is inaccurate, it is nevertheless suitable for such moderators as graphite, for which \(\xi\) is small. For neutrons slowing down in water, the form of the distribution is far from a Gaussian curve, since even in a single collision with a hydrogen atom a neutron may lose a large part of its initial energy. In this case the picture of continuous slowing down proves unsuitable.

10. SOLUTIONS OF THE REACTOR EQUATION

In a reactor with a chain reaction operating at a constant level of neutron density, there must exist a balance between the production of neutrons, on the one hand, and leakage and absorption, on the other. This condition must be satisfied for every element of the reactor volume and for every energy interval \(dE\).

Thus, we must solve an infinite number of equations, which may be written in the form

\[ \frac{\lambda_t}{3}\Delta n(E)v - n(E)v\Sigma_a + Q(E)=0 \tag{10,1} \]

for all \(E\) from \(E_0\) to \(E_t\).

It is important not to confuse the slowing-down density \(q(E)\) and \(Q(E)dE\)—the number of neutrons (per second per \(\text{cm}^3\)) which compensates leakage and absorption in the energy interval \(dE\). As was shown above, in deriving equation (10,1),

$$ Q(E)\,dE=dq=\frac{\partial q}{\partial E}\,dE. $$

When calculating \(Q(E_T)\), it must be remembered that when the neutrons reach thermal energies, they subsequently cease to slow down. Therefore, in this case the number of neutrons passing in \(1\) sec through \(1\ \mathrm{cm}^3\) into thermal energies is simply equal to the slowing-down density calculated for an energy slightly exceeding the thermal one:

$$ Q(E_T)=q(E_T). \tag{10,2} $$

In the presence of absorption during the slowing-down process, we must modify this expression by means of formula (9,9):

$$ Q(E_T)=q'(E_T)=q(E_T)p(E_T), \tag{10,2a} $$

where \(p\) is the probability with which the neutron avoids resonance capture.

It seems convenient to collect all neutrons of thermal energy into one group, even if they have somewhat different energies. We can then write a separate equation for the flux of thermal neutrons \(nv_T\):

$$ \frac{\lambda_{\mathrm{tr}}}{3}\Delta nv_T-nv_T\Sigma_{aT}+Q(E_T)=0. \tag{10,3} $$

For the flux of fast neutrons we use equation (9,4):

$$ \Delta q+\frac{\partial q}{\partial \tau}=0, $$

equivalent to equation (10,1), but expressed in terms of the slowing-down density \(q(E,x,y,z)\).

For the time being we still do not know \(Q(E_T)=q(E_T)\), which must be substituted into equation (10,3). It must appear as the solution of equation (9,4). We can, however, obtain \(q(E_0)\) in the following way. If we know the fission cross section \(\Sigma_f\), the total absorption cross section \(\Sigma_a\) for the boiler material, and \(\nu\)—the average number of neutrons emitted per fission event—then we can calculate the multiplication constant \(k_T\), defined as the number of fission neutrons appearing for each absorbed thermal neutron. The usual definition of the multiplication constant differs slightly from this; it is as follows:

If the process begins with the absorption of one thermal neutron, then the probability that it will produce fission is equal to \(\Sigma_f/\Sigma_a\). We then obtain \(\nu\Sigma_f/\Sigma_a\) fast neutrons for each absorbed thermal neutron. Therefore

$$ k_T=\frac{\Sigma_f}{\Sigma_a}\cdot \nu. $$

We can further write, for the density of slowing-down at the initial energy with which neutrons are produced in the fission process, \(E_0\) (age \(\tau_0\)), the expression

\[ q(\tau_0)=q(E_0)=nv_{\tau}\Sigma_{a_\tau}k_\tau. \tag{10,4} \]

Thus, our problem consists in solving the stationary boiler equations, which we rewrite once more:

\[ \frac{\lambda_\tau}{3}\Delta nv_\tau - nv_\tau\Sigma_{a_\tau}+Q_{th}=0, \tag{10,3} \]

\[ \Delta q+\frac{\partial q}{\partial \tau}=0. \tag{9,4} \]

These equations must be solved under the following boundary conditions:

\[ nv_\tau=0 \text{ at the boundary,} \tag{10,5} \]

\[ q(E)=0; \]

\[ q(E_0)=q(\tau_0)=nv_\tau\Sigma_{a_\tau}\cdot k_\tau. \tag{10,4} \]

Suppose that

\[ q(E)=nv_\tau\Sigma_{a_\tau}\cdot k_\tau f(\tau), \]

where \(f(\tau)\) is a function of \(\tau\), on which we impose the condition \(f(\tau_0)=1\), so that equation (10,4) is satisfied for \(\tau=\tau_0\). Substituting this expression into equation (9,4), we obtain, after division by \(nv_\tau\Sigma_{a_\tau}k_\tau f\):

\[ \frac{\Delta nv_\tau}{nv_\tau}-\frac{f'(\tau)}{f(\tau)}=0. \tag{10,6} \]

Since the first term is a function only of \(x, y,\) and \(z\), while the second term is a function only of \(\tau\), the equation is satisfied only in the case when each of its terms is equal to a constant.

We shall therefore write:

\[ \frac{\Delta nv_\tau}{nv_\tau}=A \tag{10,7} \]

and

\[ \frac{f'(\tau)}{f(\tau)}=-A, \tag{10,8} \]

where \(A\) is a constant. Integrating equation (10,8) from \(\tau\) to \(\tau_0\), we obtain:

\[ f(\tau)=f(\tau_0)e^{A(\tau_0-\tau)}=e^{A(\tau_0-\tau)}, \tag{10,9} \]

since \(f(\tau_0)=1\). Passing then from \(f\) back to the expression for \(q(E)\), we find:

\[ q(E)=nv_\tau\Sigma_{a_\tau}k_\tau e^{A(\tau_0-\tau)}. \tag{10,10} \]

Thus, the term describing the production of thermal neutrons takes the form:

\[ Q(E_\tau)=q(E_\tau)\cdot p = nv_\tau \Sigma_{a_\tau} k_\tau \cdot p e^{A\tau_0}, \tag{10,11} \]

where \(\tau=0\) corresponds to \(E=E_\tau\).

The product \(k_\tau\cdot p \equiv k\) is what is usually called the multiplication constant, denoting the number of thermal neutrons appearing per absorbed thermal neutron; it is assumed here that during slowing down there is no leakage of fast neutrons. In reality, the number of fast neutrons appearing per absorbed thermal neutron is equal to \(k_\tau \cdot p e^{\tau_0 A}\), where the factor \(e^{\tau_0 A}\) determines the fraction of fast neutrons that have avoided leakage.

Let us now substitute this expression into equation (10,3). After using equation (10,7), we obtain:

\[ \frac{\lambda_\tau}{3} A n v_\tau - nv_\tau \Sigma_{a_\tau} + nv_\tau \Sigma_{a_\tau} k e^{\tau_0 A}=0. \tag{10,12} \]

In this equation the function \(nv_\tau\) enters as a factor in each of the terms. If we divide by \(nv_\tau \Sigma_{a_\tau}\) and substitute

\[ L^2=\frac{\lambda_\tau}{3\Sigma_{a_\tau}} \]

(formula 7,5), then after regrouping the terms we obtain:

\[ k-\left(1-L^2 A\right)e^{-\tau_0 A}=0. \tag{10,13} \]

This transcendental equation for \(A\) is called the critical equation. The quantity \(A\) is a function only of the reactor constants \(k, L^2, \tau_0\), each of which, in principle, can be computed from nuclear constants. In order to emphasize that this quantity is a property of the materials forming the reactor, and does not depend on the shape or size of the reactor, we shall denote it by \(A_m\), where the subscript “\(m\)” refers to the reactor material. \(A_m\) is defined as the solution of equation (10,13).

If the multiplication constant \(k\) is greater than unity, then equation (10,13) has one real negative solution, as can be seen by elementary considerations.

It now remains to find the solution of the equation

\[ \Delta nv_\tau = A nv_\tau, \tag{10,7} \]

satisfying the boundary conditions \(nv_\tau=0\) on the extrapolated boundary of the reactor. Let us note that this is, in essence, the same equation whose solutions (for \(a^2<0\)) were given in Table II of Section 7.

For given dimensions and shape of the reactor, the solution of equation (10,7), which vanishes on the boundary, corresponds to a definite

to the value \(A\) found. This quantity depends only on the shape and size of the boiler. We shall denote this quantity by \(A_{\mathrm{g}}\) (the subscript “g” indicates that \(A\) depends on the geometry). The criticality condition of the boiler is now written in the form:

\[ A_{\mathrm{g}}=A_{\mathrm{m}} . \]

Only in this case is \(nv_r\) a solution of the stationary boiler equation. \(A_{\mathrm{g}}\) and \(A_{\mathrm{m}}\) are defined so that below we shall be able to consider, in general form, the kinetics of the boiler. If the power level of the boiler increases, we say that the boiler is supercritical, and then \(|A_{\mathrm{g}}|<|A_{\mathrm{m}}|\). If the power of the boiler decreases, the boiler is subcritical, and \(|A_{\mathrm{g}}|>|A_{\mathrm{m}}|\). The quantity \(-A_{\mathrm{g}}\) basically determines the curvature of the neutron-flux distribution. If this curvature is too large (\(-A_{\mathrm{g}}>-A_{\mathrm{m}}\)), which corresponds to a very small boiler, the boiler will be subcritical.

In the following section we shall give explicit expressions for \(A_{\mathrm{g}}\), as functions of the boiler dimensions, for boilers of various shapes.

11. CRITICAL DIMENSIONS FOR BOILERS OF THE SIMPLEST SHAPES

In this paragraph we shall present solutions of the equation

\[ \Delta nv = A nv, \tag{10,7} \]

satisfying the condition \(nv=0\) on the boundary.

The solutions are taken directly from Table II.

In addition, we also give the critical dimensions for boilers of the simplest shapes.

I. Plane boiler (a layer infinite in the directions of the axes \(y\) and \(z\); thickness \(t\) in the \(x\) direction).

Boundary:

\[ \text{plane } x=\pm \frac{1}{2}t. \]

Equation:

\[ \frac{d^2 nv}{dx^2}=A_{\mathrm{g}}\,nv. \]

Solution:

\[ nv=C\cos\frac{\pi x}{t}; \qquad A_{\mathrm{g}}=-\frac{\pi^2}{t^2}. \]

Critical dimension:

\[ A_{\mathrm{g}}=A_{\mathrm{m}}\longrightarrow t_0=\frac{\pi}{\sqrt{-A_{\mathrm{m}}}}. \]

Critical volume:

\[ V_0\sim\infty. \]

II. Spherical boiler (radius \(R\)).

Boundary:

sphere of radius \(R\).

Equation:

\[ \frac{d^2 n v}{d r^2}+\frac{2}{r}\frac{d n v}{d r}=A_r n v. \]

Solution:

\[ n v=\frac{C}{r}\sin\frac{\pi r}{R};\qquad A_r=-\frac{\pi^2}{R^2}. \]

Critical size:

\[ A_r=A_m \to R_0=\frac{\pi}{\sqrt{-A_m}}. \]

Critical volume:

\[ V_0=\frac{4\pi}{3}R_0^3=\frac{4\pi^4}{3(A_m)^{3/2}}=\frac{130}{(-A_m)^{3/2}}. \]

III. Cylindrical boiler (radius \(R\), height \(H\)).

Boundary: circular cylinder \(r=R\), planes \(z=\pm \frac{1}{2}H\).

Equation

\[ \frac{\partial^2 n v}{\partial r^2}+\frac{1}{r}\frac{\partial n v}{\partial r}+\frac{\partial^2 n v}{\partial z^2}=A_r n v. \]

Solution:

\[ n v=C\cos\frac{\pi z}{H}J_0\left(\frac{2.405\,r}{R}\right);\quad A_r=-\left(\frac{\pi^2}{H^2}+\frac{(2.405)^2}{R^2}\right). \]

Critical dimensions corresponding to the smallest volume:

\[ A_r=A_m \to \begin{cases} R_0=\dfrac{2.405\sqrt{\dfrac{3}{2}}}{\sqrt{-A_m}}=\dfrac{2.495}{\sqrt{-A_m}},\\[1.2em] H_0=\dfrac{\pi\sqrt{2}}{2.405}R_0=1.847R_0, \end{cases} \qquad \begin{cases} H_0=\dfrac{\pi\sqrt{3}}{\sqrt{-A_m}}=\dfrac{5.441}{\sqrt{-A_m}}. \end{cases} \]

Smallest critical volume:

\[ V_0=\pi R_0^2H_0=\frac{148.2}{(-A_m)^{3/2}}. \]

IV. Boiler in the form of a rectangular parallelepiped—edges \(a\), \(b\), and \(c\) in the directions \(x\), \(y\), \(z\).

Boundary:

planes \(x=\pm \frac{1}{4}a,\ y=\pm \frac{1}{2}b,\ z=\pm \frac{1}{2}c\).

Equation:

\[ \frac{\partial^2 n v}{\partial x^2}+\frac{\partial^2 n v}{\partial y^2}+\frac{\partial^2 n v}{\partial z^2}=A_r n v. \]

Solution:

\[ nv=C\cos\frac{\pi x}{a}\cos\frac{\pi y}{b}\cos\frac{\pi z}{c}; \qquad A_r=-\left(\frac{\pi^2}{a^2}+\frac{\pi^2}{b^2}+\frac{\pi^2}{c^2}\right). \]

The critical dimensions corresponding to the smallest volume are:

\[ A_r=A_m \to \]

\[ a_0=b=c=\frac{\pi\sqrt{3}}{\sqrt{-A_m}}=\frac{5.44}{\sqrt{-A_m}}. \]

The smallest critical volume is:

\[ V_0=a_0^3=\frac{161}{(-A_m)^{3/2}}. \]

Suppose, for example, that the value of \(A_m\), obtained by solving the critical equation (10.13), is equal to:

\[ A_m=-10^{-4}\ \mathrm{cm}^{-2}. \]

Then the critical dimensions will be:

For	Critical dimensions	Volume
plane boiler	\(t=314\ \mathrm{cm}\)
spherical boiler	\(R_0=314\ \mathrm{cm}\)	\(1.30\cdot10^8\ \mathrm{cm}^3\)
cylindrical boiler of smallest volume	\(R_0=294\ \mathrm{cm}\), \(H_0=544\ \mathrm{cm}\)	\(1.48\cdot10^8\ \mathrm{cm}^3\) (14% greater than for the sphere)
cubic boiler	\(a_0=544\ \mathrm{cm}\)	\(1.61\cdot10^8\ \mathrm{cm}^3\) (24% greater than for the sphere)

The relative dimensions of these boilers are shown to scale in Fig. 17.

Fig. 17.

In Fig. 18 are shown graphs of three functions used in the theory of boilers. These functions are as follows:

\[ \cos(\pi/2)u, \]

where

\[ u= \begin{cases} 2x/t & \text{for a slab},\\ 2z/H & \text{for a cylinder}. \end{cases} \]

\[ (1/\pi u)\sin\pi u,\quad \text{where } u=r/R \text{ for a sphere}, \]

\[ J_0(2.405\,u),\quad \text{where } u=r/R \text{ for a cylinder}. \]

It should be noted that the curves representing these three functions differ very little from one another; in some rough calculations one may, without allowing a large error, use in all three cases one and the same formula with a cosine.

12. NEUTRON CYCLE

In a critical pile the production of neutrons is exactly balanced by absorption and leakage. Suppose that we introduce \(N\) fission neutrons (energy \(E_0\)) into such a pile and follow them during one cycle. These neutrons will

Fig. 18.

slow down; some of them will leave the pile while still fast, while others will be absorbed at energies above thermal. The remaining neutrons will become thermal inside the pile.

The number of neutrons that have reached the energy \(E + dE\) inside the pile is equal to:

\[ N p(E_0, E + dE)\, e^{[\tau_0-\tau(E+dE)]A}, \]

whereas the number of neutrons that have reached the energy \(E\) is equal to:

\[ N p(E_0, E)\, e^{[\tau_0-\tau(E)]A}. \]

The difference of these two expressions represents the loss of neutrons in the energy interval \(dE\), due to absorption and leakage; it is equal to:

\[ N\left[e^{[\tau_0-\tau(E)]A}\frac{dp(E_0,E)}{dE}\,dE\right]+ \]

\[ +Np(E_0,E)e^{[\tau_0-\tau(E)]A}\left(-A\frac{d\tau(E)}{dE}\,dE\right), \]

where the first term represents the number of initial neutrons \(N\) that are absorbed in the energy interval \(dE\), and the second—the number of neutrons leaving the pile in the energy interval \(dE\).

The number of neutrons that reach thermal energies is equal to

\[ Ne^{\tau_0 A}p(E_0,E_T). \]

Of these neutrons, a fraction equal to \(\dfrac{1}{1-L^2A}\) is absorbed while thermal, whereas the remaining fraction \(\dfrac{-L^2A}{1-L^2A}\) leaves the pile, also while thermal.

This can be seen from the following simple reasoning: for any elementary volume of the pile (and thus also for the pile as a whole), the ratio of thermal leakage to thermal absorption is:

\[ \frac{L}{A} = \frac{-\dfrac{\lambda_t}{3}\Delta n v}{nv\Sigma_a} = \frac{-\lambda_t A}{3\Sigma_a} = -L^2A. \tag{12,1} \]

Thus the fraction of absorbed thermal neutrons will be equal to:

\[ \frac{A}{L+A} = \frac{1}{1+\dfrac{L}{A}} = \frac{1}{1-L^2A}, \]

whereas the fraction of neutrons leaving the pile is equal to:

\[ \frac{L}{L+A} = \frac{\dfrac{L}{A}}{1+\dfrac{L}{A}} = \frac{-L^2A}{1-L^2A}. \]

We therefore know how many neutrons out of the initial number \(N\) are absorbed in the pile and at what energies. The number of neutrons equal to

\[ N\frac{e^{\tau_0 A}p(E_0,E_T)}{1-L^2A}, \]

is absorbed at thermal energies, and the number

\[ N\int_{E_T}^{E_0} e^{[\tau_0-\tau(E)]A}\,\frac{dp(E_0,E)}{dE}\,dE \]

is absorbed at energies above thermal. Let us denote by \(k(E)\) the average number of fission neutrons produced per neutron absorbed at energy \(E\). Then:

\[ k(E)=\frac{\Sigma_f(E)}{\Sigma_a(E)}\nu(E). \]

By \(k_T(E)\) we have denoted the value of \(k(E)\) for thermal energies.

The neutron cycle is thus completed. The final result is as follows. If the process began with \(N\) fission neutrons, then after one cycle we obtain

\[ \left[ \frac{e^{\tau_0 A}p(E_0,E_T)k_T}{1-L^2A} + \int_{E_T}^{E_0} e^{[\tau_0-\tau(E)]A}\, \frac{dp(E,E_0)}{dE}\, k(E)\,dE \right] \]

neutrons. This expression may be written in the form:

\[ N\left[ \frac{e^{\tau_0 A}p_T k_T}{1-L^2A} + e^{[\tau_0-\overline{\tau(E)}]\overline{A}}\,(1-p_T)\,\overline{k(E)} \right], \]

where \(\overline{\tau_0-\tau(E)}\) is the mean “age” from fission energy to the energy at which the neutrons are absorbed (during slowing down), \(\overline{k(E)}\) is the mean multiplication constant for neutrons absorbed at energies above thermal, and \(p_T=p(E_0,E_T)\).

If the reactor is critical, then the number of fission neutrons at the end of the cycle must be the same as at the beginning.

Thus, the critical equation for the general case has the form:

\[ \frac{e^{\tau_0A}p_Tk_T}{1-L^2A} + e^{\overline{[\tau_0-\tau(E)]}\,\overline{A}}\,(1-p_T)\,\overline{k(E)} =1. \tag{12,2} \]

A reactor in which fission occurs only owing to the absorption of thermal neutrons can be realized in two ways: either \(p_T=1\) (neutrons are absorbed only at thermal energies), or \(p_T<1\), but \(\overline{k(E)}\) vanishes. The latter corresponds to the case of natural uranium reactors, in which there is considerable resonance absorption (not leading to fission) by \(U^{238}\) nuclei at energies exceeding thermal energies.

In such a pile the critical equation has the form:

\[ \frac{k e^{-\tau_0 A}}{1-L^2 A}=1 \tag{12,3} \]

or

\[ k-(1-L^2 A)e^{-\tau_0 A}=0, \tag{12,4} \]

where we have denoted:

\[ k=k_T p_T=\frac{\Sigma_f(E_T)}{\Sigma_a(E_T)}\,\nu p_T. \tag{12,5} \]

The multiplication constant \(k\) may also be written in the form

\[ k=f\eta p_T, \tag{12,6} \]

where \(f\) is the ratio of the absorption cross section of the fissile substance at \(E_T\) to the absorption cross section of all materials of the pile at \(E_T\), i.e.

\[ f=\frac{\Sigma_a^{\mathrm{fis}}(E_T)}{\Sigma_a(E_T)} \]

is the fraction of absorbed thermal neutrons accounted for by absorption in the fissile substance, and \(\eta\) is \(\nu\) multiplied by the ratio of the fission cross section of the fissile substance at \(E_T\) to the absorption cross section of the fissile substance at \(E_T\), i.e.

\[ \eta=\frac{\Sigma^{\mathrm{fis}}(E_T)}{\Sigma_a^{\mathrm{fis}}(E_T)}\,\nu \]

is nothing other than the number of fission neutrons produced per thermal neutron absorbed in the fissile material.

In a pile operating not on thermal neutrons \(p_T=0\), i.e. not a single neutron reaches thermal energies; this condition can be realized by creating a high concentration of fissile material, or by making the moderator concentration low, or, finally, by providing both. In this case the critical equation may be written in the following form:

\[ \overline{k(E)}\,e^{[\overline{\tau_0-\tau(E)}]A}=1. \tag{12,7} \]

For any of these cases—thermal, nonthermal, or mixed—the critical size of a pile without a reflector is obtained by setting \(A_r=A_m\), where \(A_m\) is the material root of the critical equation (12,2), (12,4), or (12,7), and is determined by the dimensions and shape of the pile (see Section 11).

If the dimensions of the pile are noncritical, then \(A=A_r\) will be greater or less than \(A_m\), and the number of neutrons at the end of the cycle will differ from the number of neutrons at the beginning. This means that the neutron density in the pile will change with time. If there is a relative increase in the number of neutrons in one cycle,

then for mixed noncritical reactors we obtain:

\[ \frac{e^{\tau_0 A} p_T k_T}{1-L^2A} + e^{[\tau_0-\tau(E)]A}(1-p_T)\,\bar{k}(E) =1+\gamma . \tag{12,8} \]

If \(\gamma>0\), the reactor is supercritical and the neutron density increases; if \(\gamma<0\), the reactor is subcritical, and the neutron density decreases.

For a reactor operating on thermal neutrons, \(\gamma\) is given by the formula:

\[ \gamma=\rho\,\frac{e^{\tau_0 A}}{1-L^2A}, \tag{12,9} \]

where:

\[ \rho\equiv k-(1-L^2A)e^{-\tau_0 A} \tag{12,10} \]

is called the reactivity of the reactor. The reactivity of the reactor is an important quantity in the theory of reactor kinetics.

In the preceding consideration of the neutron cycle we assumed that all fission neutrons arise with one and the same energy \(E_0\). In order to take into account the experimentally obtained continuous spectrum of fission neutrons, it should be assumed that the quantities \(\tau_0\) and \(\tau_0-\tau(E)\) in equations (12,2), (12,7), and (12,8) are quantities averaged over the fission spectrum.

13. REACTOR WITH A REFLECTOR; MULTIGROUP METHOD

Up to now we have considered only homogeneous reactors. It is very useful, however, to surround the reactor with a reflector having a high albedo. Since in this way the leakage of both fast and slow neutrons can be reduced, such a reactor will have a smaller critical volume. It turns out that if a very good reflector is used, the overall size of the critical reactor together with the reflector can be made smaller than the size of a “bare” reactor. Thus, with the aid of a reflector, a considerable saving can be achieved in the initial loading of the reactor with fissile material. An additional gain also occurs due to the flattening of the neutron distribution in a reactor with a reflector. In a “bare” reactor the fissile material near the edge of the reactor is used with very low efficiency, since the neutron flux at that point is very small. Owing to the reflector, the neutron flux is made more uniform over the reactor volume, and therefore the averaged flux may prove to be considerably larger than in a reactor without a reflector. If, for example, the reactor power is limited by the temperature at the center of the reactor, then in a reactor with a reflector a higher power can be obtained. The same result can also be obtained by means of a nonuniform distribution of fissile material over the reactor volume.

Unfortunately, it proved impossible to obtain an exact solution of the pile equations for a pile with a reflector in the interesting case when the properties of the reflector differ from those of the reactor. To solve this problem it is necessary to use approximate methods.

A very general class of methods is known as the multigroup method. In this method the neutrons in the pile and the reflector are distributed into groups with identical (averaged) energies; with respect to each of these neutron groups, the materials of the pile and reflector are characterized by certain averaged properties. This method proves very useful for calculating the action of control rods, which strongly absorb thermal neutrons but have little effect on fast neutrons.

For the neutron flux of each energy group one can write an equation of the form (10,1). The equations turn out to be coupled with one another because neutrons slow down; the disappearance of a neutron from some group is accompanied by its appearance in a group of lower energy. As before, the equations for the groups of highest and lowest energy turn out to be coupled with the fission process.

The simplest theory of this type is the one-group theory, in which we assume that absorption and production occur at one and the same neutron energy. This approximation proves not too poor for a fast-neutron pile, as, for example, for a pile in which both absorption and production occur at neutron energies close to those at which they are produced in fission. For a thermal-neutron pile the theory would correctly describe the phenomenon taking place if the neutrons born in fission were thermal, and not fast. To illustrate this method we shall apply it to the cases of plane and spherical piles with reflectors; in the following paragraph we shall consider these problems by the method of two-group theory. It should be noted that the critical size of a “bare” pile (i.e., a pile without a reflector) cannot be correctly calculated by either of these two methods.

However, it is convenient to “adjust” the pile constants so that the correct critical size would be obtained in both the one-group and the two-group approximations.

We shall begin with equation (6,1), which we write as an equation for the neutron flux \(\varphi = nv\). This quantity inside the pile, characterized by the constants \(\lambda_t^0\), \(\Sigma_a^0\), and \(k\), satisfies the relation:

\[ -\frac{\lambda_t^0}{3} A\varphi + \Sigma_a^0\varphi - k\Sigma_a^0\varphi = 0. \tag{13,1} \]

Within the reflector, characterized by the constants \(\lambda_t^1,\ \Sigma_a^1,\ k=0\) (there is no birth of neutrons), this equation has the form:

\[ -\frac{\lambda_t^1}{3}\Delta\varphi+\Sigma_a^1\varphi=0. \tag{13,2} \]

These equations may be written in the form:

\[ \begin{cases} \Delta\varphi+K_0^2\varphi=0 & \text{(in the boiler)}\\ \Delta\varphi-K_1^2\varphi=0 & \text{(in the reflector)} \end{cases} \tag{13,3} \tag{13,4} \]

where:

\[ \begin{cases} K_0^2=(k-1)\dfrac{3\Sigma_a^0}{\lambda_t^0}=\dfrac{k-1}{L_0^2},\\[6pt] K_1^2=\dfrac{3\Sigma_a^1}{\lambda_a^1}=\dfrac{1}{L_1^2}. \end{cases} \tag{13,5} \tag{13,6} \]

They must be solved under the boundary conditions:

I. \(\varphi=0\) on the outer boundary of the reflector.

II. \(\varphi\) and \(\dfrac{\lambda_t}{3}\varphi'\) must be continuous at the boundary between the reactor and the reflector.

For example, consider a plane boiler of thickness \(2a\), covered on both sides by a reflector of thickness \(t\).

The solutions of equations (13,3) and (13,4) are:

\[ \varphi_0=\cos K_0x\quad (|x|\leq a), \tag{13,7} \]

\[ \varphi_1=c\,\operatorname{sh}[K_1(t+a-x)]\quad (a\leq |x|\leq a+t), \tag{13,8} \]

where boundary condition I is satisfied, since we have chosen as the solution the hyperbolic sine, which vanishes at the boundary of the reflector. Condition II at \(x=A\) gives:

\[ \lambda_t^0K_0\operatorname{tg}K_0a=\lambda_t^1K_1\operatorname{cth}K_1t. \tag{13,9} \]

This transcendental equation relates the critical half-thickness of the boiler to the thickness of the reflector. As \(t\to0\), \(\operatorname{cth}K_1t\to\infty\), \(K_0a=\pi/2\); the critical half-thickness is expressed by the same formula as for the “bare” boiler (see section 11), namely: \(a_0=\pi/2K_0\). These formulas, however, are not identical, since \(K_0\), defined by equation (13,5), does not coincide with \(\sqrt{-A_m}\), defined as the solution of (10,13).

Fig. 19.

Let us introduce the reflector saving \(S=a_0-a\), i.e. the amount by which the critical half-thickness of the boiler is decreased (thanks to the reflector) in comparison with the half-thickness of the “bare” boiler. Obviously, \(S\) depends on the thickness

reflector, as well as on the constants of the boiler and the reflector. From equation (13,9) we obtain directly:

\[ S \equiv a_0-a=\frac{1}{K_0}\left[\frac{\pi}{2}-\operatorname{arctg}\left(\frac{\lambda_t^1 K_1}{\lambda_t^0 K_0}\operatorname{cth} K_1 t\right)\right]. \tag{13,10} \]

In the limiting case of a small reflector thickness \((K_1 t \ll 1)\), equation (13,10) will have the form:

\[ S=a_0-a=\frac{\lambda_t^0}{\lambda_t^1}\,t. \tag{13,11} \]

If the reflector scatters so effectively that the transport length in it is less than the transport length in the boiler, then the economy of the reflector proves to be greater than its thickness. This means that the size of the boiler plus the reflector will be smaller than the dimensions of the “bare” boiler.

Analogous calculations can be carried out for the more realistic case of a spherical boiler of radius \(R\), covered with a reflecting layer of thickness \(t\). The transcendental equation corresponding to equation (13,9) will have the form:

\[ \operatorname{ctg} K_0 R=\frac{1}{K_0 R}\left(1-\frac{\lambda_t^1}{\lambda_t^0}\right)-\frac{\lambda_t^1 K_1}{\lambda_t^0 K_0 t}\operatorname{cth} K_1 t, \tag{13,12} \]

and the critical radius of the “bare” boiler is equal to \(R_0=\pi/K_0\). The economy of the reflector \(S=R_0-R(t)\) in this case can be determined by means of a numerical solution of the equation for specified values of the boiler constants.

These quantities were calculated for spherical and plane boilers having the dimensions given in Section 11. We assume that slowing down in the boiler is accomplished by means of graphite and that the boiler is surrounded by a water reflector. We take the following rounded values:

\[ A=-10^{-4}=-K_0^2,\quad K_0=10^{-2}\ \mathrm{cm}^{-1}. \]

In addition, we use the following values of the constants:

\[ \lambda_t^0=2.7\ \mathrm{cm},\quad K_1=\frac{1}{2.8}\ \mathrm{cm}^{-1},\quad \lambda_t^1=0.4\ \mathrm{cm}. \]

The result, shown in Fig. 20, gives the dependence of the economy of the reflector \(S\) on its thickness in these cases. The general character of both curves is the same. The economy of the reflector reaches “saturation,” lying between 18 and 20 cm for reflector thicknesses exceeding 10 cm. As is seen from the curve, a spherical graphite boiler which, without a reflector, has a critical radius of 314 cm can be made critical at a radius of 294 cm if a water reflector 10 cm thick is used.

For a water reflector of infinite thickness (thickness greater than 10 cm), the saving \(S\) turned out to be equal to:

\[ \begin{array}{ll} \text{for a sphere} & 19.7\ \text{cm}\\ \text{for a plane layer} & 18.8\ \text{cm}\\ \text{for an infinite circular cylinder} & 19.7\ \text{cm} \end{array} \]

The absolute numerical values are very crude; however, the general course of the curves differs not very greatly from the course of the curves determined by a more exact theory. It was shown that the value of the saving \(S\), computed with the aid of one-group theory, is underestimated. More exact calculations by the method

Fig. 20.

of two groups give larger values, since the albedo of a reflector for fast neutrons exceeds the albedo for thermal neutrons.

In the case when the diffusion lengths in the reactor and in the reflector are the same \((\lambda_f^0=\lambda_t^1)\), equations (13, 12) for the sphere and (13, 9) for the plane have the same form; therefore the reflector saving, referred to the half-thickness of the plane reactor and to the radius of the spherical reactor, is one and the same.

14. THE TWO-GROUP METHOD

Owing to the fact that the properties of the nuclei of the boiler material and of the reflector with respect to fast neutrons usually differ greatly from their properties with respect to slow neutrons, calculation by the one-group method cannot give exact results. The next order of approximation consists in the fact that

consider separately the fluxes of fast \(\varphi_f\) and thermal neutrons \(\varphi_\tau\). Despite the fact that fast neutrons correspond to a wide range of energies, we shall combine them approximately into one group and assign to the substances in which they diffuse appropriately chosen averaged characteristics.

Inside the reactor \(\varphi_\tau\) and \(\varphi_f\) satisfy the equations:

\[ \frac{\lambda_f}{3}\Delta\varphi_f-\Sigma_{a_f}\varphi_f+k\Sigma_{a_\tau}\varphi_\tau=0, \tag{14,1} \]

\[ \frac{\lambda_\tau}{3}\Delta\varphi_\tau-\Sigma_{a_\tau}\varphi_\tau+\Sigma_{a_f}\varphi_f=0. \tag{14,2} \]

These are equations of the same type as our original equations. The three terms of each of them describe leakage, absorption, and production of neutrons, respectively. We have introduced the fictitious absorption cross section of fast neutrons \(\Sigma_{a_f}\), which takes into account not the actual absorption of fast neutrons, which we neglect, but the loss of neutrons from the fast group, occurring when, as a result of slowing down, they enter the thermal group. The term \(\Sigma_{a_f}\varphi_f\) therefore represents both the source of thermal neutrons (the number of neutrons which become thermal, referred to one cubic centimeter per second) and the loss in the fast group. The source of fast neutrons (fission neutrons), on the other hand, is proportional to the thermal flux and is equal to \(k\Sigma_{a_\tau}\varphi_\tau\), if \(k\equiv k_\tau p_\tau\) is the number of thermal neutrons arising as a result of the absorption of one thermal neutron, and \(p_\tau=1\).

Inside the reflector \(\varphi_f\) and \(\varphi_\tau\) satisfy analogous equations, with the exception that the multiplication coefficient \(k=0\), since we assume that there are no fissile substances in the reflector. The remaining constants in the equation will, generally speaking, also have other values in the reflector.

It is necessary to choose such values of \(\Sigma_{a_f}\) and \(\lambda_f\) as to obtain the same mean square slowing-down length \(6\tau_0\) as is given by the exact theory. This can be achieved in the following way: we try to solve equation (14,1), assuming that at the point \(r=0\) there is a point source of neutrons (fission). Let \(K_f^2=3\Sigma_{a_f}/\lambda_f\); then the equation can be written in the form:

\[ \Delta\varphi_f-K_f^2\varphi=0\quad (r>0) \tag{14,3} \]

and the solution will be (see Section 7):

\[ \varphi_f=A\frac{e^{-K_f r}}{r}, \tag{14,4} \]

For $\overline{r^2}$ we obtain, as in Section 7:

\[ \overline{r^2} = \frac{\int r^3 \varphi_f\, dV}{\int \varphi_f\, dV} = \frac{\int_0^\infty r^2 \frac{e^{-K_f r}}{r}\,4\pi r^2 dr} {\int_0^\infty \frac{e^{-K_f r}}{r}\,4\pi r^2 dr} = \frac{6}{K_f^2}. \tag{14,5} \]

Therefore we choose $K_f^2=\dfrac{1}{\tau_0}$; with this choice, equation (14,5) agrees with equation (9,12). The approximation consists in substituting the function $\dfrac{1}{r}e^{-K_f r}$ for the exact Gaussian function

\[ \frac{e^{-r^2/4\tau_0}}{\sqrt{4\pi\tau_0}}. \]

Another way of solving the problem consists in a detailed consideration of the slowing-down process. The mean number of collisions experienced by a neutron during slowing down from $E_0$ to $E_T$ is equal to
$N=\dfrac{1}{\xi}\ln\dfrac{E_0}{E_T}$. Then, on average, one act of “absorption” will occur for every $N$ collisions. The ratio of the fictitious absorption cross section $\Sigma_{af}$ to the averaged scattering cross section $\overline{\Sigma}_s$ will therefore be equal to

\[ \frac{\Sigma_{af}}{\overline{\Sigma}_s} = \frac{1}{N} = \frac{\xi}{\ln \dfrac{E_0}{E_T}}, \]

whence we obtain:

\[ \Sigma_{af} = \frac{\xi \overline{\Sigma}_s}{\ln \dfrac{E_0}{E_T}}. \tag{14,6} \]

Then:

\[ K_f^2 = \frac{3\Sigma_{af}}{\lambda_f} = \frac{3\xi\overline{\Sigma}_s}{\lambda_f\ln \dfrac{E_0}{E_T}} = \frac{1}{\tau_0}, \tag{14,7} \]

if $\lambda_f$ is defined as the correspondingly averaged transport length for fast neutrons.

Now we shall find solutions of equations (14,1) and (14,2) satisfying the conditions:

\[ \left. \begin{aligned} \Delta\varphi_f &= A\varphi_f,\\ \Delta\varphi_T &= A\varphi_T, \end{aligned} \right\} \tag{14,8} \]

where in both equations the constant $A$ is one and the same. After sub-

substitution into equations (14.1) and (14.2), we obtain a system of two linear equations:

\[ \left(\frac{\lambda_f}{3}A-\Sigma_{a_f}\right)\varphi_f+k\Sigma_{a_T}\varphi_T=0, \tag{14,9} \]

\[ \Sigma_{a_f}\varphi_f+\left(\frac{\lambda_T}{3}A-\Sigma_{a_T}\right)\varphi_T=0. \tag{14,10} \]

This system has solutions not equal to zero only in the case when the determinant formed from its coefficients vanishes, i.e., if

\[ \left(\frac{\lambda_f}{3}A-\Sigma_{a_f}\right)\left(\frac{\lambda_T}{3}A-\Sigma_{a_T}\right)-k\Sigma_{a_T}\Sigma_{a_f}=0. \tag{14,11} \]

Dividing by \(\Sigma_{a_T}\Sigma_{a_f}\) and introducing \(L^2=\lambda_T/3\Sigma_{a_T}\) and \(\tau_0=\lambda_f/3\Sigma_{a_f}\), we obtain:

\[ k-(1-L^2A)(1-\tau_0 A)=0. \tag{14,12} \]

This is the critical equation of the two-group theory. It is analogous to equation (10,13) and reduces to this equation if we put \(e^{-\tau_0 A}=1-\tau_0 A\); this is a good approximation if \(\tau_0 A\ll 1\), i.e., if the mean square slowing-down length is small in comparison with the cross-sectional area of the critical boiler. It differs from equation (10,13) in that the equation is quadratic in \(A\), and not transcendental.

For the case when the condition \(\tau_0 A\ll 1\) is not satisfied, the critical size of the “bare” boiler, determined by equation (14,12), is too small (\(A\) is too large); however, one can change the definition of the age \(\tau_0\) in the theory so as to obtain the correct critical size by means of equation (14,12).

In the limit, when the number of groups increases without bound, the critical equation for the multigroup theory passes into the equation of the continuous theory (10,13). This can be shown in the following way. The critical equation for \(n\) groups (in addition to the thermal group) can be written, by analogy with equation (14,12), in the following form:

\[ k-(1-L^2A)(1-\tau_1 A)(1-\tau_2 A)\ldots(1-\tau_n A)=0, \tag{14,12a} \]

where \(\tau_1,\tau_2,\ldots\) are to be understood as the “partial ages” of the neutrons in the various groups. It is clear that, as the neutron is continuously slowed down from fission energy to thermal energies, it will pass through all these groups, and we shall have:

\[ \tau_0=\tau_1+\tau_2+\tau_3+\ldots+\tau_n=\sum_1^n \tau_i. \]

To compute the product

\[ f_n=(1-\tau_1 A)(1-\tau_2 A)(1-\tau_3 A)\ldots(1-\tau_n A). \]

as \(n \to \infty\), take the logarithm of both sides; then we obtain:

\[ \ln f_n=\sum_1^n \ln (1-\tau_i A). \]

Since \(n \to \infty\) and each of the \(\tau_i\) becomes small, we may use the approximate equality: \(\ln(1-x)\simeq -x\). Then

\[ \ln f_\infty=\lim_{n\to\infty}\ln f_n =\lim_{n\to\infty}\sum_{i=1}^n(-\tau_i A) =-A\sum_{i=1}^n \tau_i =-A\tau_0, \]

\[ f_\infty=e^{-\tau_0 A}. \]

Equation \((14,12a)\) then takes the form of the critical equation of the continuous theory, namely:

\[ k-(1-L^2 A)e^{-\tau_0 A}=0. \]

The critical equation may be written in the following form:

\[ A^2-\left(\frac{1}{\tau_0}+\frac{1}{L^2}\right)A-\frac{k-1}{L^2\tau_0}=0, \tag{14,13} \]

Its two solutions will be:

\[ \left. \begin{aligned} A_1&=\frac12\left(-\frac{1}{\tau_0}+\frac{1}{L^2}\right) -\frac12\sqrt{\left(\frac{1}{\tau_0}+\frac{1}{L^2}\right)^2+4\frac{k-1}{L^2\tau_0}},\\ A_2&=\frac12\left(-\frac{1}{\tau_0}+\frac{1}{L^2}\right) +\frac12\sqrt{\left(\frac{1}{\tau_0}+\frac{1}{L^2}\right)^2+4\frac{k-1}{L^2\tau_0}}. \end{aligned} \right\} \tag{14,14} \]

If \(k>1\), the solutions will have opposite signs. The negative root \(A_1\) is approximately equal to

\[ A_1\simeq -\,\frac{k-1}{L^2+\tau_0}, \tag{14,15} \]

and the positive root, if \(k-1\) is small,

\[ A_2\simeq \frac{1}{\tau_0}+\frac{1}{L^2}. \tag{14,15a} \]

For each of the two possible values of \(A\), the ratio \(\varphi_\tau/\varphi_f\) is determined by equations \((14,9)\) and \((14,10)\). Let \(s_1\) denote the ratio corresponding to \(A_1\), and \(s_2\) the ratio corresponding to \(A_2\). Then, if we denote \(K_f^2=1/\tau_0\), the expressions for \(s_1\) and \(s_2\) will have the form:

\[ s_1=-\frac{\lambda_f}{\lambda_\tau}\frac{A_1-K_f^2}{kK_\tau^2}, \qquad s_2=\frac{\lambda_f}{\lambda_\tau}\frac{A_2-K_\tau^2}{kK_\tau^2}. \tag{14,16} \]

If \(A_1<0\) and \(A_2>0\), then \(s_1>0\) and \(s_2<0\). The general solution of equations (14.1) and (14.2) for the fluxes of fast and thermal neutrons is a linear combination of terms corresponding to the two values \(A_1, A_2\), determined by equation (14.14).

In the case of a plane boiler we have, inside it, the solutions:

\[ \left. \begin{aligned} \varphi_f &= \cos \mu_1 x + C \operatorname{ch} \mu_2 x,\\ \varphi_\tau &= s_1 \cos \mu_1 x + C s_2 \operatorname{ch} \mu_2 x, \end{aligned} \right\} \tag{14.17} \]

where the positive quantities \(\mu_1^2=-A_1,\ \mu_2^2=A_2\) have been introduced. The equations for the fluxes of fast and thermal neutrons inside the reflector can be written in the same form as equations (14.1) and (14.2), only with \(k=0\). The critical equation then coincides with equation (14.12) for \(k=0\). Therefore

\[ (1-L'^2 A')(1-\tau_0' A')=0, \tag{14.18} \]

where the primes indicate that these quantities refer to the reflector. The solutions of equation (14.18) will be:

\[ \left. \begin{aligned} A_1' &= \frac{1}{\tau_0'} = K_f'^2,\\ A_2' &= \frac{1}{L'^2} = K_\tau'^2. \end{aligned} \right\} \tag{14.19} \]

In this case the ratios of the fluxes of thermal and fast neutrons corresponding to the two values of \(A\) are:

\[ r_1 = -\frac{\Sigma_{af}'}{\dfrac{\lambda_\tau'}{3}A_1' - \Sigma_a'} = \frac{\lambda_f'}{\lambda_\tau'} \frac{K_f'^2}{K_f'^2-K_\tau'^2}, \]

\[ r_2=\infty. \tag{14.20} \]

The last condition simply means that, in the reflector, the fast-neutron flux corresponding to the second solution vanishes. Then, for a reflector of thickness \(t\), the solution vanishing at the outer boundary \(x=\pm(a+t)\) will be:

\[ \varphi_f' = A\,\operatorname{sh} K_f'(a+t-|x|), \tag{14.21} \]

\[ \varphi_\tau' = r_1 A\,\operatorname{sh} K_f'(a+t-|x|) + B\,\operatorname{sh} K_\tau'(a+t-|x|). \tag{14.22} \]

There remain four boundary conditions that must be satisfied at the interface between the reactor and the reflector at \(x=\pm a\). They express the continuity of the fluxes of fast and thermal neutrons and the continuity of the total current of fast and thermal neutrons. The four equations arising from these four

of the boundary conditions are, in the general case, incompatible, since we have only three arbitrary constants \(A, B, C\). These equations have solutions only for some definite value of the reactor half-width \(a\). Thus, the critical size of the reactor \(a\) is determined as an implicit function of the reflector thickness \(t\), which enters the equations as a parameter. In order to obtain \(a\) as a function of \(t\), it is necessary to resort to rather lengthy numerical calculations. This has been done for certain problems of practical interest, and may be found in the reports connected with the project. The results obtained for definite combinations of reactor and reflector show that the economy of the reflector calculated by the two-group method usually proves to be several percent higher than according to the simplified one-group theory.

15. CONTROL OF THE PILE

Under practical conditions the dimensions of the pile always exceed the critical dimensions, in order to make it possible to control the pile. In thermal-neutron piles, for this purpose use is made of control rods made of a material having a large absorption cross section for thermal neutrons (cadmium, boron steel). Thus, instead of changing the dimensions of the pile, one changes its effective constants \(k\), \(L^2\), and \(\tau_0\) by introducing foreign substances into the construction of the pile or removing them from it.

Another way of considering the problem consists in assuming that the control rods change the boundary conditions which the neutron flux must satisfy. For example, if a control rod is a “black” body with respect to thermal neutrons, then the neutron flux must vanish at the extrapolated boundary of the control rod.

It is convenient to introduce again the quantity \(\rho\) (equation (12,10)), called reactivity, which is defined in the continuous theory by the formula

\[ \rho = k - (1 - L^2 A)e^{-\tau_0 A} \tag{15,1} \]

or by the analogous expression in the one-group method:

\[ \rho = k - (1 - L^2 A). \tag{15,1a} \]

If \(\rho = 0\), then the pile is exactly critical, since equation (10,13) is then satisfied. If \(\rho\) is positive, then the neutron flux increases with time; if \(\rho\) is negative, then it decreases.

The pile and the control rods must be designed so that there is a considerable excess of reactivity, which

could be controlled. As the fissionable materials present in the reactor are consumed, the reactivity decreases, because the number of fissioning atoms that absorb neutrons decreases, and also because the fragments formed in fission absorb neutrons.

Thus useful absorption is reduced and nonproductive absorption is increased. In order for the reactor to continue operating after appreciable use of the fissioning nuclei and the accumulation of harmful fission products, it is necessary to provide a sufficiently large excess of reactivity in the reactor so that the critical conditions \((\rho = 0)\) can be restored by removing the control rods from the reactor.

Another effect is connected with the fact that the reactor constants, and consequently \(A_m\), depend on the temperature of the reactor. There are several temperature effects, which are in principle associated with 1) changes in the absorption and scattering cross sections of the substances in the reactor with neutron energy and 2) changes in the physical constants (dimensions, density) of the reactor with temperature.

The reactivity, depending on the design of the reactor, may either increase with temperature or decrease. It seems desirable to have a reactor with a negative temperature coefficient, since in this case a small disturbance that increases the temperature will lead to a decrease in reactivity, which in turn will lead to a return to the initial temperature. Such a reactor will be stable, whereas a reactor with a positive temperature coefficient will be unstable. Thus, for all reactors the temperature coefficient is an important factor in their stability. In high-temperature reactors this factor is still more important, since in such a reactor the difference between the reactivity of the hot and the cold reactor is relatively large.

In an experimental reactor it is necessary to have an excess of reactivity so that the insertion of neutron absorbers (for example, for the production of radioactive isotopes) would not shut down the reactor.

In addition to these relatively large effects, which change the reactivity of the reactor and require a coarse regulating mechanism, there also exist fluctuations of reactivity with short periods, which must also be regulated in order for the reactor to operate at a constant power level. For example, the reactivity of a reactor may change somewhat because of a change in the temperature of the coolant.

The excess reactivity provided for in the design of the reactor may exceed several times the reactivity associated with delayed neutrons (see Section 16).

Generally speaking, it may be asserted that control of reactors operating on thermal neutrons is simpler, since in them there are nuclei with a large absorption cross section—

cesses for thermal neutrons. The problem of controlling boilers operating with fast neutrons is considerably more difficult, since the absorption cross sections for fast neutrons are always significantly smaller.

In order to clarify for ourselves the influence of a change in the reactivity of the boiler on the neutron flux, we must consider solutions of the boiler equations that depend on time. This will be done in the following section.

16. BOILER EQUATIONS DEPENDING ON TIME

In deriving the boiler equations in Section 10 we made the assumption that there is an exact balance between neutron leakage and absorption, on the one hand, and the production of neutrons, on the other. We must now investigate the case in which there is no such balance. If more neutrons arise in the boiler than are lost through leakage and absorption, then we should expect the neutron flux to increase with time. Indeed, if we carry out this reasoning for \(1\ \mathrm{cm}^3\) of boiler volume, then we obtain:

\[ \text{production} - (\text{leakage} + \text{absorption}) = \text{rate of change of the density of thermal neutrons} \]

or

\[ Q-\left(-\frac{\lambda_t}{3}\Delta nv+nv\Sigma_a\right)=\frac{dn}{dt}, \tag{16,1} \]

i.e.

\[ \frac{\lambda_t}{3}\Delta nv-nv\Sigma_a+Q=\frac{dn}{dt}. \tag{16,2} \]

For simplification we shall first consider the kinetics of the boiler by the one-group method. In this case we put:

\[ Q=nv\Sigma_a k. \tag{16,3} \]

Introduce the quantity

\[ K^2=\frac{3\Sigma_a}{\lambda_t}(k-1)=\frac{k-1}{L^2}, \tag{16,4} \]

put the neutron lifetime equal to \(l=\lambda_a/v=1/v\Sigma_a\), and divide the equation by \(v\). Then equation (16,2) becomes

\[ \Delta n+K^2 n=\frac{l}{L^2}\frac{\partial n}{\partial t}. \tag{16,5} \]

In order to separate the variables, put

\[ n(x,y,z,t)=n_0(x,y,z)\,f(t). \tag{16,6} \]

After substituting this expression into equation (16.5) and dividing by \(n_0 f\), we obtain:

\[ -\frac{\Delta n_0}{n_0}+K^2=\frac{l}{L^2}\frac{1}{f(t)}\frac{df}{dt}. \tag{16.7} \]

Since the left-hand side is a function only of \(x, y\), and \(z\), while the right-hand side is a function only of \(t\), each part must be equal to a constant, which we shall denote by \(\frac{l}{L^3}\Lambda\).

Then

\[ \frac{1}{f}\frac{df}{dt}=\Lambda . \tag{16.8} \]

Integrating from \(t=0\) to \(t\), we obtain:

\[ f\sim e^{\Lambda t}, \tag{16.9} \]

where \(\Lambda\) is determined by the equation

\[ -\frac{\Delta n_0}{n_0}+K^2=\frac{l}{L^3}\Lambda , \tag{16.10} \]

obtained by substituting (16.9) into equation (16.7). If we denote

\[ -\frac{\Delta n_0}{n_0}=A, \tag{16.11} \]

then, using equations (15.1a) and putting \(n=n_0 e^{(\rho/l)t}\), we obtain:

\[ \Lambda=\frac{L^3}{l}(A+K^2)=\frac{1}{l}\,[k-(1-L^2A)]=\frac{\rho}{l}. \tag{16.12} \]

The interpretation of equation (16.12) is evident. The term \(-L^2A\), according to equation (12.1), is the ratio of leakage to absorption, i.e. the number of thermal neutrons leaving the boiler per one thermal neutron absorbed in the boiler.

Since, as a result of fission, \(k\) neutrons arise per absorbed neutron and one of these neutrons must be left to maintain the chain reaction, the quantity \(\rho=k-(1-L^2A)\) represents the total increase (or decrease) in the number of neutrons per absorbed neutron. Then, on the average, every \(l\) sec (where \(l\) is the mean lifetime of the neutrons) there will occur an increase (or decrease) in the number of neutrons in \(1\ \mathrm{cm}^3\) by the amount \(\rho n\). The rate of change of \(n\) will therefore be:

\[ \frac{\partial n}{\partial t}=\frac{\rho}{l}\,n. \]

This equation can be integrated immediately. Integration gives: \(n=n_0 e^{(\rho/l)t}\)—an expression coinciding with the one just obtained.

The quantity \(T=\frac{1}{\Lambda}=\frac{l}{\rho}\) is called the “period of the boiler.” It may be very long if \(\rho\) is very small compared with \(l\). In a boiler operating on thermal neutrons, \(l=1/v\Sigma_a\).

ELEMENTARY THEORY OF A REACTOR

may be of the order of \(10^{-3}\) sec, whereas in a reactor operating on fast neutrons this quantity is many orders of magnitude smaller. Thus it turns out that, if the preceding arguments described the whole picture completely, then it would be very difficult to control even reactors operating on thermal neutrons. For example, considering the extreme case, if the value of \(k\) for a critical reactor operating on thermal neutrons increases, upon rapid removal of a control rod, by \(0.1\), so that \(\rho\) becomes equal to \(0.1\), then the activity of the reactor will increase in one second by the factor

\[ e^{\rho/l} = e^{0.1/10^{-3}} = e^{100} = 10^{43}. \]

This catastrophic result should be regarded as a warning in the sense that unexpected large changes in the reactivity of a reactor must not be allowed.

In actual cases all reactors are equipped with special emergency devices that prevent such a catastrophe. It is evident that operation with reactors having a large excess of reactivity presents a great danger. One can say almost with certainty that ordinary control rods cannot serve for controlling a reactor with a period of \(1/100\) sec. However, in this way one can control a reactor with \(T = 1/10\) sec. This means that the control rods must be withdrawn at such a rate that the excess reactivity would never exceed \(1\%\). As we shall now see, some of these conclusions must be somewhat modified if one takes into account the role played in the reactor by delayed neutrons.

A favorable factor helping in reactor control has proved to be the presence, in addition to prompt neutrons, of a small fraction (\(0.76\%\)) of neutrons emitted in the fission process with a delay. Six groups of delayed neutrons have been found, with half-life periods from \(0.05\) to \(55\) sec. It has been shown that two of these periods correspond to definite nuclei—fission fragments—from which the neutrons are emitted after \(\beta\)-decay. It was found that a definite number \(\nu_{\mathrm{m}}\) (prompt neutrons) is emitted immediately after fission (in a time of the order of \(10^{-12}\) sec), whereas a small part \(\nu_{\mathrm{з}}\) (delayed neutrons) is emitted considerably later. Each group of delayed neutrons has its own characteristic period of the \(\beta\)-decay that precedes it. Thus the total number of neutrons emitted per fission is equal to the sum \(\nu_{\mathrm{m}}\) and \(\nu_{\mathrm{з}}\). The energy of delayed neutrons is somewhat less than the average energy of prompt neutrons. It is evident that the difference between delayed and prompt neutrons is important only in nonstationary states.

Data on delayed neutrons accompanying the fission of \(U^{235}\) are collected in Table III.

Table III

Half-life period $\tau_{1/2}$ (sec)	Decay constant $\lambda_i=\dfrac{0.693}{\tau_{1/2}}\ \text{sec}^{-1}$	Fraction (%) of all fission neutrons $\beta_i$	Neutron emitter
0.05	14	0.029	?
0.43	1.61	0.084	?
1.52	0.456	0.24	?
4.51	0.154	0.21	?
22.0	0.0315	0.17	Xe(137)
55.6	0.0125	0.026	Kr(87)
Total . . .	Total . . .	$\beta=0.76\%$	Total . . .

The effect of delayed neutrons consists in a considerable increase in the effective mean lifetime of neutrons, provided only that the change in reactivity is not so large (as in the catastrophic case described above) that the pile becomes supercritical on prompt neutrons alone. If $\beta$ is the fraction of delayed neutrons, and $l_3$ is the mean delay time, then the effective mean lifetime of a neutron is equal to: $l_{\mathrm{eff}}\sim l+\beta l_3$. Using the approximate values $\beta=1/100$ and $l_3\sim10$ sec, we find:

\[ l_{\mathrm{eff}}\sim 10^{-3}+\frac{1}{10}\sim\frac{1}{10}\ \text{sec}. \]

17. KINETICS OF A PILE. TRANSIENT REGIME

In the preceding section the transient regime of a pile was investigated by means of the one-group method, without direct consideration of the action of slowed-down neutrons. In this section we shall consider the kinetics of a “bare” homogeneous pile with delayed neutrons on the basis of the more exact method of Section 10.

We shall denote by $\lambda_1,\ \lambda_2,\ldots$ (or, in the general case, $\lambda_i$) the decay constants of the various groups, and the number of emitters of delayed neutrons arising per one fission event by $\nu\beta_i$. Then, if $\nu$ is the total number of neutrons emitted per one fission, then $\nu\sum_i \beta_i=\nu\beta=\nu_3$ represents the number of delayed neutrons per one fission, and $(1-\beta)\nu$—the number of prompt-

given neutrons per fission. Similarly, we may split the multiplication constant into two terms:

\[ k = k_{з}+k_{м}, \]

where

\[ k_{з}=\frac{\Sigma_f}{\Sigma_a}\beta \nu \tag{17,1} \]

and

\[ k_{м}=\frac{\Sigma_f}{\Sigma_a}(1-\beta)\nu . \]

The nonstationary equation for the flux of thermal neutrons in the boiler can be written, by analogy with equation (16,1), as:

\[ \frac{\lambda_t}{3}\Delta n v_\tau - n v_\tau \Sigma_a + \left(n v_\tau \Sigma_a k_{м}+\lambda c\right)e^{\tau_0 A} = \frac{dn}{dt}. \tag{17,2} \]

Here \(c\) is the number of fission fragments in \(1\ \mathrm{cm}^3\) emitting delayed neutrons. For simplicity we shall consider the case when there is only one kind of delayed-neutron emitter, and only later indicate how all six groups may be treated simultaneously. Let \(\lambda\) be the decay constant (the probability of emission in \(1\ \mathrm{sec}\)) for the delayed-neutron emitter. Then each second \(\lambda c\) delayed neutrons will be emitted. In equation (17,2) the source will be represented by the terms:

\[ \left(n v_\tau \Sigma_a k_{м}+\lambda c\right)e^{\tau_0 A}. \]

The concentration of such fission fragments \(c\) at a given moment depends on the number of fission events occurring during a small interval of time preceding that moment.

For the rate of change of \(c\) we have:

\[ \frac{dc}{dt}=n v \Sigma_a k_{з}-\lambda c, \tag{17,3} \]

where the first term gives the number of neutrons produced, and the second the number of decays in \(1\ \mathrm{cm}^3\) per \(1\ \mathrm{sec}\). Using the relations

\[ l_{м}=1/v\Sigma_a \quad \text{and} \quad k_{з}=\beta k, \]

we obtain:

\[ \frac{dc}{dt}=\frac{n\beta k}{l_{м}}-\lambda c. \tag{17,4} \]

In the case of equilibrium \(dc/dt=0\), and therefore

\[ c=\frac{n\beta k}{l_{м}\lambda} \tag{17,5} \]

is the concentration of delayed-neutron emitters after prolonged operation of the boiler, when their formation and decay are exactly balanced.

We must solve the system of differential equations (17,2) and (17,4). We assume the time dependences in the form:

\[ \left. \begin{aligned} n&=n_0 e^{\Delta t},\\ c&=c_0 e^{\Delta t}. \end{aligned} \right\} \tag{17,6} \]

Then \(\partial n/\partial t=\Lambda n\) and \(\partial c/\partial t=\Lambda c\). Substituting this into equations (17,2) and (17,4), we obtain:

\[ \lambda c=\frac{n\beta k}{l_m}\frac{\lambda}{\Lambda+\lambda}, \tag{17,7} \]

\[ \frac{\lambda_t}{3}Anv-nv\Sigma_a+ \left[ nv\Sigma_a(1-\beta)k+\frac{n\beta k}{l_m}\frac{\lambda}{\Lambda+\lambda} \right]e^{-\tau_0 A} = \Lambda n, \tag{17,8} \]

where \(\Delta nv/nv\) has been set equal to the constant \(A\), and expression (17,5) for \(\lambda c\) has been used. If we now multiply this equation by \(e^{-\tau_0 A}\), divide by \(nv\Sigma_a=n/l_m\), and use the fact that \(L^2=3\lambda_t/\Sigma_a\), we obtain:

\[ e^{-\tau_0 A}(L^2A-1)+k-\beta k+\beta k\frac{\lambda}{\Lambda+\lambda} = \Lambda l_m e^{-\tau_0 A}. \tag{17,9} \]

We define the reactivity, as before, by the relation

\[ \rho=k-(1-L^2A)e^{-\tau_0 A}. \tag{15,1} \]

Then equation (17,9) takes the form:

\[ \rho+\beta k\frac{\lambda}{\Lambda+\lambda} = \Lambda l_m e^{-\tau_0 A}. \tag{17,10} \]

This is a quadratic equation for \(\Lambda\); it may be written in the form:

\[ \Lambda= \frac{\rho/k}{ \displaystyle \frac{l_m}{k}e^{-\tau_0 A}+\frac{\beta}{\Lambda+\lambda} }. \tag{17,11} \]

In the case when several groups of neutrons are considered, characterized by constants \(\beta_i,\lambda_i\), equation (17,11) is replaced by the equation

\[ \Lambda= \frac{\rho/k}{ \displaystyle \frac{l_m}{k}e^{-\tau_0 A}+\sum_i\frac{\beta_i}{\Lambda+\lambda_i} }. \tag{17,11a} \]

When all six groups of delayed neutrons are taken into account, equation (17,11a) will be of order \(6+1\) with respect to \(\Lambda\), and therefore it has 7 roots: \(\Lambda_0\Lambda_1\ldots\Lambda_6\). The general solution can be expressed as a sum of exponentials of the form \(\sum a_i e^{\Lambda_i t}\).

Equation (17,11) can be written in the form:

\[ \Lambda^2 l+\Lambda(l\lambda+\beta-\rho/k)-\rho\lambda/k=0, \]

where \(l\) is abbreviated notation for \(\dfrac{l_m}{k}e^{-\tau_0 A}\).

Since the last term is negative, the two roots have opposite signs. Let \(\Lambda_1>0\) and \(\Lambda_2<0\). The solutions of the equations are:

\[ \Lambda_1 = -\frac{\beta+l\lambda-\rho/k}{2l} + \frac{1}{2l} \sqrt{ \left(\beta+l\lambda-\frac{\rho}{k}\right)^2+\frac{4l\rho\lambda}{k} }, \tag{17,12} \]

\[ \Lambda_2 = -\frac{\beta+l\lambda-\rho/k}{2l} - \frac{1}{2l} \sqrt{ \left(\beta+l\lambda-\frac{\rho}{k}\right)^2+\frac{4l\rho\lambda}{k} }. \tag{17,13} \]

The quantity \(\dfrac{\rho}{k} - \beta = \dfrac{1}{k}(\rho - k_3)\) is proportional to the excess of reactivity over the value still compatible with the possibility of controlling the boiler when delayed neutrons are taken into account, and is equal to the excess of reactivity, multiplied by \(1/k\), over the quantity at which the boiler becomes critical on prompt neutrons alone. Suppose that the reactivity is sufficiently small, so that the boiler is not critical on prompt neutrons alone. Then \(\beta > \rho/k\), and since \(l\lambda\) is small, formulas (17,12) and (17,13) may be written approximately in the form:

\[ \Lambda_1 \sim \frac{\rho/k}{l + \dfrac{\beta - \rho/k}{\lambda}}. \tag{17,14} \]

\[ \Lambda_2 \sim -\frac{\beta - \rho/k}{l}. \tag{17,15} \]

The general solution then consists of two terms:

\[ n = a_1 e^{\Lambda_1 t} + a_2 e^{\Lambda_2 t}. \tag{17,16} \]

The first represents a slow rise, and the second a transient effect, rapidly dying out with time.

As an example, let us consider the values:

\[ \frac{\rho}{k} = 0.005,\qquad \beta = 0.01,\qquad \lambda = \frac{1}{10}\ \mathrm{sec}^{-1},\qquad l = 10^{-3}\ \mathrm{sec}. \]

Then

\[ \Lambda_1 \sim 0.1\ \mathrm{sec}^{-1}, \]

\[ \Lambda_2 \sim -5\ \mathrm{sec}^{-1}. \]

The arbitrary constants are determined by the initial conditions at \(t = 0\), which we shall take as:

\[ n(0) = n_0 = \mathrm{const}; \qquad \frac{dc}{dt}(0) = 0 \tag{17,17} \]

and from equations (17,2), (17,5):

\[ \frac{dn}{dt}(0) = \frac{\rho/k}{l}\, n_0, \tag{17,18} \]

since the fact that some neutrons are delayed does not affect the initial slope of the curve. The constants \(a_1\) and \(a_2\), entering into the solution (17,16), are then determined with the aid of the equations:

\[ a_1 + a_2 = n_0 \tag{17,19} \]

and

\[ \Lambda_1 a_1 + \Lambda_2 a_2 = \frac{\rho n_0}{kl}. \tag{17,20} \]

the solutions of which are:

\[ a_1=\frac{\Lambda_2-\rho'kl}{\Lambda_2-\Lambda_1}\,n_0 \simeq \frac{\beta}{\beta-\rho/k}\,n_0 \tag{17,21} \]

and

\[ a_2=\frac{\rho'kl-\Lambda_1}{\Lambda_2-\Lambda_1}\,n_0 \simeq -\,\frac{\rho/k}{\beta-\rho/k}\,n_0 . \tag{17,22} \]

Collecting all the results obtained, we can write the general solution:

\[ n(t)=\left[\frac{\Lambda_2-\rho'kl}{\Lambda_2-\Lambda_1}e^{\Lambda_1 t} +\frac{\rho'kl-\Lambda_1}{\Lambda_2-\Lambda_1}e^{\Lambda_2 t}\right]n_0 . \]

Substituting numerical values, we obtain:

\[ a_1\sim 2n_0, \]

\[ a_2\sim -\,n_0 . \]

Thus, for this case the approximate solution will be:

\[ n \simeq n_0\left(2e^{0.1t}-e^{-5t}\right). \]

This solution is represented schematically in Fig. 21.

Fig. 21.

It is not difficult to see that the neutron density doubles in the first half-second. This means that at first the boiler “does not know” that it must “wait for” delayed neutrons. Subsequently the density increases more slowly, and the curve acquires the steepness characteristic of delayed neutrons. In the case when \(\rho/k>\beta\), the initial steep rise does not cease, since in this case the boiler is supercritical on prompt neutrons alone.