Probability and differential equations

Some atomic nuclei are energetically unstable and can undergo spontaneous radioactive decay to form a more stable nucleus. The process of radioactive decay is purely random and, in a given fixed time period, a nucleus has a fixed probability of decaying. The process is also memoryless. This means that it doesn't matter whether the nucleus has gone one time period without decay or one hundred time periods, the probability of decay in the next time period remains the same.

A time period of particular use is the time for which the probability of decay is 0.5. Such a probability means that over the time period, half of the nuclei will decay and half will not. This time period is called the half-life of the nucleus. If I start will 100 g of a radioactive element, then after one half-life there will be 50 g remaining. Remember the decay is memoryless so that waiting another half-life will reduce the 50 g to 25 g as there is a fixed probability of decay of 0.5 for each nucleus and, as there are half as many nuclei to start with in this time period, the total number that will decay is also half as many.

The number of nuclei remaining afer $$N$$ half-lives, $$A_N$$ is given mathematically by $$A_N = \frac{A_0}{2^N},$$ where $$A_0$$ is the initial number. If you are only considering the nuclide that is decaying the above equation also applies to the mass. Care must be taken with mass however! Consider 100 g of polonium-210, $$^{210}Po$$. After one half-life (138.38 days) there will be 50 g left. However, there will not be 50 g of its stable decay daughter lead-206, $$^{206}Pb$$, as each nucleus of lead has approximately 2% less mass than each nucleus of polonium, leading to there being about 49 g of lead. The missing mass is largely accounted for by the alpha particles (a helium nucleus consisting of two protons and two neutrons) emitted during the decay process.

The interactive element below simulates the radioactive decay of red nuclei to yellow nuclei. Three different half-life values can be selected and line markers are applied at every half-life interval. Notice how in each interval, the number of red nuclei halves, following the equation given above.

The plot in the interactive element above is drawn using time intervals that are much smaller than the half-life so a smooth curve is drawn. This has an impact on the probability of decay value that should be used. To see why, consider the following equation for how the number of nuclei changes over a time period: $$A_{t} = A_0 - pA_0 = A_0 \left( 1-P \; \right),$$ where $$P \;$$ is the probability of decay over the time period. If the time period is one half-life then, by definition, $$P = 0.5$$. Now consider taking two time steps, each of half a half-life. A first guess is that as the time-step is half as long, the probability of decay should be half as large, so that after the first half time-step the number of nuclei is $$A_{\frac{t}{2}} = A_0 \left( 1-\frac{P}{2} \; \right),$$ and after the second half time-step the final value is $$A_{t} = A_{\frac{t}{2}} \left( 1-\frac{P}{2} \; \right) = A_0 \left( 1-\frac{P}{2} \; \right) ^ 2.$$ Using $$P = 0.5$$ in the above will give the wrong answer - it will not yield a reduction by half over the half-life. To achieve this a value of $$P = 0.586$$ (3 d.p.) is required. This exercise can be extended to ever smaller sub-divisions of the half-life. When this is done the value of $$P \;$$ converges to the natural logarithm of 2 (about 0.693). The above showed that when the half-life interval is divided into two sub-intervals, the number of nuclei is prediced by $$A_{t} = A_0 \left( 1-\frac{P}{2} \; \right) ^ 2,$$ and as, by definition, $$A_{t} = \frac{1}{2}A_0$$, a value of $$P \approx 0.586$$ is required to satisfies this. The above can be extended, for example if the half-life interval is divided into four sub-intervals, the number of nuclei is prediced by $$A_{t} = A_0 \left( 1-\frac{P}{4} \; \right) ^ 4,$$ and generally when the interval is sub-divided into $$N$$ substeps $$A_{t} = A_0 \left( 1-\frac{P}{N} \; \right) ^ N,$$ and a value of $$P \;$$ is required such that $$\left( 1-\frac{P}{N} \; \right) ^ N = \frac{1}{2},$$ as $$N \rightarrow \infty$$. To solve this, write $$y = log_e \left( 1-\frac{P}{N} \; \right) ^ N = N log_e \left( 1-\frac{P}{N} \; \right)=\frac{log_e \left( 1-\frac{P}{N} \; \right)}{\frac{1}{N}}.$$ Allowing $$N \to \infty$$ in the above yields $$\frac{0}{0}$$. To avoid this L'Hopital's rule is used: $$\lim_{x\to c} \frac{f(x)}{g(x)}=\lim_{x\to c} \frac{f^\prime (x)}{g^\prime (x)}$$ Writing $$log_e \left( 1-\frac{P}{N} \; \right) = log_e \left( \frac{N-P}{N} \; \right) = log_e \left( N-P\right) - log_e \left( N\right)$$ the top and bottom expressions of $$y$$ are differentiated to give $$\frac{f^\prime (x)}{g^\prime (x)} = \frac{\frac{P}{N(N-P)}}{-\frac{1}{N^2}} = \frac{-PN^2}{N^2-NP},$$ and $$\lim_{N\to \infty} y = \lim_{N\to \infty}\frac{-PN^2}{N^2-NP} = -P.$$ Recall $$y = log_e \left( 1-\frac{P}{N} \; \right) ^ N$$, and as it required that $$\left( 1-\frac{P}{N} \; \right) ^ N = \frac{1}{2}$$, it follows that in the limit $$N\to \infty$$ $$-P = log_e{\frac{1}{2}} \Rightarrow P=log_e(2).$$ This result is as expected from a differential equation model of radioactive decay. The process is modelled by $$\frac{dA}{dt} = - \lambda A,$$ where $$\lambda$$ is the rate of change constant. The solution to this is $$A(t) = A_0 e^{-\lambda t}.$$ If the time period is one half life then when $$t = 1$$, $$A = 0.5A_0$$, so that $$e^{-\lambda} = \frac{1}{2} \Rightarrow -\lambda = log_e(\frac{1}{2}) \Rightarrow \lambda = log_e(2).$$