Chain Reaction
Philipp K. Janert
May/June 2006
Abstract:
We look at the decay of a system of unstable particles.
Depending on the system size and the size of the particles
a chain reaction may, or may not, occur. We derive a scaling
relation which simplifies the problem significantly.
Consider a crystal, i.e. a regular arrangement of atoms on a lattice.
The nucleus of one of the atoms decays spontaneously, emitting two
neutrons. One of the neutrons passes so close to one of the remaining
atoms that it is captured by it. This nucleus, in turn, decays,
emitting two neutrons. If the entire crystal is large enough, there
is a good chance that one or both of these neutrons will get caught,
leading to a chain reaction of decaying atoms. On the other hand, if
the crystal is too small, the neutrons will leave the crystal before
striking any other atom.
For simplicity's sake, let's assume a square lattice in two dimensions,
where the atoms are located at the integer positions: ( m, n ) with
m, n both integer.
Now the atom located at ( m1, n1 ) emits a neutron, which travels
in a straight line until it is either captured by another atom or
leaves the crystal. The direction of the neutron can be described
through a vector
of unit length, where
is the angle with the horizontal.
Figure 1:
An atom decays sending out a neutron (red). The atom passes
the second atom in a distance beta (green).
 |
How close does the neutron traveling in this direction come to an
arbitrary other atom located at ( m2, n2 )? In vector notation,
we can write (also cf. Figure 1):
where we have made use of the fact that in two dimensions the
vector ( -
y,
x ) is orthogonal to (
x,
y ) (in three or more
dimensions, one has to work a little harder). The coefficients

and

are real numbers which control the lengths
of the vectors and are as yet undetermined.
We can eliminate
by realizing that the above vector
equation is equivalent to the following two scalar equations:
Now, multiply the first equation by
and the second
equation by
and add (making use of the well-known
trigonometric identity
) to find:
A neutron starting at (
m1,
n1 ) and traveling in direction

will pass a second atom at (
m2,
n2 ) in a distance
which is given by the absolute value of

according to the
previous formula.
The neutron will be captured by the first atom it passes at a distance
shorter than than the critical scattering length
. In other words,
the neutron will be captured by the first atom it passes for which
.
In principle then, the program is very simple: we start out with a
single pair of neutrons, choose random directions for both, and then
compare each neutron to all the other atoms in the crystal (which
have not yet decayed) to see whether the above relation holds. If
yes, the neutron is captured and the capturing atom will decay at
the next time step. (We will discuss a better implementation below.)
We let the reaction continue until there are no more free neutrons
around and every atom which had captured a neutron has decayed. Once
the reaction has come to an end, we calculate the fraction f of
atoms which have survived.
The results are shown in Figure 2.
Figure 2:
Fraction of surviving atoms f as function of the system
size in N, where N is the number of atoms along each edge of
the crystal. The scattering length
of the atoms in the
crystal varies from
to
.
 |
The graph shown demonstrates what we should have suspected earlier:
that there are really two parameters in this system. One is
the system size (measured in the number N of atoms along each edge
of the crystal) and the other is the scattering length
of
each atom, so that the fraction f of surviving atoms apparently
is a function of both parameters:
For sufficiently small N and/or sufficiently small
,
the reaction dies out quickly and almost all of the atoms survive.
But for larger crystals and/or larger scattering lengths, the
reaction continues until most of the atoms will decay. The
transition between these two regimes is rather sharp.
The interesting fact about the graph is that it suggests that
we can change either the crystal size N or the
scattering length
to achieve the same effect. The
question is whether we can make this intuition more precise.
Here is a simple, intuitive argument. Consider an atom at the
center of the crystal. Part of the
(or
)
``view'' from this atom will be obscured by the other atoms in
the crystal. It seems reasonable to assume that the fraction
f of surviving atoms depends on the fraction
of the
``horizon'' which is occupied by other atoms in this way. This
fraction clearly depends on both the size of the crystal
(measured in N) and the size of the atoms (measured in
).
To estimate the fraction
of the total ``view'' or
``horizon'' occupied by other atoms, let's first ask ourselves
how large another atom at a distance R will appear. A single
atom will obscure an angle
, such that
. For
(i.e. when the
atom is relatively far away, which is always the case, since we
assume that
is much smaller than the spacing of the
atoms),
is small and we can approximate:
Now, the fraction

of the total ``horizon'' occupied
by atoms at a distance
R is:
where
k is the number of atoms at a distance
R from the
center of the crystal. The number
k is proportional to the
length of the circumference of a circle of radius
R
(

), since the atoms are spaced at a uniform
distance from one another. Therefore:
where
C is a constant of order unity.
(The symbol

indicates that the left-hand side and the right-hand
side are proportional to one another, i.e. they are equal to within a
constant factor.) It may come as surprise that

is independent
of
R, but this can be understood easily: each atom appears smaller by
a factor 1/
R, but there are
R-times as many on the circumference, so
the effects cancel.
This expression tells us the fraction of the horizon occupied
by atoms at a distance R from the center. To find the
total fraction of the horizon occupied, we have to sum
over all atoms in the crystal, that is we have to integrate
over R:
The total size of the crystal
Rmax is proportional
to the length of each ``edge'' of the crystal:

.
With this we can find the final result:
This expression expresses the total fraction of the ``horizon'' covered
by other atoms in a crystal of size
N, with atoms of diameter

.
If we assume that the fraction f of surviving atoms depends only on
this quantity, then we would expect that f will stay the same for
different N, provided we arrange
such that
stays
constant:
If this hypothesis was true, we would expect all values for

to fall on a
single curve when plotted
as function of

. The result is shown in Figure 3,
and the collapse of the data onto a single curve is immediately
evident. The scaling relation derived above clearly holds!
This is a very interesting result, because we have been able to
simplify the problem significantly, even without having a detailed
understanding of the dependence of f on N and
. What
we have been able to achieve is to reduce a function of two variables
(i.e.
) to a function of only one variable,
by showing that the apparently independent arguments N and
in reality only occur in the combination
in f.
A function of a single variable is always easier to understand than
one of several -- the single curve in Figure 3 is easier to understand
than the family of curves in Figure 2!
Figure 3:
Same data as in Fig. 2, but plotted as function of the
fraction of the horizon covered
. Note the
almost perfect collapse of the data.
 |
Note that we have been able to make this simplification without
knowing anything about the form of f -- we have no
explicit formula for f as function of
. This is typical
for situations allowing the application of scaling relations: a
scaling argument based on physical intuition allows to simplify
the problem significantly, but a closed form for the result is
still missing.
Based on the discussion so far, the implementation is rather
straightforward. There is one optimization, however, which
should be mentioned.
The naive implementation, which compares each neutron to all
atoms in the crystal to find the closest distance is horribly
inefficient:
(since we
have an
process for all
decaying atoms in the crystal). Since
the neutrons move in straight lines, and the crystal is a regular
arrangement of points, we should be able to find an algorithm which
only looks at
atoms to
determine neutron capture.
This can be achieved in several ways, all of which require a little
bookkeeping. The algorithm used in the code sample first determines
whether the
- or
-component of the neutron's direction is larger.
Let us assume that it is the
-component which
is larger. Then we can describe the motion of the neutron as consisting
of steps of unit length in
-direction, with
correspondingly smaller steps in
-direction.
At each step, the neutron will be positioned between two atoms: one
exactly ``above'' and one exactly ``below''. Only these two atoms are
possible capturing partners, all other atoms in the crystal being further
away. If we continue stepping in
-direction
until the neutron leaves the crystal, we evaluate all atoms which
have a chance of capturing the neutron.
If the
-component of the neutron's direction is
the larger one, the same discussion applies, but with the
- and
-components
interchanged.
Simple as this system is, it has a surprisingly rich phenomenology
and the discussion in this article has just started to scratch the
surface. Here are some additional questions one may ask:
- We have only considered the fraction of surviving atoms, but other
quantities can be studied, for instance the number
of time steps it takes before the reaction peters out. This time
also depends on N and
-- can a similar scaling relation
be found?
- Here, we have only concentrated on the final state when the reaction
has run its course. We can also study the dynamic behavior of the
system. What is the number of free neutrons as a function of time?
Clearly, this number first grows, then falls to zero. What else can
we know about it?
- We studied the system's dependence on two variables: the
system size N and the size of the atoms
. In reality,
there is (at least) a third parameter, namely the multiplication
factor of neutrons. We have kept this fixed, so that each decaying
atom emits exactly two neutrons. But we can make this number a
random variable as well, with an arbitrary distribution and only
fix its average. In this case, the multiplication factor does not
even have to be an integer. How does f scale with this multiplication
factor?
- What is the structure of the remains of the crystal? Are all the
decayed atoms clustered towards the center or are they spread
throughout the crystal? To study this questions, I suggest fixing
at roughly 1.5, and then selecting N to be as large
as possible (
) while keeping
small.
- Finally, we may ask whether it is possible to maintain a steady-state
chain reaction by replenishing decayed atoms at the appropriate rate.
We can do this (unphysically) by randomly ``re-growing'' decayed
atoms, or (more realistically, but requiring more programming skill)
by moving two independent lattices towards one another at the
required rate. Alternatively, we can also consider inserting an
inert material to absorb excess neutrons to throttle the chain
reaction at the desired rate. What is necessary to make this
process stable?
The idea for this problem was taken from A. K. Dewdney's
Computer Recreations column in Scientific American (March 88).
This column, together with many others has been reprinted in
book form: The Magic Machine by A. K. Dewdney
(Freeman and Company, 1990).
Scaling relations as described in the text have become very
popular in physics and applied mathematics in the last twenty
years. Two books devoted exclusively to them are Scaling
by Grigory Isaakovich Barenblatt (Cambridge University Press, 2003)
and Scaling, Self-similarity, and Intermediate Asymptotics:
Dimensional Analysis and Intermediate Asymptotics by the same
author (Cambridge University Press, 1996). The former book is
more elementary, the latter more advanced.
- - chainreaction.perl
- A Perl script to study the chain reactions in a 2-dimensional crystal.
©
2006 by Philipp K. Janert.
All rights reserved.
www.toyproblems.org