Nuclei consist of neutrons and protons bound together. But protons are positively charged and hence repel one another. Thus large nuclei consist mostly of neutrons and some, having large numbers of protons, are quite unstable.

Some of them, from time to time, give off some energy and an electron or helium nucleus, and change themselves into some other nucleus. Each of these changes is described by what is called a "half life"; which is the time it takes for half of the nuclei of the given kind to decay.

We model this by the assertion that in a population of nuclei of a given radioactive isotope, in any small finite time interval, , a certain proportion, , of them will decay.

If the population of these nuclei is represented as the function p, then p(t) will obey

which corresponds to the differential equation

This is a differential equation whose solution we know, because it states that the function p is a constant multiple of itself.

This is precisely the property of p(0) * exp(-ct), which is its solution.

In the case of an isotope that only undergoes a single decay, this is the whole story. The only question to raise is, what is the relation between the constant c here and the "half-life" of this decay?

Recall that the half-life of the process is the length of time needed for the population to diminish to half its original size. The half-life, T, of this process then obeys

Upon dividing by p(0) and taking logarithms of both sides we obtain

A more interesting situation occurs when the decay product in the decay of the isotope you start with is itself radioactive, and decays itself with its own half life.

We can model this situation by defining p0(t) as the population of the original isotope, p1(t) the population of its decay product, and p2(t) as the population of the decay product of the decay product.

Then we have the same model as before for p0(t):

In the previous case, the population of the decay product was p(0) minus p(t). Now we must write a differential equation for it:

Here the first term represents the increase in population of the first decay product coming from decay of the original isotope, while the second term represents the effect of its own decay. (The decay constants associated with these two decays are here c1 and c2).

We do not need an additional equation for p2, though we could write one with the same approach,, because the sum, p0(t) + p1(t) + p2(t), will be p0(0), and we can deduce p2(t) once we have calculated the others, assuming that we started with a pure sample of the starting isotope, at least approximately.

Thus, in this case we have two dependent variables, p1 and p2, and two differential equations to solve.

We are studying modeling in this section so will not solve these equations in detail. However they present little difficulty. The equation for p0(t) can be solved exactly as the one for p(t) in the first case, with the same solution. This solution can be substituted for p0 in the second equation, and it is straightforward to solve the resulting equation. We will solve this equation when we study such things later on.

Please notice that these models only make sense when the initial population is quite large, since the only possible changes in population here are integers. In other words, a given nucleus either decays or doesn't, so at any time the number in each population is an integer. Thus, if we make dp smaller than 1, we leave the realm in which the model makes sense, and enter an artificial mathematical domain. Thus, keeping dp finite, as we do in numerical calculations, rather than letting it go to 0, as we do in formal differentiation is more compatible with the model.

Fortunately, we get essentially the same answer using finite or "infinitesimal" differences, in problems of this kind.

Exercise 10.1 Suppose isotope 0 decays to isotope 1 which decays to isotope 2 which decays to isotope 3. (Radium decays to lead by a still longer sequence). Set up a model of variables, decay constants, and differential equations to describe this system.