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.

An Unstable System

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 ( m₁, n₁ ) 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 $( \cos \phi, \sin \phi )$ of unit length, where $\phi$ 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 ( m₂, n₂ )? In vector notation, we can write (also cf. Figure 1):

$\begin{displaymath} \begin{pmatrix}m_2 \ n_2 \end{pmatrix} - \begin{pmatrix}m_... ... + \beta \begin{pmatrix}-\sin \phi \ \cos \phi \end{pmatrix}, \end{displaymath}$

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 $\alpha$ and $\beta$ are real numbers which control the lengths of the vectors and are as yet undetermined.

We can eliminate $\alpha$ by realizing that the above vector equation is equivalent to the following two scalar equations:

$\begin{align*} \alpha \cos \phi - \beta \sin \phi & = m_2 - m_1 \\ \alpha \sin \phi + \beta \cos \phi & = n_2 - n_1 \end{align*}$

Now, multiply the first equation by $- \sin \phi$ and the second equation by $\cos \phi$ and add (making use of the well-known trigonometric identity $\sin^2 \phi + \cos^2 \phi = 1$ ) to find:

$\begin{displaymath} \beta = - \sin \phi ( m_2 - m_1 ) + \cos \phi ( n_2 - n_1 ). \end{displaymath}$

A neutron starting at ( m₁, n₁ ) and traveling in direction $\phi$ will pass a second atom at ( m₂, n₂ ) in a distance which is given by the absolute value of $\beta$ 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 $\sigma$ . In other words, the neutron will be captured by the first atom it passes for which $\vert\beta\vert \leq \sigma$ .

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 $\sigma$ of the atoms in the crystal varies from $\sigma =0.01$ to $\sigma =0.15$ .

A Scaling Argument

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 $\sigma$ of each atom, so that the fraction f of surviving atoms apparently is a function of both parameters:

$\begin{displaymath} f = f( N, \sigma ) \end{displaymath}$

For sufficiently small N and/or sufficiently small $\sigma$ , 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 $\sigma$ 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 $360^\circ$ (or $2 \pi$ ) ``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 $\gamma$ 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 $\sigma$ ).

To estimate the fraction $\gamma$ 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 $\psi$ , such that $\tan \psi = 2 \sigma / R$ . For $\sigma \ll R$ (i.e. when the atom is relatively far away, which is always the case, since we assume that $\sigma$ is much smaller than the spacing of the atoms), $\psi$ is small and we can approximate:

$\begin{displaymath} \psi \approx \tan \psi = \frac{2 \sigma}{R}. \end{displaymath}$

Now, the fraction $\gamma$ of the total ``horizon'' occupied by atoms at a distance R is:

$\begin{displaymath} \gamma( R ) = k \frac{\psi}{2 \pi} = k \frac{2 \sigma / R}{2 \pi}, \end{displaymath}$

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 ( $k \sim 2 \pi R$ ), since the atoms are spaced at a uniform distance from one another. Therefore:

$\begin{displaymath} \gamma( R ) = C \sigma \sim \sigma, \end{displaymath}$

where C is a constant of order unity. (The symbol $\sim$ 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 $\gamma(R)$ 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:

$\begin{displaymath} \Biggl. \gamma = \int_0^{R_{\mathrm{max}}} \gamma( R ) dR \sim \int_0^{R_{\mathrm{max}}} \sigma dR \Biggr. \end{displaymath}$

The total size of the crystal R_max is proportional to the length of each ``edge'' of the crystal: $R_{\mathrm{max}} \sim N$ .

With this we can find the final result:

$\begin{displaymath} \gamma \sim N \sigma. \end{displaymath}$

This expression expresses the total fraction of the ``horizon'' covered by other atoms in a crystal of size N, with atoms of diameter $\sigma$ .

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 $\sigma$ such that $\gamma$ stays constant:

$\begin{displaymath} f( N_1, \sigma_1 ) = f( N_2, \sigma_2 ) \quad \textrm{if} \quad \gamma_1 = N_1 \sigma_1 = N_2 \sigma_2 = \gamma_2. \end{displaymath}$

If this hypothesis was true, we would expect all values for $f(N, \sigma)$ to fall on a single curve when plotted as function of $\gamma = N \sigma$ . 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 $\sigma$ . What we have been able to achieve is to reduce a function of two variables (i.e. $f(N, \sigma)$ ) to a function of only one variable, by showing that the apparently independent arguments N and $\sigma$ in reality only occur in the combination $N \sigma$ 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 $\gamma = N \sigma$ . 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 $\gamma$ . 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.

Implementation Notes

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: $\mathcal{O}(N^4)$ (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 $\mathcal{O}(N)$ 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.

Questions and Alternatives

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 $\sigma$ -- 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 $\sigma$ . 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 $\gamma$ at roughly 1.5, and then selecting N to be as large as possible ( $N \approx 100$ ) while keeping $\sigma$ 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?

ToyProblems

About