Please bear with my long winded description.

It is classical to model radioactive decay of some particles $A \rightarrow A^\ast$ by the differential equation \begin{align} \frac{\mathrm dN_A}{\mathrm dt}=-\lambda N_A \end{align} where $N_A$ refers to the number of particles in state $A$ and $\lambda$ is the rate of decay. If we initially start with $N_0$ number of particles in state $A$ then by solving the above differential equation we see that \begin{align} N_A(t) = N_0\mathrm e^{-\lambda t} \end{align} gives us an explicit number of particles in state $A$ as a function of time $t$. In particular, we also see that \begin{align} \lambda = \frac{\ln 2}{t_{1/2}} \end{align} where $t_{1/2}$ is the half-life. However, the model is completely deterministic, meaning no probability required.

However, quoting Wikipedia

Radioactive decay is a stochastic (i.e. random) process at the level of single atoms, in that, according to quantum theory, it is impossible to predict when a particular atom will decay, regardless of how long the atom has existed.

But collectively, the mean number of particles in state $A$ is well modeled by the above deterministic model. However, the difficulty here is to articulate the connection between the deterministic model and the stochastic model in a rigorous manner.

Let $T$ be a continuous positive random variable (RV) that denotes the lifetime of a single (microscopic) particle in state $A$. It is natural to make the assumption that in the stochastic model the decay of a particle in state $A$ is independent of the time that the particle spends in state $A$. With this assumption and a little bit of math (which is not important for the discussion), one can show that $T$ has an exponential distribution with a parameter $\lambda$, i.e. the PDF of $T$ is given by \begin{align} p_\lambda (t) = \begin{cases} \lambda\mathrm e^{-\lambda t} & \text{ if } t \geq 0,\\ 0 & \text{ if } t<0. \end{cases} \end{align} and the cumulative distribution function (CDF) is given by \begin{align} P(t\geq T) = 1-\mathrm e^{-\lambda t} \end{align} for $t \geq 0$ and 0 otherwise. In particular, we see that the probability of the particle in the decayed state $A^\ast$ at time $t$ is precisely $1-\mathrm e^{-\lambda t}$ and the probability in of the particle in the state $A$ is $\mathrm e^{-\lambda t}$. Hence at each fixed time $t$ we have a Bernoulli distribution. If we make another assumption that a collection of the radioactive particles decay independently, then we see that given $N_0$ particles at time $t$ the probability that $N$ particles survives and $N_0-N$ particles decay is given by the binomial distribution \begin{align} P(N_t=N) = \frac{N_0!}{N!(N_0-N)!}(1-\mathrm e^{-\lambda t})^{N_0-N}\mathrm e^{- \lambda N t} \end{align} By standard computations, we also see that the mean of the binomial distribution is $\operatorname{E}[N_t]= N_0\mathrm e^{-\lambda t}$ and the variance is $\operatorname{Var}[N_t]= N_0\mathrm e^{-\lambda t}(1-\mathrm e^{-\lambda t})$ which leads to the relative standard deviation \begin{align} \frac{\sqrt{\operatorname{Var}[N_t]}}{\operatorname{E}[N_t]} = \frac{1}{\sqrt{N_0}}\left(\mathrm e^{\lambda t}- 1\right)^{1/2}. \end{align}

I have two questions.

Question 1: Looking at the relative standard deviation, it seems to me that for very large time $t$ the deterministic model seems to be a poor approximation of the stochastic model since the relative standard deviation grows exponentially in time. Is my assertion valid?

I actually disagree with the above calculation but have no valid argument against it.

Question 2: My interest in the above analysis comes from my desired to give a rigorous justification for using Monte-Carlo simulation of radioactive decay to arrive at the above deterministic model. However, I'm stuck after the above analysis. Could someone provide me with some reference or explanation to rigorous justify the connection?

• For the Monte Carlo approach look at the Gillespie method (J. Phys. Chem. v 81, p2340, 1977) and his book 'Markov Processes' publ. Academic Press. May 23 '18 at 8:12

The average disintegration rate is given by $\bar n=\sum_{n=0}^{N_0} nP(n)$ where $P$ is the binomial distribution. If $M$ is the average number of events in time $t$ then $\displaystyle M=N_0(1-e^{-\lambda t})$. For times smaller than the lifetime $1/\lambda$ the $\displaystyle e^{-\lambda t} \to 1-\lambda t$ then $M=N_0\lambda t$ and as the disintegration rate is $M/t$ this leads to $\displaystyle \frac{dN}{dt}=-\lambda N$. Thus provided the measuring time-bins are small compared to the lifetime this rate law should hold.

As the decay proceeds the number of counts becomes small, perhaps $10^6$ at short times to just 1 or 2 at long times then of course the signal to noise becomes huge, zero counts in some time bins for example and 1 or 2 counts in others. (This becomes an issue in analysing data, see Anal. Chem., 2003, 75 (16), pp 4182–4187).

The expected variance is, as you quote, $V=\sigma^2=N_0e^{-\lambda t}(1-e^{-\lambda t})$ and can be rewritten as $V=Me^{-\lambda t}$ where $M$ is the mean number of events recorded in time $t$. As the time bins are normally short so that $\lambda t <<1$ then $\sigma=\sqrt{M}$ so that $\displaystyle \frac{\sigma}{M}=\frac{1}{\sqrt{M}}$

• Does not look like a good method to me, the method is not a general method that can be done quickly with a pocket calculator. When the half life is very high when compared with the time of observation such as decays from 1 mg of radium-226 during one hour. The count rate is given by lambda N. But for dealing with a system where the observation time is similar or longer than one half life then the Bateman equations are normally the most common method used in the nuclear / radioactivity world. The method of using time bins is normally only used for things like multicompartment animal models. May 23 '18 at 16:05
• This answer has nothing to do with the Bateman equations. By 'time bins' I mean the total time divided by the number of data points that are measured in an experiment or number calculated if that is what you are doing. Photon counting the emission from molecules is entirely similar to radioactive decay and a few thousand points is normally sufficient. Of course it does depends on how many orders of magnitude in intensity you need to follow the decays over, 4 decades would be typical. The matrix method I commented on (see below) is equivalent to the Bateman but based on linear algebra. May 23 '18 at 16:59
• Have you read the masters thesis from the USAF, I am not very interested in the case of the single radionuclide with a moderate half life or hours or days. It is rather simple to measure its halflife. The interesting cases are mixtures of two or more radionuclides or complex chains of radionuclides. The time bin method is a poor first choice as it is often so computationally demanding, I see it as a last resort method. I would always try with the Bateman equations first, if they fail I would then try a "time bins" method May 23 '18 at 19:15
• I had a quick look at the thesis, it describes standard known methods, nothing new; the moments method is known to be very poor, for example, and the matrix method effectively the same as in my comment. I think that you misunderstand what I mean by time-bins: you have to calculate data at discrete time points these are 'time bins' as they are discrete separated by $\Delta t$ May 24 '18 at 7:13
• Having run excel based simulations of radioactive decay many times, the problem I see is that anything other than the Bateman equation method requires a lot of cells to calculate a value. Can I ask you are you doing the simple system of a radionuclide decaying to a stable product or are you doing decay chains, it is not clear to me which of the two are you trying to deal with. Out of interest have you ever done radioactivity work ? May 24 '18 at 16:29

There are several methods of dealing with the problem of radioactive decay calculations. I rather like the Bateman equations as a method. They can be used to calculate the amount of a radionuclide at time t, given the amounts at zero time and the half lives of the radionuclides in question.

This includes decay chains such as A --> B --> C --> D --> E --> F etc etc

The weakness of the Bateman equation method is that the dynamic range of a digital computer causes errors to appear. The core of the problem (or nucleous of the problem) is that in a Bateman equation we have terms such as (Lambda1-Lambda2).

When the two lambda values are very close to each other then error due to the way in which excel (or even more specialist maths software) subtracts one number from another causes some errors to appear.

I have on occasion tried to use Bateman equations to model things like the conversion of nuclides into each other in a reactor core. Such as

U238 --> U239 --> Np239 --> Pu239 --> Pu240 --> Pu241 --> Am241 --> Am242

But what happens is that the further we get along the chain the worse the errors become, after a while a spreadsheet starts to give an unreasonable output such as predicting an almost infinite amount of a later radionuclide. For the early radionuclides in the chain the spreadsheet works perfectly.

In general the greater the difference between the halflives the better an excel sheet using Bateman equations works. For example if we were to try to consider the system.

Pu241 --> Am241 --> Np237

As the half lives of the three nuclides are very different the lambda values are very different. In case anyone does not know what a lambda value is, it is the decay constant used in the equation.

A = Ao exp (-lambda t)

We get a lambda value by dividing ln(2) by the half life.

Here in our Pu241, Am241, Np237 system it would work very nicely.

If on the otherhand we were to attempt to consider

Y-95 --> Zr-95 --> Nb-95

As the half lives of the Zr and Nb nuclides are rather similar then excel is less able to make a good prediction.

There is a thesis written on this subject which was done by a USAF officer (Captain Logan J. Harr) who considered the evolution of mixtures of fission products in the seconds after a short duration fission event (Nuclear bomb detonation). This is an interesting but rather complex problem as the decay chains are typically very long.

The Monte Carlo method has been considered by others both on Stack exchange and elsewhere. I am sure that with plenty of computational power that the Monte Carlo method will work well, the disadvantage of it I think is that it is just a brute force attack on a problem. It is a rather clever bute force attack but it is still a brute force attack.

Another method would be to copy the method used to deal with multicompartment models for things like plutonium and lead in a human. Here if we treat a human as being four compartments, we can assume a series of first order reactions take plutonium (or lead) from one compartment to another. If we set up a spreadsheet such that we calculate the amount of metal in each compartment by assuming that the flow is constant over a series of short time slices, then it is possible with short time slices to get a good prediction.

For example here is one paper on the subject.

• For a sequence of species such as you give as we usually are not interested in what the exact equations actually are, it is easier to use a matrix method to solve the equations. If $K$ is the matrix of rate constants and $C_0$ the initial values then $C=Xe^{-t\lambda}X^{-1}C_0$ where $X$ are eigenvectors, $\lambda$ are the eigenvalues of $K$ , and $e^{-t\lambda}$ a diagonal matrix of these values with zeros elsewhere. May 23 '18 at 8:04