physor2018-scopatz.tex

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass[letterpaper]{physor2018}
%
%  various packages that you may wish to activate for usage
\usepackage{tabls}
\usepackage{cites}
\usepackage{epsf}
\usepackage{appendix}
\usepackage{ragged2e}
\usepackage[top=1in, bottom=1.in, left=1.in, right=1.in]{geometry}
\usepackage{enumitem}
\setlist[itemize]{leftmargin=*}
\usepackage{caption}
\captionsetup{width=1.0\textwidth,font={bf,normalsize},skip=0.3cm,within=none,justification=centering}
\usepackage{xspace}
\usepackage{amsmath}

%\usepackage[justification=centering]{caption}

\newcommand{\ith}{i$^{\mathrm{th}}$\xspace}
\newcommand{\jth}{j$^{\mathrm{th}}$\xspace}
\newcommand{\pth}{p$^{\mathrm{th}}$\xspace}
\newcommand{\qth}{q$^{\mathrm{th}}$\xspace}
\newcommand{\Cth}{C$^{\mathrm{th}}$\xspace}

%
% Define title...
%
\title{BINARY FORMULATION OF DECAY EQUATIONS \\
       WITH LIMITING CASES \& CRAM COMPARISON}
%
% ...and authors
%
\author{%
  % FIRST AUTHORS j
  %
  \textbf{Anthony M. Scopatz$^1$ and Aaron Meurer$^1$}\\
  $^1$University of South Carolina \\
  541 Main Street \#011, \\
  Columbia, SC 29208 \\
  \\
  \url{scopatz@cec.sc.edu}
}
%
% Insert authors' names and short version of title in lines below
%
\newcommand{\authorHead}      % Author's names here use et al. if more than 3
           {Scopatz and Meurer}
\newcommand{\shortTitle}      % Short title here (Shorten to fit all into a single line)
           {Limiting Cases of the Binary Formulation of Decay Equations}
\renewcommand{\thetable}{\Roman{table}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   BEGIN DOCUMENT
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\maketitle
\justify

\begin{abstract}
This paper presents a binary formulation of the Bateman equations
for radioactive decay. This work includes various limiting cases
that arise in real-world calculations due to ill-behaved nuclear data
or floating point arithmetic issues. A binary formulation is advantageous
due to it speed, often orders of magnitude faster than a traditional solver.
After the description of the new binary method,
a verification study that compares the binary solver
to the Chebyshev Rational Approximation Method (CRAM). The comparison
between CRAM and the binary solver demonstrates that the method is
well within floating point precision for most nuclides and times.
\end{abstract}
\keywords{bateman, decay, cram, pyne}

\section{INTRODUCTION}
\label{sec-intro}
The Bateman equations \cite{bateman1910solution} governing radioactive decay are an
important subexpression
of the generalized transmutation equations. In many cases, it is desirable to compute
decay on its own, outside of the presence of any neutron or photon field. In this
case, radioactive decay is a function solely of intrinsic physical parameters:
half-lives, branch ratios, and fission product yields. This document recasts the
Bateman equations into a form that is better suited for computation than
the traditional expression. Unlike previous work along these lines
\cite{scopatz2015decay}, the formulation presented here aims to remove the
assumption that the nuclear data is well-behaved.

The fundamental nuclear data is not well-behaved. There are degeneracies in
half-lives. There are cases where prominent nuclides (such as U-238) are
almost stable relative to floating point precision. Since these issues
cannot be tackled in the data itself, they must be handled by the
mathematics. The standard Bateman equations are acceptable for the
vast majority of nuclide decays. However, this work presents modifications
to the standard equations for those limiting cases where usual solution
falls apart.

In \S\ref{sec-normal-method}, this paper derives the Bateman equation
limiting cases using the traditional decay constant formulation.
\S\ref{sec-binary} continues on to recast these in a binary formulation.
The binary formulation has computation speed advantageous over the
traditional decay constant method. \S\ref{sec-impl-spec-approx} presents implementation
specific approximations that are made with the equations found in
\S\ref{sec-binary}. The assumptions listed in \S\ref{sec-impl-spec-approx}
have been reduced in number from previous work \cite{scopatz2015decay}.
\S\ref{sec-verify} verifies the binary methodology by comparing its
results to two alternative solvers. Finally, \S\ref{sec-conc} summarizes
the work and the results presented here.

\section{BATEMAN EQUATIONS WITH LIMITING CASES}
\label{sec-normal-method}
The canonical expression of the Bateman equations for a decay chain
proceeding from a nuclide $A$ to a nuclide $Z$ at time
$t$ following a given pathway may be seen in Eq. (\ref{canon-bateman}) below \cite{CETNAR2006640}:
\begin{equation}
\label{canon-bateman}
N_C(t) = \frac{N_1(0)}{\lambda_C} \cdot \gamma \cdot \sum_{i=1}^C \lambda_i c_{i} e^{-\lambda_i t}
\end{equation}
The symbols in Eq. (\ref{canon-bateman}) are described in
Table \ref{nomenclature} in \S\ref{sec-nomen}. Additionally, this expression requires that
$c_{i}$ be defined as in Eq. (\ref{c_i-def}):
\begin{equation}
\label{c_i-def}
c_i = \prod_{j=1,i\ne j}^C \frac{\lambda_j}{\lambda_j - \lambda_i}
\end{equation}
Furthermore, the total chain branch ratio $\gamma$ is defined as the product of the
branch ratio between any two species \cite{harr2007precise}:
\begin{equation}
\label{gamma-total}
\gamma = \prod_{i=i}^{C-1} \gamma_{i \to i+1}
\end{equation}

Minor modifications to Eq. (\ref{canon-bateman}) are needed for terminal species:
the first nuclide of a decay chain and the ending stable species. For the first nuclide in a chain,
setting $C=1$ will reduce the Bateman equations to the simple form see in Eq. (\ref{n-single}):
\begin{equation}
\label{n-single}
N_C(t) = N_1(0) e^{-\lambda_1 t}
\end{equation}
For stable species, the appropriate equation is derived by taking the limit
such that the decay constant of the stable nuclide ($\lambda_C$) goes to
zero. Note that every $c_i$ contains exactly one $\lambda_C$
in the numerator which cancels with the $\lambda_C$ in the denominator
in front of the summation:
\begin{equation}
\label{stable-lim}
\lim_{\lambda_C \to 0} N_C(t) = N_1(0)  \gamma \left[e^{-0t} + \sum_{i=1}^{C-1} \lambda_i \left(\frac{1}{0 - \lambda_i} \prod_{j=1,i\ne j}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_i} \right) e^{-\lambda_i t} \right]
\end{equation}
\begin{equation}
\label{n-stable}
N_C(t) = N_1(0)  \gamma \left[1.0 - \sum_{i=1}^{C-1} \left(\prod_{j=1,i\ne j}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_i} \right) e^{-\lambda_i t} \right]
\end{equation}

From here, certain chains have intermediate nuclides which are \emph{almost} stable
(within the limits of floating point precision). For example, decaying
from Es-254 to Po-210 goes through U-238, which is very close to stable relative to all of the
other nuclides in the chain. This can trigger some terms of the Bateman equation to
underflow, overflow, or generate \texttt{NaN} values. Such situations must be avoided,
if at all possible. To handle these floating point issues, call $p$ the index of the nuclide
that is almost stable. The Bateman equations can be reduced by the
observation that $\lambda_p \ll \lambda_{i\ne p}$ after the $p$-term is separated out
from the summation:
\begin{equation}
\label{almost-stable-0}
\frac{N_C(t)}{N_1(0)} = \frac{\gamma}{\lambda_C}\sum_{i\ne p}^C \left[\lambda_i \frac{\lambda_p}{\lambda_p - \lambda_i}
                                                    \left(\prod_{j\ne i,p}^C \frac{\lambda_j}{\lambda_j - \lambda_i}\right)
                                                    e^{-\lambda_i t}\right]
                       + \frac{\gamma}{\lambda_C} \lambda_p \left(\prod_{j\ne p}^C \frac{\lambda_j}{\lambda_j - \lambda_p} \right) e^{-\lambda_p t}
\end{equation}
\begin{equation}
\label{almost-stable-2}
   \frac{N_C(t)}{N_1(0)} = \frac{\gamma}{\lambda_C}\sum_{i\ne p}^C \left[\lambda_i \frac{\lambda_p}{- \lambda_i}
                                                        \left(\prod_{j\ne i,p}^C \frac{\lambda_j}{\lambda_j - \lambda_i}\right)
                                                        e^{-\lambda_i t}\right]
                           + \frac{\gamma}{\lambda_C} \lambda_p \left(\prod_{j\ne p}^C \frac{\lambda_j}{\lambda_j}\right) e^{-\lambda_p t}
\end{equation}
\begin{equation}
\label{almost-stable-3}
   \frac{N_C(t)}{N_1(0)} = \frac{-\gamma\lambda_p}{\lambda_C}\sum_{i\ne p}^C \left[
                                                        \left(\prod_{j\ne i,p}^C \frac{\lambda_j}{\lambda_j - \lambda_i}\right)
                                                        e^{-\lambda_i t}\right]
                           + \frac{\gamma\lambda_p}{\lambda_C} e^{-\lambda_p t}
\end{equation}
Eq. (\ref{almost-stable-3}) is valid when the last
nuclide in the chain is itself unstable. However, when the last nuclide is stable,
both the \pth (almost stable nuclide) and
the \Cth (last and stable nuclide) must be removed, split off from
the summation, and handled separately. As previously, take $\lambda_C \to 0$ and
$\lambda_p \ll \lambda_{i\ne p,C}$.
\begin{equation}
\label{last-stable-and-almost-stable-0}
\begin{split}
   \frac{N_C(t)}{N_1(0)} = & \frac{\gamma}{\lambda_C}\sum_{i\ne p}^{C-1} \left[\lambda_i \frac{\lambda_C}{\lambda_C - \lambda_i} \frac{\lambda_p}{\lambda_p - \lambda_i}
                                                        \left(\prod_{j\ne i,p}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_i}\right)
                                                        e^{-\lambda_i t}\right] \\
                           & + \frac{\gamma}{\lambda_C} \lambda_p \frac{\lambda_C}{\lambda_C - \lambda_p} \left(\prod_{j\ne p}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_p} \right) e^{-\lambda_p t} \\
                           & + \frac{\gamma}{\lambda_C} \lambda_C \frac{\lambda_p}{\lambda_p - \lambda_C} \left(\prod_{j\ne p}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_C} \right) e^{-\lambda_C t}
\end{split}
\end{equation}
\begin{equation}
\label{last-stable-and-almost-stable-1}
\begin{split}
   \frac{N_C(t)}{N_1(0)} = & \gamma\sum_{i\ne p}^{C-1} \left[\frac{\lambda_i \lambda_p}{(\lambda_C - \lambda_i)(\lambda_p - \lambda_i)}
                                                        \left(\prod_{j\ne i,p}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_i}\right)
                                                        e^{-\lambda_i t}\right] \\
                           & + \frac{\gamma\lambda_p}{\lambda_C - \lambda_p} \left(\prod_{j\ne p}^{C-1} \frac{\lambda_j}{\lambda_j} \right) e^{-\lambda_p t}\\
                           & + \frac{\gamma\lambda_p}{\lambda_p - \lambda_C} \left(\prod_{j\ne p}^{C-1} \frac{\lambda_j}{\lambda_j} \right) e^{-\lambda_C t}
\end{split}
\end{equation}
\begin{equation}
\label{last-stable-and-almost-stable-2}
   \frac{N_C(t)}{N_1(0)} = -\gamma\sum_{i\ne p}^{C-1} \left[\left(\prod_{j\ne i}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_i}\right) e^{-\lambda_i t}\right]
                           + \frac{\gamma\lambda_p}{\lambda_C - \lambda_p} e^{-\lambda_p t}
                           + \frac{\gamma\lambda_p}{\lambda_p - \lambda_C} e^{-\lambda_C t}
\end{equation}
\begin{equation}
\label{last-stable-and-almost-stable-3}
   \frac{N_C(t)}{N_1(0)} = -\gamma\sum_{i\ne p}^{C-1} \left[\left(\prod_{j\ne i}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_i}\right) e^{-\lambda_i t}\right]
                           + \frac{\gamma\lambda_p}{\lambda_C - \lambda_p} \left(e^{-\lambda_p t} - e^{-\lambda_C t}\right)
\end{equation}
\begin{equation}
\label{last-stable-and-almost-stable-4}
   \frac{N_C(t)}{N_1(0)} = -\gamma\sum_{i\ne p}^{C-1} \left[\left(\prod_{j\ne i}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_i}\right) e^{-\lambda_i t}\right]
                           -\gamma e^{-\lambda_p t} + \gamma
\end{equation}

Lastly, there remains degenerate case where two nuclides in a chain have the same exact half-lives.
This unfortunate situation arises out of the fundamental nuclear data. For example,
\begin{enumerate}
\item Lu-153 $\alpha$-decays to Tm-149, both nuclides have half-lives of 0.9 seconds.
\item Ir-172 $\alpha$-decays to Re-168, both nuclides have half-lives of 4.4 seconds.
\item Pt-183 $\alpha$-decays to Os-179, both nuclides have half-lives of 390 seconds.
\end{enumerate}
Call the degenerate nuclides the \pth and \qth species. To prevent underflow, overflow,
and \texttt{NaN} values, we must separate these nuclides out of the summation
and then take the limit as $\lambda_q \to \lambda_p$.
\begin{equation}
\label{pq-same-0}
\begin{split}
   \frac{N_C(t)}{N_1(0)} = & \frac{\gamma}{\lambda_C}\sum_{i\ne p,q}^{C} \left[\lambda_i \left(\prod_{j\ne i}^{C} \frac{\lambda_j}{\lambda_j - \lambda_i}\right) e^{-\lambda_i t}\right] \\
                           & + \frac{\gamma}{\lambda_C} \lambda_p \frac{\lambda_q}{\lambda_q - \lambda_p} \left(\prod_{j\ne p,q}^{C} \frac{\lambda_j}{\lambda_j - \lambda_p} \right) e^{-\lambda_p t} \\
                           & + \frac{\gamma}{\lambda_C} \lambda_q \frac{\lambda_p}{\lambda_p - \lambda_q} \left(\prod_{j\ne p,q}^{C} \frac{\lambda_j}{\lambda_j - \lambda_q} \right) e^{-\lambda_q t}
\end{split}
\end{equation}
\begin{equation}
\label{pq-same-1}
   \frac{N_C(t)}{N_1(0)} = \frac{\gamma}{\lambda_C}\sum_{i\ne p,q}^{C} \left[\lambda_i \left(\prod_{j\ne i}^{C} \frac{\lambda_j}{\lambda_j - \lambda_i}\right) e^{-\lambda_i t}\right]
                           + \frac{\gamma\lambda_p^2}{\lambda_C} \left(\prod_{j\ne p,q}^{C} \frac{\lambda_j}{\lambda_j - \lambda_p} \right)
                             \lim_{\lambda_q\to\lambda_p}\frac{e^{-\lambda_p t} - e^{-\lambda_q t}}{\lambda_q - \lambda_p}
\end{equation}
\begin{equation}
\label{pq-same-2}
   \frac{N_C(t)}{N_1(0)} = \frac{\gamma}{\lambda_C}\sum_{i\ne p,q}^{C} \left[\lambda_i \left(\prod_{j\ne i}^{C} \frac{\lambda_j}{\lambda_j - \lambda_i}\right) e^{-\lambda_i t}\right]
                           + \frac{\gamma\lambda_p^2}{\lambda_C} \left(\prod_{j\ne p,q}^{C} \frac{\lambda_j}{\lambda_j - \lambda_p} \right) t e^{-\lambda_p t}
\end{equation}
Eq. (\ref{canon-bateman})-(\ref{pq-same-2}) thus represent a complete system of equations that
handle the relevant limiting cases which arise in normal computation. The first is to minimize

\section{BINARY FORMULATION OF BATEMAN EQUATIONS WITH LIMITING CASES}
\label{sec-binary}
There are two main strategies which may be used to construct a version of these equations that
is better suited to computation on modern processors. The first strategy is to group and
precompute all known constants. The second is to use more efficient operations.

The first strategy aims to minimizing the number of
operations that must be performed to achieve the same result. This can be achieved
by grouping constants together and pre-calculating them. Such a pre-calculation is possible
because all of the nuclear data is known at compile time. This saves the computer from
having to perform the same operations at run time every execution.
The Bateman equations may thus be expressed as a simple sum of exponentials,
\begin{equation}
\label{sum-of-exp}
    N_C(t) = N_1(0) \sum_{i=1}^C k_{i} e^{-\lambda_i t}
\end{equation}
where the coefficients $k_i$ are defined as follows:
\begin{itemize}
\item \textbf{Single Nuclide in Chain:}
\begin{equation}
\label{k-single}
    k_i = 1
\end{equation}
\item \textbf{Last Nuclide Unstable:}
\begin{equation}
\label{k-last-unstable}
    k_i = \frac{\gamma}{\lambda_C} \lambda_i \prod_{j\ne i}^C \frac{\lambda_j}{\lambda_j - \lambda_i}
\end{equation}
\item \textbf{Last Nuclide Stable:}
\begin{equation}
\label{k-last-stable-0}
\begin{split}
    k_{i\ne C} = & -\gamma \prod_{j\ne i}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_i} \\
    k_C = & \gamma
\end{split}
\end{equation}
\item \textbf{Last Nuclide Unstable and \pth Almost Stable:}
\begin{equation}
\label{k-last-unstable-p-almost-stable-0}
\begin{split}
    k_{i\ne p} = & -\frac{\gamma\lambda_p}{\lambda_C} \prod_{j\ne i,p}^C \frac{\lambda_j}{\lambda_j - \lambda_i} \\
    k_p = & \frac{\gamma\lambda_p}{\lambda_C}
\end{split}
\end{equation}
\item \textbf{Last Nuclide Stable and \pth Almost Stable:}
\begin{equation}
\label{k-last-stable-p-almost-stable-0}
\end{equation}
\begin{equation}
\label{k-last-unstable-p-almost-stable-1}
\begin{split}
    k_{i\ne p,C} = & -\gamma \prod_{j\ne i}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_i} \\
    k_p = & -\gamma \\
    k_C = & \gamma
\end{split}
\end{equation}
\item \textbf{Half-life Degeneracy Between \pth and \qth Nuclides:}
\begin{equation}
\label{k-degen-pq-0}
\begin{split}
    k_i = & \frac{\gamma}{\lambda_C} \lambda_i \prod_{j\ne i}^C \frac{\lambda_j}{\lambda_j - \lambda_i} \\
    k_p = & t \frac{\gamma\lambda_p^2}{\lambda_C} \prod_{j\ne p,q}^C \frac{\lambda_j}{\lambda_j - \lambda_p} \\
    k_q = & 0
\end{split}
\end{equation}
\end{itemize}
If the $k_i$ above are computed at run time, then this expression results in much more
computational effort that than the original Bateman equations since $\gamma/\lambda_C$
are brought into the summation. However, when $k_i$ are pre-calculated,
many floating point operations are saved by avoiding explicitly computing $c_i$.
Note that the one exception to the precomputation is that $k_p$ in Eq. (\ref{k-degen-pq-0})
has an additional $t$ term which must be multiplied at run time.

The second strategy is to note that computers are much better at computing powers of
2 then then any other base, even $e$. Thus the \texttt{exp2(x)} function, or $2^x$,
is faster than the natural exponential function \texttt{exp(x)}, $e^x$.
This results in a run time savings between 20-25\% on modern processors.
Since the core of the Bateman equations
are exponentials, it is worthwhile to recast these equations if possible. Luckily, the decay constant
provides an intrinsic mechanism to convert to base-2:
\begin{equation}
\label{n-base-2-0}
    N_C(t) = N_1(0) \sum_{i=1}^C k_{i} e^{-\lambda_i t}
\end{equation}
\begin{equation}
\label{n-base-2-1}
    N_C(t) = N_1(0) \sum_{i=1}^C k_{i} e^{\frac{-\ln(2)\cdot t}{t_{1/2,i}}}
\end{equation}
\begin{equation}
\label{n-base-1}
    N_C(t) = N_1(0) \sum_{i=1}^C k_{i} 2^{\frac{-t}{t_{1/2,i}}}
\end{equation}
This expression can be further collapsed by defining $a$ to be the precomputed
exponent values:
\begin{equation}
\label{a-def}
    a_i = \frac{-1}{t_{1/2,i}}
\end{equation}
Thus, the final form of the binary representation of the Bateman equations are
as follows:
\begin{equation}
\label{final-form-a}
    N_C(t) = N_1(0) \sum_{i=1}^C k_{i} 2^{a_i t}
\end{equation}
where the $k_i$ are as listed above.  However, for practical purposes, it is better to
compute the $k_i$ from half-lives rather than decay constants.  This is because they
provide less floating point error, fewer opportunities to underflow, overflow, or be \texttt{NaN} or
$\infty$ valued. Furthermore, half-lives provide a better mechanism for detecting stability.
Thus, alternatively, the $k_i$ may be computed as,
\begin{itemize}
\item \textbf{Single Nuclide in Chain:}
\begin{equation}
\label{k-bin-single}
    k_i = 1
\end{equation}
\item \textbf{Last Nuclide Unstable:}
\begin{equation}
\label{k-bin-unstable}
    k_i = \gamma t_{1/2,i}^{C-2} t_{1/2,C} \prod_{j\ne i}^{C} \frac{1}{t_{1/2,i} - t_{1/2,j}}
\end{equation}
\item \textbf{Last Nuclide Stable:}
\begin{equation}
\label{k-bin-stable-0}
\begin{split}
    k_i = & -\gamma t_{1/2,i}^{C-2} \prod_{j\ne i}^{C-1} \frac{1}{t_{1/2,i} - t_{1/2,j}} \\
    k_C = & \gamma
\end{split}
\end{equation}
\item \textbf{Last Nuclide Unstable and \pth Almost Stable:}
\begin{equation}
\label{k-bin-unstable-p-almost-0}
\begin{split}
    k_{i\ne p} = & -\frac{\gamma t_{1/2,C}}{t_{1/2,p}} t_{1/2,i}^{C-2} \prod_{j\ne i,p}^C \frac{1}{t_{1/2,i} - t_{1/2,j}} \\
    k_p = & \frac{\gamma t_{1/2,C}}{t_{1/2,p}}
\end{split}
\end{equation}
\item \textbf{Last Nuclide Stable and \pth Almost Stable:}
\begin{equation}
\label{k-bin-stable-p-almost-0}
\begin{split}
    k_{i\ne p,C} = & -\gamma t_{1/2,i}^{C-2} \prod_{j\ne i}^{C-1} \frac{1}{t_{1/2,i} - t_{1/2,j}} \\
    k_p = & -\gamma \\
    k_C = & \gamma
\end{split}
\end{equation}
\item \textbf{Half-life Degeneracy Between \pth and \qth:}
\begin{equation}
\label{k-bin-pq-degen-0}
\begin{split}
    k_i = & \gamma t_{1/2,i}^{C-2} t_{1/2,C} \prod_{j\ne i}^{C} \frac{1}{t_{1/2,i} - t_{1/2,j}} \\
    k_p = & t \gamma\ln(2) t_{1/2,p}^{C-4} t_{1/2,C}  \prod_{j\ne p,q}^C \frac{1}{t_{1/2,p} - t_{1/2,j}} \\
    k_q = & 0
\end{split}
\end{equation}
\end{itemize}
Again, the one exception to the precomputation is that $k_p$ in Eq. (\ref{k-bin-pq-degen-0})
has an additional $t$ term which must be multiplied at run time.
Still, with precomputed $k$, $a$, and the \texttt{exp2()} function, this
formulation minimizes the number of floating point operations while completely
preserving physics. No assumptions were made aside from the Bateman equations
themselves.
Note that it is not possible to reduce the number of operations further.

\section{IMPLEMENTATION SPECIFIC ADJUSTMENTS}
\label{sec-impl-spec-approx}
The formulation presented in \S\ref{sec-binary} holds generally for any decay chain.
However, certain modifications are used in practice to reduce the number of chains and terms
that are calculated. The two modifications within PyNE's solver are:
\begin{enumerate}
\item Decay alphas are not treated as He-4 production.
\item The $k_i$ and $a_i$ are filtered to reject terms where,
   \begin{equation}
   \label{filter-cond}
   \frac{|k_i|}{\max(|k_i|)} < 10^{-16}
   \end{equation}
   This filtering prevents excessive
   calculation from species which do not significantly contribute to
   end atom fraction. The threshold $10^{-16}$ was chosen as
   because it is a reasonable na\"ive estimate of floating point error after
   many operations. Note that we may filter on only the $k_i$ because
   $2^{a_i t} \le 1$.  That is, the exponential component can only
   reduce the magnitude of a term, not increase it.
\end{enumerate}
In principle, these two statements are reasonable. However, they
may preclude desired behavior by users. In such a situation, these
assumptions should be revisited.


\section{METHOD VERIFICATION}
\label{sec-verify}
This section contains two comparisons of the binary Bateman equations
as presented in Eq. (\ref{final-form-a})-(\ref{k-bin-pq-degen-0}). These
comparisons verify that the method described here is correct up to
floating point error. The first comparison simply verifies that the
initial nuclide decays properly. The second and more involved comparison
uses the Chebyshev Rational Approximation Method (CRAM)
for transmutation \cite{pusa2010computing,pusa2012correction} as the point
of comparison.

In the comparisons here, 3491 nuclides were decayed from an initial unit
concentration. Decays took place on a log-uniform time grid from $10$ to $10^{20}$
seconds. Including $t=0$, there were 40 points in the time grid. This grid
spans time domains that are generally of interest to nuclear
engineers. This nuclide \& time comparison thus provides a set of 139640
decayed compositions to compare.

\subsection{$N_1(t)$ Comparison}
\label{subsec-n1}
For this first comparison, note that the solution for the initial nuclide is
the relatively simple $N_1(t) = e^{-\lambda_1 t}$, as $N_1(0) = 1$. Since this is analytic and
easy to compute, the binary Bateman solver ought to compute this same
value to within floating point error.  Call $b(t)$ the concentration of the
same initial nuclide as computed by the binary Bateman solver.

\begin{figure}[!htb]
  \centering
  \includegraphics[scale=0.80]{./n1-diff.pdf}
  \caption{The difference between known decay of initial nuclide and the
           binary Bateman solution for all 3491 nuclides. The red line
           is the average of the absolute value of difference with
           error bars representing 1 standard deviation. The gray lines
           are curves for a sub-selection of example nuclides which comprise
           the mean.}
  \label{fig-n1}
\end{figure}

Fig.~\ref{fig-n1} shows the absolute value of the difference between
the traditional Bateman solution $e^{-\lambda_1 t}$ and the binary
solver described above. The gray lines show the difference as a function
of time for individual, representative nuclides in the set analyzed.
The red curve in Fig.~\ref{fig-n1} displays the mean difference over all
times, along with one standard deviation error bars.

As can be seen from Fig.~\ref{fig-n1}, the mean difference tends to zero for
higher decay times.  Moreover, this mean difference curve itself seems to
follow an exponential decay. Importantly, all values in Fig.~\ref{fig-n1}
fall within what could be reasonably expected for floating point error for
\texttt{double} precision computation. Additionally, the gray curves
demonstrate that for some individual nuclides, the difference has a minimum
modal bound. These can be as high as $2.2\times 10^{-16}$. These higher
difference modes persist at long time scales, even as the mean difference
tends to zero.

\subsection{CRAM Comparison}
\label{subsec-cram}
Recent versions of PyNE (v0.5.4+)
now include a both a binary Bateman solver and a CRAM solver.
The CRAM solver is code generated using SymPy \cite{10.7717/peerj-cs.103},
while the binary Bateman solver is code generated using a custom algorithm.
For the CRAM solver, a decay-only matrix can be constructed from the
fundamental nuclear data. This decay-only matrix, thus, does not include
any transmutation from neutron interactions. This usage is distinct from
most CRAM applications, which are typically focused on depletion calculation.

To create a fair comparison, both solvers were generated from the same
nuclear data. Half-lives and branch
ratios came from the May 1$^{\mathrm{st}}$, 2017 release of ENSDF
\cite{Bhat1992,ensdfmaintained}. Fission product yield data was obtained
from published WIMS-D results \cite{aldama2003wims}. Therefore, any differences
in the results between the binary Bateman solver and the CRAM
solver can be ascribed to differences in methodology, not differences in data.

However, it should be noted that CRAM is an approximation and the Bateman
equations are the analytic solution. Therefore, we should expect differences
between these two methods up to the approximation order of CRAM. In practice,
the CRAM order is chosen to be at the level of floating point precision.
Since we generally are dealing with \texttt{double} data (64-bit), rather than
\texttt{float} data (32-bit), the floating point error is expected to be
approximately $10^{-16}$ for an initial unit atom density $N_1(0) = 1$.
Thus the CRAM approximation order is chosen as 16 in the comparison here.

\begin{figure}[!htb]
  \centering
  \includegraphics[scale=0.80]{./diff-hist.pdf}
  \caption{A histogram of the absolute value of the differences between the
           binary Bateman solution $b(t, i\to j)$
           and the CRAM solution  $r_{16,16}(t, i\to j)$ for all decays
           from an initial \ith nuclide to a \jth nuclide for all times.}
  \label{fig-diff-hist}
\end{figure}

Call $b(t, i\to j)$ the binary Bateman concentration at time $t$ of the \jth nuclide
coming from an initial unit concentration.
Also denote $r_{16,16}(t, i\to j)$ as this same value, but instead as computed
by a $16^{\mathrm{th}}$ order CRAM approximation. Fig.~\ref{fig-diff-hist}
displays a histogram of the absolute value of the difference between
$b(t, i\to j)$ and $r_{16,16}(t, i\to j)$ for all $i \to j$ for all times.
This histogram is binned by decade from $10^{-32}$ to $1$.

As can be seen in Fig.~\ref{fig-diff-hist}, the overwhelming majority of
decay compositions which are computed fall within floating point precision
($\le 10^{-16}$). For those times and nuclides which are not within precision,
many of the remaining differences are within what might be considered
expected floating point error after many arithmetic operations ($\le 10^{-8}$).

This still leaves a handful of nuclides and times (around 600)
where the difference is beyond acceptable floating point errors. Since differences
in data may be discounted, these errors are likely caused by CRAM solves.
CRAM happens to have relatively inaccurate eigenvalues for the decay problem. Additionally,
CRAM performs
far more arithmetic operations when inverting the decay matrix than the binary Bateman
solver performs. This occasionally leads to catastrophic floating point cancellation
issues and may also result in the high errors seen here.

Still, even though certain cases have high errors, no error exceeds 0.1.
Additionally, the nuclides which have differences in excess of $10^{-4}$
all come from relatively negligible decay chains. For example, this
list includes transitions from Tl-177 to W-168, Es-254 to
Pb-206, and Bh-267 to Pb-207. These initial species are rarely relevant
to traditional nuclear engineering applications.


\section{CONCLUSIONS}
\label{sec-conc}

In summary, the binary formulation of the traditional Bateman equations
is the mathematical basis for a fast and correct solver for radioactive decay problems. Because this
method minimizes run time arithmetic operations and takes advantage of
\texttt{exp2()} speed over \texttt{exp()}, the solver here can be orders of
magnitude faster than alternatives. Furthermore, in comparison to CRAM used
for the same problem, the binary solver is correct to within floating point
error for almost all cases. Those differences with CRAM which are
unreasonably high stem from the approximations that CRAM makes. These
approximations involve inverting many matrices. CRAM, thus, requires many
more arithmetic operations and much longer run times.

The advantage of deriving the limiting cases for the binary Bateman solver
seen in Eq. (\ref{final-form-a})-(\ref{k-bin-pq-degen-0}) is that in enables
a reduction in the number of assumptions made by previous work
\cite{scopatz2015decay}. This in turn lead to a drastic reduction in the
error of the previous solver. Properly handling degenerate half-lives and
almost-stable nuclides improves the quality of the solution for the whole
decay chain in which these issues arise

\section*{NOMENCLATURE}
\label{sec-nomen}

Variables used throughout this paper have are described in Table \ref{nomenclature}.

\begin{table}[!htb]
\centering
\caption{Symbols, Units, \& Meaning}
\label{nomenclature}
\begin{tabular}{|l|l|l|}
\hline
\textbf{Symbol} & \textbf{Units} & \textbf{Meaning} \\
\hline
$C$                  & unitless  & Length of the decay chain \\
$i$                  & unitless  & Index for \ith species, on range $[1, C]$ \\
$j$                  & unitless  & Index for \jth species, on range $[1, C]$ \\
$p$                  & unitless  & Index for \pth species, which is special \\
$q$                  & unitless  & Index for \qth species, whose half-life is the same as $p$ \\
$t$                  & seconds   & Time \\
$t_{1/2,i}$          & seconds   & Half-life of the \ith species \\
$\lambda_i$          & 1/seconds & Decay constant of \ith species, $\ln(2)/t_{1/2,i}$ \\
$N_i(t)$             & 1/cm$^3$  & Number density of the \ith species at time t \\
$\gamma_{i \to i+1}$ & unitless  & Branch ratio between \ith species and its child \\
$\gamma$             & unitless  & Total branch ratio for a chain \\
$k_i$                & unitless  & Coefficients allowing for sum-of-exponentials \\
$a_i$                & 1/seconds & Negative inverse of half-life, $-1/t_{1/2,i}$ \\
\hline
\end{tabular}
\end{table}

\section*{ACKNOWLEDGEMENTS}
\label{sec-ack}
This work was funded in part by the University of South Carolina. The authors
wish to express our gratitude.


% Bibliography
\setlength{\baselineskip}{12pt}
\bibliographystyle{physor}
\bibliography{physor}

\end{document}