Decay & Ingrowth

Decay

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.”

Ingrowth

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:

where each S_i denotes a state in the chain, with decay constant λ_i. The origi­nal an­ces­tor is state S₁ and the final decay prod­uct of in­ter­est is state S_n. If the final state is stable, then λ_n = 0 s⁻¹. 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 S₁ at time 0 will be in state S_n after time t has elapsed. The prob­abil­ity of this event (E_n) can be cal­cu­lated from the dis­tri­bu­tions of the decay times, T₁, T₂, …, T_n, of the n states. Each time T_i 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 T₁ + T₂ + ⋯ + T_n 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 f₁(t), f₂(t), …, f_n(t), we can write (f₁ ∗ f₂ ∗ ⋯ ∗ f_n)(t) without ambiguity. Commutativity implies that (f ∗ g)(t) = (g ∗ f)(t).

An atom initially in state S₁ will be in state S_n 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 S_n till time t without decay­ing further. So, the event E_n occurs iff T₁ + T₂ + ⋯ + T_n − 1 ≤ t and T_n > t − (T₁ + T₂ + ⋯ + T_n − 1). So, the atom ingrowth factor IF_atom is equal to the prob­abil­ity Pr[E_n], given by:

Then the activity ingrowth factor IF_activity is equal to

Notice that if S_n is stable, then λ_n is zero and IF_activity = 0, but IF_atom ≠ 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.

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 λ = (λ₁ + λ₂) / 2, and b = (λ₂ − λ₁) / 2.

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 e^x − 1 ac­cu­rately. The con­vo­lu­tion equa­tion can be writ­ten using expm1 as follows:

where d = λ₂ − λ₁. 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 ²²²Rn in­growth from ²³⁸U.)

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 ²³⁸U 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 h_k(x₁, x₂, …, x_n) 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 h_k(x₁, x₂, …, x_n) is also easy to evalu­ate incrementally.

•  h₀(x₁, x₂, …, x_n) = 1
•  h_k + 1(x₁) = x₁^k + 1 = x₁ h_k(x₁)
•  h_k + 1(x₁, x₂, …, x_i + 1) = h_k + 1(x₁, x₂, …, x_i) + x_i + 1 h_k(x₁, x₂, …, x_i + 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 ²³⁴U → ²³⁰Th → ²²⁶Ra, 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 ²²⁶Ra 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.

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 c₁, c₂, …, c_n, which may be either real or com­plex, then its in­verse Laplace trans­form is

where Res_{s = ci}{e^st F(s)} means the residue of e^st F(s) at the pole c_i.

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 m_i 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 c₁, c₂, …, c_n of orders m₁, m₂, …, m_n, re­spec­tively, then

In our problem the function G(s) has poles −μ₁, −μ₂, …, −μ_n. So,

where A_i, k is given by:

The expression for A_i, 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 m • x 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 (∗):

Of course Equation (∗∗) has the same round-off issues as Equation (∗) when the μ_it are clus­tered together, but it works well enough when the μ_it 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 μ_it are large.

Other Series

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

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 k_i are non-negative. Then

This series is essentially a multi­vari­able Taylor series in the vari­ables λ₁, λ₂, …, λ_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 μ = (λ₁ + λ₂ + ⋯ + λ_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

where

and where for any N and any x₁, x₂, …, x_p,

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 C_{r + 1, k}, C_{r + 2, k}, …, C_n, k and after each C_i, k is used, replace it with C_{i, 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 φ₀(t; x₁, x₂, …, x_p) = 1. For N ≥ 0,

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

So, the time required to re­cal­cu­late φ_N(t; x₁, x₂, …, x_p) 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

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(2ⁿ), but it can be imple­mented more effi­ciently than that.

Define

Then

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

Notice that the value of F(x₁, x₂, …, x_n) does not depend on the order of the argu­ments. So, re­arrange the argu­ments if nec­es­sary to ensure

The Siewers algo­rithm now calcu­lates F(x₁, x₂, …, x_n) using the recursion

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

where c = (x₁ + x_n) / 2, and where the Maclaurin series is used for F(x₁ − c, x₂ − c, …, x_n − c).

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(x_i, x_i + 1, …, x_j) 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 x_i + s − x_i > δ or x_{i + s − 1} − x_i − 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 x_{i +s −1} − x_i > δ, then

Otherwise, calculate

where c = (x_i + x_{i + s − 1}) / 2, using the Maclaurin series.

Finally, calculate and output  f₁⁽ⁿ⁾, which equals F(x₁, x₂, …, x_n).

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

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

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: ²³⁸U → ²²²Rn — to be completed

Application to mean decay or ingrowth during counting — not com­pleted, but sim­ple (include D₀(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 Counting” — 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., ²²²Rn progeny). For a Power­Point summary see Non-Poisson Count­ing Statistics, or “What’s This J Factor All About?”.

Comments/suggestions

top

3 Feb 2010 Keith