Calculation of the critical radius of a fissionable device
Comments on a note in the "Los Alamos Primer" (by Robert Serber and Richard Rhodes) and an article on neutron multiplication by Rudolf Peierls.
A Mathematica document by Sören Molander
Change notes:
2005-09-09: Fixed a typo in Eq 1''. Minor changes in text and elaborated an integral.
2005-09-10: Fixed a minus sign in the Eigenvalue calculation of Peierls formula. The plot of the Lf/f Eigen function is now correct. Also compared Peierls approximation with a more accurate one and made color coded plots.
2005-09-11: Checked the solution of Peierls forumla for β→α which checks much better with the Primer version. Rewrote some parts of the text
2005-09-14: Made numerical values modular for swift evaluation of the document for either U235 or PU239.
2005-09-17: Changed notation from total transport cross section to the more correct transport cross section.
2005-10-09: Fixed a factor of 1/2 missing in eq. 1''.
Introduction
One can say that modern nuclear physics started with Werner Heisenberg who applied quantum mechanics to the atomic nucleus 1932, the same year that James Chadwick discovered the neutron. In 1934 Enrico Fermi discovered unknown radioactive elements when bombarding a substance with neutrons and in 1938 Strassman and Hahn performed the first controlled fission experiment. In 1939 Lise Meitner and Otto Frisch used Niels Bohr's liquid drop model to conclude that Uranium was transmuted into Barium under the release of large amounts of energy. In 1939 an explosion of theoretical and experimental work was done on nuclear fission and the same year the notion of critical mass first saw the day, at least in published form. Leo Szilard and and Enrico Fermi were among the first to realize that nuclear fission could trigger a chain reaction, but they decided not to publish their findings to prevent Germany from using their results. It has been disputed if Werner Heisenberg tried delibrately to slow down the wartime German efforts towards a nuclear bomb, what is known is that the efforts were behind those of the allies.
When I did nuclear physics in undergraduate school I was interested in the physics of neutron multiplication and how it gave rise to the notion of "critical mass". Naturally I never expected to find a realistic calculation using detailed models (which is rightly classified information), but was more interested in a "back of the envelope" argument that highlighted the basic physics. However, I had never found any information on the internet or in books on this subject, so I decided to compile this document. It is an attempt to fill in some of the details of the critical mass calculations based on the article by Rudolf Peierls published 1939 [1] and the notes in the book The Los Alamos Primer[2]. The Primer was a set of 5 lectures given by Robert Serber in April 1943 for the physicists at Los Alamos on the basic physics and the related engineering problems of building a nuclear bomb. Much of the material had been worked out at a small group of people at Berkely headed by Robert Oppenheimer and including names like Hans Bethe, Edward Teller and Richard Tolman. The Primer was originally declassified in 1965 and has now been published in an extended book form by Robert Serber and Richard Rhodes, where additional material has been included.
Many of the physical properties of Uranium were not well known at the time when the document was written down and here I have tried to use the modern values of cross sections and densities (I am not sure about the value for the absolute geometric cross sections). The effect of an external shell (or "tamper", an outer shell which mirrors neutrons makes the effective critical mass even smaller) has been disregarded, readers are referred to the Primer and the seminal book by Richard Rhodes The Making of the Atom Bomb[3] for more details. One of the first notes on the critical mass and neutron multiplication was published by Francis Perrin [4]; I have not read this paper so I can't comment on it here.
First let us define some important physical properties of fissionable materials, the most central ones being the "cross section" (see the Endnotes in the Primer). The geometric cross section is the area taken up by a nucleus of radius R. Assuming R = 
cm one gets σ=π
= 3 
. Assume a thin foil of thickness d and area A with a number density of n = 
(
is Avogadro's number 6.022 
) neutrons per unit volume and M is the mass in atomic units. For U235 n ∼ 4.84 
(see below). To get some idea on how sparse the atoms are scattered in the foil, let's calculate the average distance between atoms. This is simply 1/
∼ 2.7 
or 27000 times the size of an atom and from the point of view of a small particle entering the foil at high speed, it sees mostly empty space. The total area filled with nuclei is thus nAdσ and a neutron passing through the foil has a chance of hitting a foil nucleus with probablity 
= ndσ. In a similar way, the probability of fission is given as a cross section 
and the fraction of collisions that give rise to fission is 
. The number of fissions in the foil is 
=nσd
=n
d. The time it takes a neutron to pass through the foil is d=vτ, where v is the velocity which for 1 Mev neutrons can be approximated by E=
and is about 1.4 
cm/s. This means that the rate of fissions is 
/τ and the inverse gives the time per fission, τ=
.
The cross sections here are given below for the relevant energies involved, around 1 Mev (the values here are given in the extended version of the Primer). Measurements made by several groups before 1939 (including Niels Bohr) showed that U235 had a significant fission cross section over a wide range of energies, which was important since only fast neutrons could make a significant contribution to a possible chain reaction in a bomb (unlike a nuclear fission reactor, where the opposite is true). Early experiments showed that the fission cross section of Uranium 238 drops sharply at energies < 1 Mev so natural Uranium has to be enriched to extract the isotope 235. Plutonium is a more efficient medium (higher fission cross section), it was discovered in 1942 by Glenn Seaborg and the isotope PU239 was used in the Trinity explosion 1945.
Technical notes
The Mathematica document can evaluated using values for U235 or for PU 239, see below. The first time the notebook is evaluated there will be some minor warning messages mostly related to the naming of variables (they will disappear the 2:nd time). The evaluation of the integral in eq. 12 takes a long time, please be patient. Many of the physical properties of PU239 are classified, so the numerical results for this case should be taken with a grain of salt, especially the total cross section. Mathematica consistently outputs numerical values to 15 decimals in the html-document. This is an artefact and the real accuracy is in the order of 1-2 decimal places.
 |
1.22*10^^-24 |
 |
4.*10^^-24 |
| Average number of neutrons per fission ν |
2.52 |
 |
18.9 |
| Velocity of 1 Mev Neutron [cm/s] |
1.38276*10^^9 |
| Time between fissions [s] |
1.22394*10^^-8 |
| Mean free path of U 235 Neutron [cm] |
5.16185 |
Critical mass
In kinetic gas theory the mean free path is defined as ℓ=
where n = number of particles per volume and 
is the transport cross section. The transport cross section is defined as the elastic scattering 
× 1-cos(θ) where θ is the scattering angle. If the neutrons were scattered homogenously in all directions 
= 
but in practice it is smaller and has to determined by laboratory measurements. An important relation is the one between the diffusion constant D and the mean free path, D=
where v is the average velocity and the constant D is a part in Fick's law of diffusion which states that the current J= -D∇N(r). As can be seen above, this gives an approximate mean free path of Uranium 235 of 5.1 cm. Since the mean free path turns out to be in the same order as the critical radius to be calculated, standard diffusion theory is not really applicable. This was the approach of Peierls in his paper of 1939 and he derives an integral equation for the number of neutrons N as a function of time and position. The reasoning is as follows: The number of scattered and secondary neutrons from an element at (x',y',z') in time dt is β vdt' N(x',y',z',t')dx'dy'dz' where v = neutron velocity and β is the number of neutrons emerging after each fission. When these neutrons have travelled a distance r they are reduced by a factor of 
(1/α is the mean free path) and they will fill a shell with radius r and thickness vdt .The density of neutrons after a time t' = t-
is therefore 
N(x',y',z',t') 
dx'dy'dz'. Integrating over the interior of the sphere leads to the integral equation
' (1)
Separating the variables N(x,y,z,t)=N(x,y,x)
and setting λ= 0 for criticality,
' (1')
Here α is defined as α=n(
+
+
) - scattering, diffusion and absorption cross sections - which is the inverse of the path length, i.e. ℓ=1/α. β is defined as the number of neutrons emerging from the fission process and β=n(
+ν
). For spontaneous fission - criticality - to occur β must exceed α. Peierls evaluates this integral in two extreme cases: 1-
≪ 1 and 1-
∼1 . The former case is the one leading to standard diffusion theory, and the other leads to an Eigenvalue problem, briefly discussed below. The integrand is assumed to be slowly varying and is expanded in a Taylor series, retaining the 2:nd order terms. Changing to spherical coordinates leads to 3 angular and 2 exponential integrals. The three angular integrals all yield the same value, and the x-integral is
and one of the exponential integrals is
Collecting terms gives
N(x,y,z)=β 
= β(N(x,y,z)
+ 
N(x,y,z)) (1'')
The limit is taken to infinity since the radius of the sphere is assumed to be much larger than the mean free path ℓ. Putting this together and yields a Helmholtz equation:
(2)
The solution is a spherical Bessel function of order 0 or
N(r) = 
(3)
with the wavenumber 
= 3 
(1-
). In the Primer, Serber simply states the following equation:
+ div J = 
(4)
In the extended annotated version of the Primer there is a note to this chapter that derives the equation in full. Assuming that diffusion theory is valid Fick's law J= -D∇N(r) gives the Laplace operator term. The source term on the r.h.s. is explained as follows: In each fission, ν number of neutrons produced, and 1 neutron is always absorbed (the one that triggers a fission in each step in the chain reaction) which together with the time between fissions give the number of neutrons produced per unit time,
(4')
Separation of variables N = N1(x,y,z)
gives
N1(r) + 
N1(r) = 0 (5)
with the same solution as (3). The boundary condition in the original version of the Primer is the simplest one possible is used: N1(R) = 0. This is most certainly not true but gives a first rough estimate of the critical radius. This boundary condition is satisfied when the wavenumber k=
or when the "effective" neutron number ν'=ν-1 - 
Dτ/
. If ν' is negative, too many neutrons escape for a chain reaction to be started and so the critical radius is
= π
Here is where kinetic gas theory enters. The time between fission was shown above to be τ=
=
where 
is the transport cross section, and since
D=
= 
=
, where k = 
(6)
Plugging in numerical values for U 235 gives a first estimate of the critical radius (~ 13 cm) and mass. Since the mass ∝ 
and 
∝ ℓ ∝
, the critical mass is proportional to ~ 
. If a mechanism could be designed to push the active core to a higher density, criticality could be reached. This was one of the turning points in the Manhattan project. Problems with the high spontaneous fission rates of PU239 made the "shooting" mechanism (used for the Uranium device in Hiroshima) impossible and instead led to the design an implosion mechanism used to "squeeze" a core of Pu 239 to achieve critical mass.
| Critical radius [cm] |
13.7506 |
| Critical Mass [kg] |
205.835 |
Critical mass with improved boundary condition
In order to yield a better estimate a more accurate boundary condition has to be introduced. This was done in Peierls' paper, and he cited a result from astrophysics:
0.71
+ αN(r)=0 (7)
The number 0.71 is an approximation valid for a critical radius >> mean free path. The number is an approximation from radiative transfer theory, where a scattering integral on the right hand side over the sphere is approximated with the average value. On the left hand side, the cosine obliquity factor μ ≡cos(θ) is replaced with an "average" cosine factor
. On such choice is a root mean square average
= 
= 
. Assuming that N depends linearly on μ N ≈Cμ, the result is 
=1/
≈0.71
In the Endnote of the Primer, Serber gives a more intuitive explanation. Imagining a thin slice in the material of width 2ℓ (2 times the mean free path), the flow of neutrons going out to the left is 1/2 N1 v and to the right 1/2 N2 v. The net current is J = 1/2(N1-N2)v. The difference N1-N2 = 
2ℓ, and this gives J = -ℓv
. Generalizing this to 3 dimensions readily gives J = -
. At the boundary the neutrons have "nowhere" to go except out, so this gives them only 1 degree of freedom. This gives the following system of equations:
N v = -
(B.C.) (8)
N(r) = 
(Eq 3)
Putting the solution of the Helmholtz equation (3) in the boundary condition (8) gives
= -
(
)
x cot(x) = 1-
(9)
with b=
and x=kr.
This is identical to equation (10) given in Peierls except that α=kℓ and the factor 2/3 is replaced with 0.71. Solving (7) numerically for U235 (see below) gives x = 2.26 and
= 
= 2.2698 ℓ 
(10)
Comparing this to Eq. 6 gives 2.26/π ~ 0.72 smaller radius than the original estimate. Note that the mass is ~ 80 kg instead of 210 kg and the effect of tamper around the device will further reduce the critical radius. Serber states that equations above are still empirical and that the full problem with an exact boundary condition was solved exactly by two people - Eldred Nelson and Stan Frankel - in the Berekely theory group. The Primer gives no clue to this, but most probably the solution is based on an integral equation like in Peierls' paper. He also states that this solution only differed slightly from the one from diffusion theory with the improved boundary condition.
The root above was found using a numerical procedure, which was cumbersome in the days when numerical calculations required tables and had to be checked and double checked for errors. In the paper by Peierls a simple approximative formula for the root is derived, using a Taylor expansion of x cot(x) around Pi. Using Mathematica it's easy to derive an approximation formula for Eq. 9. Expanding around Pi to first order gives (remember that b = 2/3 k ℓ)
= 1-
Rearranging terms gives
- x (
+π) + 
(
-π) = 0, A = 
-
The solution is given below (the first root is physically relevant). Compared to the numerical solution 2.26, the error is ~ 3 %. The numerical value for b ~0.78 (U235), and the approximation gets better for b < 1 as can be seen on the graph below, or when
b = 
k ℓ = 
< 1.
Peierls' formula for critical mass
In Peirls paper a simple formula is derived for both the low and high beta regimes and he produces a graph of the critical radius valid for both regimes of approximation:
= 
x + 0.71 
= 0.552 x + 0.216
, for x = 
≪1 (from Taylor expansion of x cot(x))
= 0.78 - 1.02 (1-x) , for x ~ 1 (from Eigenvalue calculation)
(remember that β=n(
+ν
) ). The first part of the formula can be derived as follows: From eq (7) and the solution of the diffusion equation (3)
rk cot(rk) = 1-
where b = 0.71
Taylor expansion around Pi yields:
Now discard the quadratic term to get a managable polynomial equation of second order and set x → r k. Together with the RHS of the equation above this gives:
Solving for the reciprocal of r, and using the wavenumber from Eq. (2) k = 
α
= 
and Taylor expanding around x=0 yields:
Using α = 
~ β for x << 1, yields the desired forumla. Note that the approximation is valid only for x ≪ 1 typically useful in the range x = [0, 0.2]. The blue curve shows the deviation from square root solution above, and it follows closely up to values ~ 0.4.
The second part of Peirls formula (x ~1) is derived as follows. Starting from Eq. 1', in this approximation, the sphere is small compared to 1/α and the exponent is thus negliable, leading to
(11)
Assuming that there is no anglar dependence in N, using the cosine theorem for 
→ 
+
- 2 r r' cos(θ) yields
N(r) = 
=
. Putting x = r/R and f = nr
f(x) = 
This is an Eigenvalue problem, λ f = L f, with λ = 
. The kernel is not bounded, but it's square is which means that one can get the highest Eigenvalue from the Rayleigh-Ritz quotient
< 
< 
.
The idea now is to iterate the integral using eigenfunctions to get an estimate of λ. The simplest choice here is f = x, which leads to
(1-
) 
The value of λ = 
for x=0, and x=1 is in the range of [1,2]. An obvious improvement is
f = x - b
, which gives (see below)
Lf = 
(1-
) l
- b { 
(1-
) l
+ 
+ 
x} (12)
[Incidentally, there is an error in Peierls' paper: The polynomial 
(1- 
) in the second term should be 
(1-
). I verified the calculations in the paper with mathematica using the correct formula, and the numerical results are identical to those in the paper so it must have been a typesetting problem]. The constant b is determined such that Lf/f assumes the same value for x=0 and x=1. Using the fact that 
(1-x) log(1-x) = 0 gives the value for b = 0.633975 (see below). The first plot shows the Eigenfunction and the second plot shows the ratio 
for this function in the same interval. Peierls probably tabulated this function around x ~ 0 to get a robust estimate of λ. It is interesting to note that results by Mathematica tally almost exactly with the numerical values in the paper, which gives some credit to enormous work that scientists in Peierls days had to put into evaluating integrals and functions.
The eigenvalue at x = 0 is then λ=
= 1.57735. This approximation was done assuming α = 0, and this is now improved using a first order Taylor expansion of Eq. 1':
lf = Lf + M f = 
sin(θ) dθ dφ dr = Lf - 
dθ dr . The θ integral turns out to be
-2α (|r + r'| - | r - r'|) 
and so
M f = -αr
Now we turn to first order perturbation theory which states that a small perturbation of the value is given by
δ l = - 
Using the eigenfunction in Eq 12 gives
This means that
= 
-
= 0.788 - 0.398 α r . Now Peierls uses a trick: For small α use the fact that x = 
∼1 - 
and thus β = 
. Now use the value for r when α=0 on the rhs: r ∼ 
= 
and so
= 0.788 - 
= 0.788 - 1.01(1-x)
Graphing this gives the following plot (fig. 1 in the original paper). Since none of the formulae are valid in the intermediate range, some interpolation must done. Thus Peierls formula must be seen as a back-of-the-envelope estimate and more accurate values must be calculated numerically as in eq. 9
Note that the definition of the cross sections in Peierls and the Primer differ, and α contains the scattering, absorption and the disintegration cross sections. In the Primer the transport cross section is the sum of the fission c.s and an angular integral over the scattering cross sections. Thus the values given above cannot be directly used in Peierls' formula. Nevertheless, if one tries to identify alfa with the reciprocal of the average mean free path, this gives:
Identifying α and β gives
α → n 
and
or
β → 
Using the values above for the cross sections for U 235 one sees that Peierls original formula underestimates the critical radius as compared to Eq. 10, which is not surprising given the number of approximations involved. Also more accurate values for the different cross sections are probably needed to make an honest comparison beteeen the formulae
 |
0.680882 |
| Peierls critical radius[cm] |
5.81703 |
| Peierls critical Mass [kg] |
15.5831 |
| Modern value of radius (U235) |
8.91003 |
| Modern value of mass (U235) |
56. |
| Modern value of radius (PU239) |
5.09629 |
| Modern value of mass (PU239) |
11. |
When Otto Frisch and Rudolf Peierls first used the formula the results were astounding: Instead of a few tons of uranium, less than a few tens of kilos was needed. In the famous Frisch-Peierls memorandum from March 1940 the figure 2.1 cm is mentioned, which was one of the reasons the memorandum become a serious concern in the US and the UK during the war. Apart from the fact that Peierls' formula underestimates the critical mass, it turns out they used the wrong numerical values for many of variables involved, including the density, the fission cross section and the number of neutrons per fission and in reality the mass was larger. A modern value for the critical mass for a bare sphere of U235 is 56 kg which corresponds to a radius of 8.9 cm. In practice the "device" (=bomb) the critical mass was determined by measuring neutron multiplication from a neutron source surrounded by a sphere of the radioactive material. By making the sphere successively larger, the multiplication could be extrapolated and the critical mass measured (Primer, page 33).
References
[1] Rudolf Peierls, 1939. "Critical conditions in neutron multiplication". Proc. Camb. Phil. Soc. 35:610.
[2] Robert Serber, Edited with an Introduction by Richard Rhodes (1992) "The Los Alamos Primer". University of California Press, ISBN-0-520-07576-5.
[3] Richard Rhodes, 1986. "The Making of the Atom Bomb".Touchstone Book, Simon & Schuster New York, ISBN 0-648-81378-5.
[4] Francis Perrin, 1939."Calul relatif aux conditions évenetuelles de transmutation en châine de l'uranium. Comptes Rendus, 208:1394.
