We know from the laws of physics that the radioactive decay of an isotope is a random process. The probability of decay per unit time is constant, independent of time, and is given by the decay rate λ. This implies that the probability for an isotope to decay in a sufficiently small time period dt is given by λdt (when λdt << 1). There are cases for which the decay rate is not constant but they will not be elaborated here.
Suppose that we have N0 atoms of an isotope with decay rate λ at the beginning of time. The mathematical formula that describes the number of atoms present at a later time is given by
However, we will try to arrive at the same exponential decay form by using a Monte Carlo technique.
Assuming that we are looking at the number of atoms in sufficiently small time steps so that λdt << 1. The probability that one of these atoms decay within each time step is λdt so we can test if an atom have decay or not by comparing a random number generated from a uniform distribution [0,1) by this probability. This comparison is performed for all atoms, in each time step. Let us write this in pseudo code:
number of atoms = N
for each time step dt:
for each atom:
R = random number between 0 and 1
if( R < lambda * dt ): N = N-1
Save N for this time step.
So, if we generate a random number that is smaller than λdt, we say that it decays. On the average, after many samples of this kind, the fraction of decayed atoms will converge to the true fraction that would have decayed if we infinitely many atoms.
Let us look at 137Cs. It has a half life of 11018.3 ± 9.5 days, according to NIST. Its decay constant is ln(2)·T-1½ ≈ 22.977·10-3 per year. Suppose we start with 1000 atoms of 137Cs. Running the code above, coded in Javascript with a time step of one year, we get the following number of remaining atoms for the first few years:
We can plot the number of remaining atoms as a function of time:
You can compare the number of remaining atoms with the theoretical exponential decay function. The number of atoms calculated with the Monte Carlo technique is in better agreement with reality since nature is probabilistic as demonstrated above. At one half life, at about 30 years for 137Cs, we have about half of the atoms remaining.