Decay & Ingrowth


In the radio­chem­istry business one often needs to cal­cu­late decay and in­growth factors. A decay factor is the ex­pected frac­tion of a radio­nuclide re­main­ing after a speci­fied time t has elapsed, or, equiva­lently, the prob­abil­ity that an atom of that radio­nuclide will sur­vive at least until time t with­out decay­ing. Cal­cu­lating a decay fac­tor is fairly trivial. It equals eλt, where λ is the decay constant for the radio­nuclide. It’s as sim­ple as that.

In what follows it will be con­venient to have a sym­bol for the func­tion t  eλt. We will use Dλ, where the D ob­vi­ously means “decay.”



The problem becomes more in­ter­est­ing when one needs to cal­cu­late an ingrowth factor, de­fined as the ex­pected amount of a decay prod­uct that will exist at a later time be­cause of in­growth from a speci­fied an­cestor. The value of the in­growth fac­tor de­pends on whether one de­fines the amount of a nuclide in terms of the num­ber of atoms or the activ­ity. We will use the num­ber of atoms here, but the con­ver­sion to an ac­tiv­ity basis is sim­ple (because A = λN).

Consider a sim­ple non-branching decay chain:

S1 → S2 → S3 → … → Sn

where each Si denotes a state in the chain, with decay constant λi. The origi­nal an­ces­tor is state S1 and the final decay prod­uct of in­ter­est is state Sn. If the final state is stable, then λn = 0 s−1. There are many ex­am­ples of branching decay, where the final state can be reached by follow­ing any of sev­eral paths through a branch­ing decay dia­gram start­ing from the same initial state. (Math­emati­cally speak­ing the decay dia­gram is a directed acyclic graph, or DAG, because radio­nuclides decay only to lower energy states, not higher states.) The ingrowth fac­tor in this case is just a weighted sum of the in­growth fac­tors for the individual paths, with each path’s weight equal to its rela­tive prob­abil­ity of being taken (the prod­uct of the branch­ing frac­tions along the path). Each of these paths looks like the sim­ple linear chain shown above. So, we will focus on non-branching decay.

The ingrowth factor IF, when defined in terms of numbers of atoms, equals the prob­abil­ity that an atom in state S1 at time 0 will be in state Sn after time t has elapsed. The prob­abil­ity of this event (En) can be cal­cu­lated from the dis­tri­bu­tions of the decay times, T1, T2, …, Tn, of the n states. Each time Ti has an expo­nen­tial dis­tri­bu­tion, with a prob­abil­ity density func­tion (pdf) given by


The pdf for a sum of independent decay times T1 + T2 + + Tn is the con­volu­tion of the pdfs of the individual times.


where the con­vo­lu­tion oper­ator is defined for any two func­tions f (t) and g(t) by


The operator is asso­cia­tive and commuta­tive. Associa­tivity implies that if we have func­tions f1(t), f2(t), …, fn(t), we can write (f1 f2 ∗ ⋯ ∗ fn)(t) with­out ambiguity. Commutativity implies that (f g)(t) = (g f)(t).

Note: The defi­ni­tion im­plies that (f g)(0) = 0. In fact the first n − 2 de­riva­tives of a con­vo­lu­tion (f1 f2 ∗ ⋯ ∗ fn)(t) are also zero at t = 0.

An atom initially in state S1 will be in state Sn after time t has elapsed iff (if and only if) it decays through the first n − 1 states before time t and then remains in state Sn till time t without decay­ing further. So, the event En occurs iff T1 + T2 + + Tn − 1t and Tn > t (T1 + T2 + + Tn − 1). So, the atom ingrowth factor IFatom is equal to the prob­abil­ity Pr[En], given by:


Notice that if Sn is stable, then λn is zero and IFactivity = 0, but IFatom ≠ 0 un­less the in­growth time is zero.

In order to cal­cu­late either of these fac­tors, one needs to be able to cal­cu­late the con­vo­lu­tion of the ex­po­nen­tial func­tions Dλ1(t), Dλ2(t), …, Dλn(t). The theo­reti­cal so­lu­tion to this prob­lem is well-known. The dif­fi­culties arise when one tries to im­ple­ment the al­go­rithm effi­ciently for a com­puter in a man­ner that avoids severe round-off error.

If all the decay con­stants are dis­tinct (the usual case in the real world), the con­vo­lu­tion is given explicitly by


which is easily coded in a pro­gram­ming lan­guage like C, C++, C#, Java, or Fortran. (For the case where some of the λi are repeated, see below.) Equation () may look familiar to those who have dealt with the ingrowth prob­lem before.

Note: The value of the con­vo­lu­tion is zero at t = 0, as men­tioned above, although this fact may not be im­medi­ately ob­vi­ous from Equa­tion (). The func­tion’s first n − 2 de­riva­tives are also zero at t = 0.

Rounding Error

Consider the sim­plest non­trivial ex­ample, where n = 2. In this case Equa­tion () simplifies to the following:


This convolu­tion can be used when cal­cu­lating the in­growth of an im­medi­ate daughter prod­uct from a parent. Even in this sim­ple ex­am­ple there is the pos­sibil­ity of large round­ing error. Suppose t is very small. In this case the nu­mera­tor evalu­ates to the dif­fer­ence be­tween two num­bers that are nearly equal to 1. The most straight­forward im­ple­men­ta­tion of the equa­tion in a com­puter lan­guage will lead to an un­nec­es­sary loss of ac­cu­racy, pos­sibly even return­ing zero in situa­tions where the true value, although it may be small, is well within the machine’s range of posi­tive floating-point values.

An obvious fix in this sim­ple case when t is small is to use the Maclaurin series for the con­vo­lu­tion (see below). How­ever, there is an even easier fix in any pro­gram­ming lan­guage that has a li­brary func­tion for the hyper­bolic sine (sinh).


where λ = (λ1 + λ2) / 2, and b = (λ2λ1) / 2.

Warning: Do not use this trick in Micro­soft Excel. Excel, or at least the 2007 ver­sion on my home computer, uses a naïve imple­men­ta­tion of sinh. It cal­cu­lates sinh with the follow­ing equa­tion, which fails for ex­tremely small values of the argument:
The Numbers spread­sheet pro­gram by Apple uses a smarter implementation. Newer versions of Excel may also do it better.

The C language and some im­ple­men­ta­tions of its de­riva­tives also pro­vide a li­brary func­tion expm1(x), which cal­cu­lates ex − 1 ac­cu­rately. The con­vo­lu­tion equa­tion can be writ­ten using expm1 as follows:


where d = λ2λ1. If you have the expm1 func­tion in your li­brary, you should use it here.

For many of the most com­mon in­growth ex­amples in the radio­chem­istry lab, either of the two pre­ceding equa­tions will suf­fice. If there were no more com­pli­cated situa­tions than this, the prob­lem would be com­pletely solved; how­ever, there are more com­pli­cated exam­ples where n > 2. For in­stance one may want to cal­cu­late the in­growth fac­tor for a de­scend­ant of a very long-lived radio­nuclide, where some of the inter­medi­ate states are long-lived and some are short-lived. This is another situa­tion that can lead to severe round-off error. (E.g., try cal­cu­lat­ing 1 year’s worth of 222Rn in­growth from 238U.)

Equation () has the po­ten­tial to pro­duce ex­cessive round-off error when­ever there is a clus­ter of dis­tinct values λit that are too close to­gether. This can hap­pen when the values are large, but it is more com­mon when they are small, as in the sim­ple (n = 2) ex­am­ple de­scribed above when t is small, or in the 238U ex­am­ple, where sev­eral of the half-lives are very long. When­ever λit and λjt are small and almost equal — but not quite — Equa­tion () has large posi­tive and nega­tive terms that almost cancel each other — but not quite. The round­ing error in this case may be larger than the actual value.

The Maclaurin Series

When t is small enough, the Maclaurin series for the con­vo­lu­tion gives fast and ac­cu­rate results. This series can be writ­ten as follows:


where (n)k is the Pochhammer symbol, or the rising factorial, (n)k = (n)(n + 1)(n + 2)(n + k − 1), and where hk(x1, x2, …, xn) denotes the com­plete homo­geneous sym­met­ric poly­nomial of de­gree k, in n vari­ables, de­fined by


The rising factorial is easy to evalu­ate in­cre­ment­ally for suc­ces­sive values of k, be­cause (n)k + 1 = (n)k (n + k). Fortunately hk(x1, x2, …, xn) is also easy to evalu­ate incrementally.

  • h0(x1, x2, …, xn) = 1
  • hk + 1(x1) = x1k + 1 = x1hk(x1)
  • hk + 1(x1, x2, …, xi + 1) = hk + 1(x1, x2, …, xi) + xi + 1hk(x1, x2, …, xi + 1)

So, the Maclaurin series is easy to code. It also has the virtue of con­verg­ing for all values of t, but the con­ver­gence is rapid enough to be prac­tical for com­pu­ta­tion only if all the λit are some­what smaller than 1. There are many pos­sible sce­narios where some of the λit are too large to give rapid con­ver­gence of the Maclaurin series while others are so small (and there­fore close together) that they lead to un­accept­able round-off error in Equa­tion ().

When all the λit are small enough, the con­vo­lu­tion is ap­proxi­mated well by the first one or two non­zero terms of the Maclaurin series.


The (atom) ingrowth factor is then ap­proxi­mated by


An example of a scenario where this ap­proxi­ma­tion is use­ful is the cal­cu­la­tion of in­growth in the decay chain 234U → 230Th → 226Ra, where all three half-lives are fairly long. Using Equation () as writ­ten, with or­di­nary double-precision arith­metic, to cal­cu­late the in­growth of 226Ra over a time interval shorter than a year is likely to pro­duce a larger round­ing error than using just the first two terms of the Maclaurin series.

Note: By comparing the co­effi­cients in the Maclaurin series above with the de­riva­tives of Equation () evalu­ated at t = 0, one finds that

If all the λit are ap­proxi­mately equal but none is small, the fol­low­ing iden­tity allows one to com­pute the con­vo­lu­tion rapidly using a Maclaurin series (e.g., with c = λ).


However, such a situa­tion is un­likely to arise in prac­tice. So, it seems that the most ob­vi­ous approach to the prob­lem of round-off error won’t com­pletely solve it.

The Laplace Transform

Before exploring other op­tions, it will help to intro­duce a power­ful tool for deal­ing with the prob­lem: the Laplace transform. If f (t) is a func­tion of a non-negative real vari­able t, the Laplace trans­form of f (t) is the func­tion F(s) defined by


where s and F(s) are gen­erally con­sidered to be complex.

If F(s) is the Laplace trans­form of a function f (t), as shown above, and if F(s) is mero­morphic with poles c1,c2,…,cn, which may be either real or com­plex, then its in­verse Laplace trans­form is


where Ress = ci{estF(s)} means the residue of estF(s) at the pole ci.

Many useful facts about the func­tion f (t) = (Dλ1 Dλ2 Dλn)(t) can be derived by alge­braic manipu­la­tion of the trans­formed func­tion F(s) fol­lowed by ap­pli­ca­tion of the inverse trans­form. Although the same facts can be de­rived in other ways, the deri­va­tions using the Laplace trans­form tend to be easier.

Instead of applying the equa­tions shown above for the Laplace trans­form and its in­verse, it usually helps to re­mem­ber some com­mon func­tions and their trans­forms — or have access to a table of these — and to re­mem­ber that both the trans­form and its in­verse are linear.


The Laplace transform of an expo­nential func­tion Dλ(t) = eλt is par­ticu­larly simple:


By setting λ = 0, we see that


The Laplace transform of a con­vo­lu­tion (as defined above) is just the prod­uct of the trans­forms of the in­di­vid­ual functions.


where F(s) and G(s) are the Laplace trans­forms of f (t) and g(t), re­spec­tively. So, the Laplace trans­form of our con­vo­lu­tion is the rational function:


Applying the Laplace Transform

The trans­formed func­tion can be ma­nipu­lated algebraically to pro­duce new ex­pres­sions for the origi­nal con­vo­lu­tion, and some of these new ex­pres­sions lead to al­go­rithms for com­pu­ta­tion with less round-off error. Some­times the new ex­pres­sion is an in­fi­nite series, each term of which in­volves a con­vo­lu­tion of ex­po­nen­tial func­tions where each λi from a prob­lematic clus­ter has been re­placed by a sin­gle value μ, which might be one of the origi­nal λi from the clus­ter or per­haps the clus­ter av­er­age. Sub­sti­tuting a single value μ for sev­eral of the λi pro­duces a con­vo­lu­tion in which the decay con­stants are not all dis­tinct and to which Equation () doesn’t apply.

So, now we need the more gen­eral form of Equation (), which doesn’t require all the decay con­stants to be dis­tinct. For lack of a bet­ter no­ta­tion, de­fine Gμ, n(t) to be the n-fold con­vo­lu­tion of the single expo­nen­tial func­tion Dμ(t).*


The Laplace transform of Gμ, n(t) is


Our next goal is to derive ex­pres­sions that can be used to cal­cu­late the func­tion g(t) de­fined by:


where each mi is a posi­tive integer. The Laplace trans­form of g(t) is:


The function G(s) is mero­morphic, because it is a rational func­tion. In gen­eral, if F(s) is the Laplace trans­form of a func­tion f (t) and F(s) is mero­morphic with (real or com­plex) poles c1, c2, …, cn of orders m1, m2, …, mn, re­spec­tively, then


In our problem the function G(s) has poles μ1, μ2, …, μn. So,


where Ai, k is given by:


The expression for Ai, k doesn’t look very use­ful for com­pu­ta­tion, but with a little work one gets the following:


where the bullet nota­tion mx indi­cates that an ar­gu­ment x is repeated m times in the ar­gu­ment list. With another abuse of no­ta­tion, we can re­write the pre­ced­ing equa­tion as:


Then we can write the more gen­eral form of Equation ():

decay-gt2.png (∗∗)

Of course Equation (∗∗) has the same round-off issues as Equation () when the μi t are clus­tered together, but it works well enough when the μi t are well sepa­rated, and un­like the Maclaurin series, it doesn’t have slow con­ver­gence (and is less likely to over­flow) when the μi t are large.

Note: The fact that g(0) = 0 implies:

Other Series

Some relevant facts about the complete homo­geneous poly­nomials:

where each ki ≥ 0 and each zi .
where σ is any permutation of the set {1, 2, …, n}.
decay-sum-hk/png,for z < 1
decay-sum-hkn/png,if each zi < 1
decay-sum-hkmn/png,if each zi < 1

Using these facts we can easily derive the Maclaurin series for f (t).



With a few changes we can derive other series. Suppose for each λi we choose a μi, which may or may not be equal to λi. Then F(s) can be rewritten as follows:


where the ki are non-negative. Then

decay-eq3.png (∗∗∗)

This series is essentially a multi­vari­able Taylor series in the vari­ables λ1, λ2, …, λn. It can be made useful for calculation if we choose the μi care­fully. For rapid con­ver­gence the dif­fer­ences | μitλit | must be small (at least smaller than 1). The series will be sim­pler if there is only one or at most a few distinct μi for which μiλi, but the number of the μi needed to keep the round-off error small is really de­ter­mined by the number of λi clus­ters, because the μi must be well sepa­rated. If there is only one clus­ter, the result­ing series will be sim­ple enough. Say the λi are arranged in increas­ing order, and the first r of them form a clus­ter. We can let


and μi = λi for i > r. Then write f (t) as follows:


If there is more than one clus­ter of decay con­stants λi, then a more com­pli­cated series is needed, but in many situa­tions Equation (∗∗∗∗) suffices.

Given that μ = (λ1 + λ2 + + λr) / r, it is im­por­tant to note that the second term of this series (k = 1) is zero. So, con­ver­gence can­not be tested until k ≥ 2.

Implementing the series effi­ciently requires an effi­cient al­go­rithm for up­dating (Gμ, r + k Dλr + 1 Dλn)(t) for suc­ces­sive values of k. If λr + 1, λr + 2, , λn are dis­tinct, then



decay-W.png and decay-Cik.png for i = r + 1, r + 2, …, n

and where for any N and any x1, x2, …, xp,


Notice that W does not depend on k; so, it is cal­cu­lated only once. Also notice that


The algorithm can employ a sim­ple vec­tor to hold the values of Cr + 1, k , Cr + 2, k , …, Cn, k and after each Ci, k is used, replace it with Ci, k + 1. Most of the work of incre­ment­ing k is done in the re­cal­cu­la­tion of φr + k − 1, but there is an incre­mental al­go­rithm for this cal­cu­la­tion, and its com­plex­ity remains con­stant as k increases. Define


Then φ0(t; x1, x2, …, xp) = 1. For N ≥ 0,


and for i = 1, 2, …, p,


So, the time required to re­cal­cu­late φN(t; x1, x2, …, xp) for each new value of N is pro­por­tional to p but not to N.

The Siewers Algorithm

I discovered the Siewers algorithm — fortunately or unfortunately — only after do­ing a lot of work on the prob­lem my­self. This algo­rithm is based on the fact that

Siewers-basis.png, for λiλj,

where the check over a symbol (Ď) indicates a func­tion omitted from the con­vo­lu­tion. I actually knew this identity already but had not noticed that it could be used to develop an effi­cient algo­rithm. At first glance it appears to be O(2n), but it can be imple­mented more effi­ciently than that.





So, if we can calculate the function F reliably, we have solved our original problem too.

Notice that the value of F(x1, x2, …, xn) does not depend on the order of the argu­ments. So, re­arrange the argu­ments if nec­es­sary to ensure x1 x2 xn. The Siewers algo­rithm now calcu­lates F(x1, x2, …, xn) using the recursion


if the difference xnx1 is large enough (greater than some value δ). If the dif­fer­ence is too small, then


where c = (x1 + xn) / 2, and where the Maclaurin series is used for F(x1c, x2c, …, xnc).

Very small values of δ ensure rapid and accu­rate con­ver­sion of the Maclaurin series, while larger values tend to minimize round-off error in the re­cursive cal­cu­la­tion (∗∗∗∗∗). The most appro­priate choice of δ balances these two con­sidera­tions. A value between 1 and 2 is likely to give ac­cept­able results. (I get slightly bet­ter ac­cu­racy with 2 than with 1.)

The only remaining issue is how to avoid un­nec­es­sary re­cal­cula­tion of inter­mediate quan­tities of the form F(xi, xi + 1, …, xj) that are needed more than once by the algo­rithm. For s = 1 to n − 1, and for i = 1 to n + 1 − s, if either xi + sxi > δ or xi + s − 1xi − 1 > δ, then calculate


using the most appropriate method. If s = 1, then of course the most appro­priate method is given by


If s > 1, and if xi + s − 1xi > δ, then


Otherwise, calculate


where c = (xi + xi + s − 1) / 2, using the Maclaurin series.

Finally, calculate and output f1(n), which equals F(x1, x2, …, xn).

In fact a single array of length n can be used to hold the values fi(s) at each step of the iteration (s = 1 to n). The algo­rithm cal­cu­lates the follow­ing quantities as necessary:

i = 1 i = 2 i = n − 1 i = n
s = 1 f1(1) = F(x1) f2(1) = F(x2) fn − 1(1) = F(xn − 1) fn(1) = F(xn)
s = 2 f1(2) = F(x1, x2) f2(2) = F(x2, x3) fn − 1(2) = F(xn − 1, xn)
s = n − 1 f1(n − 1) = F(x1, x2, …, xn − 1) f2(n − 1) = F(x2, x3, …, xn)
s = n f1(n) = F(x1, x2, …, xn)

The quantities in each row after the first are calculated either (a) recursively, using two of the quan­tities from the pre­ceding row or (b) using the Maclaurin series. The last row holds the final result.

Final Note

Recall Equation (∗∗∗):


This equation would be most use­ful if we had an effi­cient incre­mental al­go­rithm to calculate the sum


for increasing values of k. An incre­mental al­go­rithm is easy enough to find — it’s making it effi­cient that’s hard. (Earlier we obtained effi­ciency with the assump­tion that there was only one dis­tinct value of μi for which μiλi.) There is a sim­ple incre­mental al­go­rithm to calculate


and in fact there is a very effi­cient incre­mental al­go­rithm to calculate


which follows from the some­what un­expected fact that


Recalculating this sum when increment­ing k requires only a multi­pli­ca­tion by t and a division by the new value of k.


Although this sum is easy to cal­cu­late, we un­for­tunately have no need for it. There may be a sim­ple and effi­cient al­go­rithm for the sum we do need, but this pro­gram­mer doesn’t know it yet.

The Siewers algorithm looks best.

Algorithms: A few explicit algorithms for some of the calculations described above

Example: 238U → 222Rn — to be completed

Application to mean decay or ingrowth during counting — not com­pleted, but sim­ple (include D0(t) 1 in the con­vo­lu­tion and divide the final result by t)

Uncertainty of a decay or ingrowth factor — to be completed (theo­reti­cal value ob­tained by propa­gating un­cer­tainties of time and decay con­stants would be mis­lead­ing in the presence of large and un­known round­ing error)

See also Non-Poisson CountingPDF — The deri­va­tion of equa­tions for the index of dis­per­sion (J factor) in non-Poisson count­ing meas­ure­ments makes use of ingrowth fac­tors for short-lived decay products (e.g., 222Rn progeny). For a Power­Point summary see Non-Poisson Count­ing Statistics, or What’s This J Factor All About?.

Note: (2009-11-10) It’s been more than a year and I still haven’t pro­vided many of these exam­ples.

* The symbol G is used here because the func­tion Gλ, n(t) is so closely related to the pdf for a gamma dis­tri­bu­tion. We choose to use the func­tion Gλ, n(t) instead of the gamma pdf itself so that we can allow λ to be zero, or even negative.