# 5.1.3. Method of solution

ORIGEN solves the system of ordinary differential equations (ODEs) that describe nuclide generation, depletion, and decay,

where

\(N_{i}\) = amount of nuclide

i(atoms),\(\lambda_{i}\) = decay constant of nuclide

i(1/s),\(l_{\text{ij}}\) = fractional yield of nuclide

ifrom decay of nuclidej,\(\sigma_{i}\) = spectrum-averaged removal cross section for nuclide

i(barn),\(f_{\text{ij}}\) = fractional yield of nuclide

ifrom neutron-induced removal of nuclidej,\(\Phi\) = angle- and energy-integrated time-dependent neutron flux (neutrons/cm

^{2}-s), and\(S_{i}\) = time-dependent source/feed term (atoms/s).

Note that Eq. (5.1.1) has no spatial dependence and can be interpreted as either a solution at a point in space or the spatial average over some volume. The latter interpretation is preferred here, such that \(\Phi\) is the spatially averaged neutron flux magnitude, and all energy-dependence is embedded in the one-group flux-weighted average cross sections \(\sigma_{i}\) and reaction yields \(f_{\text{ij}}\). Eq. (5.1.1) is conveniently written in matrix form as

with a \(\mathbf{A}\) commonly referred to as the “transition matrix.” The representation of the transition matrix as \(\mathbf{A = A}_{\sigma}\Phi\mathbf{+}\mathbf{A}_{\lambda}\), where \(\mathbf{A}_{\sigma}\) is the part of the transition matrix containing reaction terms and \(\mathbf{A}_{\lambda}\) is the part containing decay terms, is convenient, as the numerical solution of this system of ODEs holds the reaction, flux, and feed terms constant over step \(n\),

over time step \(t_{n - 1} \leq t \leq t_{n}.\)

Adding a continuous removal process described with rate constant \(\lambda_{i,rem}\) simply modifies the decay constant, \(\lambda_{i} \rightarrow \lambda_{i} + \lambda_{i,rem}\), whereas a continuous feed process defines a nonzero component of the \({\overrightarrow{S}}\) vector.

ORIGEN can also compute the alpha, beta, neutron, and gamma emission spectra during decay. For the “stand-alone” ORIGEN calculations described here, the transition matrix is loaded from an ORIGEN binary library file (f33), which uses sparse-matrix storage to store one or more transition matrices. The f33s may be created using COUPLE, saved from TRITON depletion calculations, or interpolated using ARP from a set of precompiled f33s distributed with SCALE.

Results from ORIGEN calculations may be stored on a binary concentration file (f71), which facilitates transfer of isotopics to other codes in SCALE. The f71 file can also store calculated emission spectra. Within ORIGEN, the f71 can be used to restart calculations from an existing set of compositions.

## 5.1.3.1. Methodology

This section describes the methodology used in performing the following main functions:

generation of problem-dependent transition matrices,

solution of the system of depletion/decay equations,

conversion from power to flux (important for reactor applications),

calculation of emission spectra, and

interpolation of pregenerated sets of transition matrices.

### 5.1.3.1.1. Generation of Problem-dependent Transition Matrix

In the transition matrix \(\mathbf{A}\) from Eq. (5.1.2),
each matrix element \(a_{\text{ij}}\) is the first-order rate constant for the
formation of nuclide *i* from nuclide *j* given below.

The transition matrix coefficients for decay and reaction transitions are stored separately and reaction transitions are always stored with \(\Phi = 1\) and later during solution of the system, depending on the step-average flux level the actual transition matrix \({\mathbf{A}_{n}\mathbf{= A}}_{\sigma,n}\Phi_{n}\mathbf{+}\mathbf{A}_{\lambda}\) is reconstructed using step-average flux, \(\Phi_{n}\).

The decay coefficients \(l_{\text{ij}}\lambda_{j}\) and \(\lambda_{i}\) are generated directly from ORIGEN decay resource data. The reaction coefficients \(f_{\text{ij}}\sigma_{j}\) and \(\sigma_{i}\) are generated using the following two-stage procedure.

Calculate all removal cross sections \(\sigma_{i}\) and yields \(f_{\text{ij}}\), including isomeric branching ratios and fission yields, by folding provided flux spectrum \(\phi^{g}\) with multigroup cross sections from the ORIGEN reaction resource and energy-dependent fission yield data from the ORIGEN yield resource.

Overwrite specific removal cross sections and yields based on a provided multigroup cross section library (see the SCALE Nuclear Data Libraries chapter) and/or user-provided one-group cross sections and yields.

The second stage is optional, but it is important for cases which there is significant self-shielding because ORIGEN’s reaction resource assumes infinite dilution for its multigroup data. The decay, reaction, and yeld resources mentioned here are described in the ORIGEN Data Resources chapter. The collapse to a one-group cross section in either stage is given by

for reaction type \(r\), nuclide \(i\), and provided multigroup flux \(\phi^{g}\). Different reaction types are recognized by their ENDF MT numbers [SCALE Cross Section Libraries chapter] on the appropriate data resource For example, MT=16 is \(\left( n,2n \right),\) and MT=107 is \(\left( n,\alpha \right)\). The removal cross section \(\sigma_{i}\) is simply calculated as the sum over all relevant reactions for a particular nuclide, \(\sigma_{i} = \ \sum_{r}\sigma_{\text{ri}}\). This type of reaction-dependent multigroup data may be contained in either the data sources available in stage 1 or 2 above. However, only two types of data are expected to be available in stage 1 reaction resource data: (1) isomeric branching and (2) fission yields.

The energy-dependent isomeric branching that describes the yield of each excited level (metastable state) of a daughter nucleus is calculated in a similar way,

where \(m\) indicates the possible metastable states and the fractions always satisfy \(\sum_{m}f_{\text{rim}}^{\ } = 1\).

Fission product yields are typically tabulated at discrete neutron energies such as thermal (0.0253 eV), fission (500 keV), and high energy (14 MeV). The yield for each fissionable nuclide is calculated in stage 1 by linearly interpolating the tabulated data using the computed average energy of fission,

where \(\sigma_{\text{fi}}^{g}\) is the multigroup fission cross
section, and \({\overline{E}}^{g}\) is the average energy in
the group (simple midpoint energy used). In addition to generating
transition data for daughter/residual nuclides, the coefficients for
byproducts such as He-4/\(\alpha\) byproducts from
\(\left( n,\alpha \right)\) reactions are also retained in the
transition matrix and associated to an appropriate nuclide in the
system: hydrogen, deuterium, tritium, ^{3}He, or ^{4}He.

### 5.1.3.1.2. Solution of the Depletion/Decay Equations

ORIGEN includes two solver kernels that can solve the depletion/decay equations of Eq. (5.1.3):

a hybrid matrix exponential/linear chains method (MATREX) and

a Chebyshev Rational Approximation Method (CRAM).

They are described in the following sections.

#### 5.1.3.1.2.1. MATREX

Referring to the system of ODEs shown in Eq. (5.1.2) and setting the external feed/source \(S\left( t \right) = 0\), there is a formal solution by matrix exponential (an analog to the solution of a single ODE of this type by exponential),

where \(\overrightarrow{N}\left( 0 \right)\) is a vector of initial nuclide concentrations, by defining the series expansion of \(exp(\mathbf{A}t\mathbf{)}\) to be

with \(\mathbf{I}\) the identity matrix. Eq. (5.1.8) and Eq. (5.1.9) describe the matrix exponential method, which yields a complete solution to the problem. However, in certain instances related to limitation in computer precision, difficulties occur in generating accurate values of the matrix exponential function. Under these circumstances, alternative procedures using either the generalized Bateman equations [ORIGEN-Bat10] or Gauss-Seidel iterative techniques are applied. These alternative procedures will be discussed in further sections.

A straightforward solution of Eq. (5.1.8) and Eq. (5.1.9) would require storage of the complete transition matrix. To avoid excessive memory requirements, a recursion relation has been developed. Substituting Eq. (5.1.9) into Eq. (5.1.8),

one may recognize a recursion relation for a particular nuclide, \(N_{i}(t)\).

where the range of indices, *j*, *k*, *m*, is 1 to *M* for matrix
\(\mathbf{A}\) of size \(M \times M\). The result is a series of terms
that arise from the successive post-multiplication of the transition
matrix by the vector of nuclide concentration increments produced from
the computation of the previous terms. Within the accuracy of the series
expansion approximation, physical values of the nuclide concentrations
are obtained by summing a converged series of these vector terms. By
defining the terms \(C_{i}^{n}\left( t \right)\) as

the solution for \(N_{i}\left( t \right)\) is given as

The use of Eq. (5.1.12) and Eq. (5.1.13) requires storage of only two vectors–\({\overrightarrow{C}}^{n}\) and \({\overrightarrow{C}}^{n + 1}\) —-in addition to the current value of the solution. However, the series summation solution in Eq. (5.1.13) is not valid until a finite limit is identified which can achieve a reasonable accuracy, i.e.,

where \(n_{\text{term}}\) is the number of terms and \(\epsilon_{\text{trunc}}\) is the truncation error. The key is to split the nuclides into two sets: those that are long-lived and permit a rapid, accurate solution via Eq. (5.1.14), and those that are short-lived and require an alternate solution.

##### Solution for Long-Lived Nuclides

This section describes the various tests used to ensure that the summations indicated in Eq. (5.1.14) do not lose accuracy due to large changes in magnitudes or small differences between positive and negative rate constants. Nuclides with large rate constants (short-lived) are removed from the transition matrix and treated separately. For example, in the decay chain \(\mathrm{A\rightarrow B \rightarrow C}\), if the decay constant for B is large, a new rate constant is inserted in the matrix for \(\mathrm{A \rightarrow C}\). This technique was originally employed by Ball and Adams [ORIGEN-BA67]. The key to determining which transitions should be removed involves calculation of the matrix norm. The norm of matrix \(\mathbf{A}\) is defined by Lapidus and Luus [ORIGEN-LL67] as being the smaller of the maximum-row absolute sum and the maximum-column absolute sum,

To maintain precision in performing the summations of Eq. (5.1.14),
the matrix norm is used to balance the user-specified time step, *t*, with
the precision associated with the word len>h employed in the machine
calculation. The constraint on the matrix norm has been chosen as

The remainder of this section shows that this constraint serves two purposes.

It allows reasonable accuracy for a reasonable number (20–60) of matrix exponential terms.

It defines what “short-lived” means over a particular time step, dictating which concentrations must be solved by alternate means.

A relationship between *m* digits of machine precision and *p*
significant digits required in all results can be stated by the
following inequality:

In this particular series, the relationship may be represented as

or alternatively,

Lapidus and Luus have shown that the maximum term in the summation for
any element in the matrix exponential function cannot exceed
\(\frac{\left( \left\lbrack \mathbf{A} \right\rbrack t \right)^{n}}{n!},\ \)where
\(n\) is the largest integer not larger than
\(\left\lbrack \mathbf{A} \right\rbrack t\). For the constraint in
Eq. (5.1.16), this yields *n*=13 and yields limit
\(\frac{\left( \left\lbrack \mathbf{A} \right\rbrack t \right)^{n}}{n!} \approx 10^{5}\).
With \(e^{\left\lbrack \mathbf{A} \right\rbrack t} \approx 10^{6}\)
and standard double precision with *m=16*, Eq. (5.1.19) evaluates
to \(10^{11} \leq 10^{16 - j}\), which indicates that five significant
figures will be maintained in values as small as 10^{–6}. The
number of terms required to converge the matrix exponential series can
be investigated by a plot of the
\(\frac{\left| \left\lbrack \mathbf{A} \right\rbrack t \right|^{n}}{n!}e^{\left\lbrack \mathbf{A} \right\rbrack t}\)
as a function of term index *n* in Eq. (5.1.19), as shown in Fig. 5.1.1

The intersection between the black line in Fig. 5.1.1 and the various curves indicates the number of terms needed to achieve \(\epsilon_{\text{trunc}} \leq 0.1\%\). For example, with \(\left\lbrack \mathbf{A} \right\rbrack t = 13.8155\), \(n_{\text{term}} = 54\) is required, and with \(\left\lbrack \mathbf{A} \right\rbrack t = 13.8155/2\), approximately \(n_{\text{term}} = 29\) is required. This behavior has been used to develop the heuristic

Thus it has been shown that the limit imposed in Eq. (5.1.16) leads to a maximum of \(n_{\text{term}} = 54\) terms with \(\epsilon_{\text{trunc}} \leq 0.1\%\).

It remains to be shown that any arbitrary system can be modified so that it does not violate Eq. (5.1.16). Because the time step \(t\) is provided and fixed, \(\left\lbrack \mathbf{A} \right\rbrack t \leq \ - 2\ \ln(0.001)\) cannot be satisfied unless the system is modified. The physical nature of the system leads to \(\max_{j}{\sum_{i}\left| a_{\text{ij}} \right|} \leq \max_{j}{2|a_{\text{jj}}|}\) based on production rates equal to loss rates when both parent and daughter nuclide are included in the system. The maximum column sum in Eq. (5.1.15) can then be bounded by twice the maximum diagonal term, \(\max_{j}{2|a_{\text{jj}}|}\). Using this upper limit as the matrix norm and substituting into Eq. (5.1.16) yields

Rearranging Eq. (5.1.21) leads to the condition

which is used to mark nuclide *j* as a short-lived nuclide for this time
step, to be solved with linear chains instead of the series-based matrix
exponential. An alternative interpretation of the short-lived condition
can be made by rewriting Eq. (5.1.22) in terms of an effective
half-life, \(t_{1/2} = \frac{\ln\left( 2 \right)}{|a_{\text{jj}}|}\), which
results in
\(t_{1/2} < \frac{{- ln}\left( 2 \right)}{\ln\left( 0.001 \right)}t \approx 0.1t\).
In other words, when a nuclide’s effective half-life (including
destruction by both decay and reaction mechanisms) is less than 10% of
the time step, it can be considered short-lived.

Finally, as a note for applications where the nuclides of interest are in long transmutation chains, it has been found that the above algorithm may not yield accurate concentrations for those nuclides near the end of the chain that are significantly affected by those near the beginning of the chain. In these applications, specifying the minimum \(n_{\text{term}}\) as

where \(\Delta Z\) is the atomic number difference and \(\Delta A\) is the mass number difference, has been found to ameliorate the issue.

##### Solution for Short-Lived Nuclides

The condition in Eq. (5.1.22) forms the basis for declaring a nuclide short-lived, and its solution is found via solution of the nuclide chain equations. In conjunction with maintaining the transition matrix norm below the prescribed level, a queue is formed of the short-lived precursors of each long-lived isotope. These queues extend back up the several chains to the last preceding long-lived precursor. According to Eq. (5.1.22), the queues will include all nuclides whose effective half-lives are less than 10% of the time interval. A generalized form of the Bateman equations developed by Vondy [ORIGEN-Von62] is used to solve for the concentrations of the short-lived nuclides at the end of the time step. For an arbitrary forward-branching chain, Vondy’s form of the Bateman solution is given by,

where \(N_{1}\left( 0 \right)\) is the initial concentration of the first precursor, \(N_{2}\left( 0 \right)\) is that of the second precursor, etc.

As in Eq. (5.1.4), \(a_{\text{ij}}\) is the first-order rate constant, and \(d_{i} = {- a}_{\text{ii}}\) which is the magnitude of the diagonal element. Bell recast Vondy’s form of the solution through multiplication and division by \(\prod_{n = k}^{i - 1}d_{n}\) and rearranged to obtain

The first product over isotopes *n* is the fraction of atoms that
remains after the *kth* particular sequence of decays and captures. If
this product becomes less than 10^{-6}, the contribution of this
sequence to the concentration of nuclide *i* is neglected. Indeterminate
forms that arise when *di=dj* or *dn=dj* are evaluated using
L’Hôpital’s rule. These forms occur when two isotopes in a chain have
the same diagonal element.

Eq. (5.1.25) is applied to calculate all contributions to the “queue end-of-interval concentrations” of each short-lived nuclide from the initial concentrations of all others in the queue described above. It is also applied to calculate contributions from the initial concentrations of all short-lived nuclides in the queue to the long-lived nuclide that follows the queue, in addition to the total contribution to its daughter products. These values are appropriately applied either before or after the matrix expansion calculation is performed to correctly compute concentrations of long-lived nuclides and the long-lived or short-lived daughters. Eq. (5.1.25) is also used to adjust to certain elements of the final transition matrix, which now excludes the short-lived nuclides. The value of the element must be determined for the new transition between the long-lived precursor and the long-lived daughter of a short-lived queue. The element is adjusted so that the end-of-interval concentration of the long-lived daughter calculated from the single link between the two long-lived nuclides (using the new element) is the same as what would be determined from the chain including all short-lived nuclides. The method assumes zero concentrations for precursors to the long-lived precursor. The computed values asymptotically approach the correct value with successive steps through time. For this reason, at least five to ten time intervals during the decay of discharged fuel is reasonable, because long-lived nuclides have built up by that time.

If a short-lived nuclide has a long-lived precursor, an additional
solution is required. First, the amount of short-lived nuclide *i* due
to the decay of the initial concentration of long-lived precursor *j* is
calculated as

from Eq. (5.1.24), assuming \(e^{- d_{i}t} \ll \ e^{- d_{j}t}\).
However, the total amount of nuclide *i* produced depends on the
contribution from the precursors of precursor *j*, in addition to that
given by Eq. (5.1.26). The quantity of nuclide *j* not accounted for in
Eq. (5.1.26) is denoted by \(N_{j}'\left( t \right)\), the
end-of-interval concentration, minus the amount that would have remained
had there been no precursors to nuclide *j*:

Then the short-lived daughter and subsequent short-lived progeny are assumed to be in secular equilibrium with their parents, which implies that the time derivative is zero,

The queue end-of-interval concentrations of all the short-lived nuclides following the long-lived precursor are augmented by amounts calculated with Eq. (5.1.25). The concentration of the long-lived precursor used in Eq. (5.1.27) is that given by Eq. (5.1.26). The set of linear algebraic equations given by Eq. (5.1.28) is solved by the Gauss-Seidel iterative technique. This algorithm involves an inversion of the diagonal terms and an iterated improvement of an estimate for \(N_{i}(t)\) through the expression

Since short-lived isotopes are usually not their own precursors, this iteration often reduces to a direct solution.

##### Solution of the Nonhomogeneous Equation

The previous sections have presented the solution of the homogeneous
equation in Eq. (5.1.8), applicable to fuel burnup, activation, and
decay calculations. However, the solution of a nonhomogeneous equation
is required to simulate reprocessing or other systems that require an
external feed term, \(S\left( t \right) \neq 0\). The nonhomogeneous
equation is given in matrix form (assumed constant over a step *n*) as

for a fixed feed or removal rate, \({\overrightarrow{S}}\). A particular solution of Eq. (5.1.30) will be determined and added to the solution of the homogeneous equation given by Eq. (5.1.10). As before, the matrix exponential method is used for the long-lived nuclides, and solution by linear chains is used for the short-lived nuclides. Assume \(\overrightarrow{C}\) an arbitrary vector with which to test a particular solution of the form

Substituting Eq. (5.1.31) into Eq. (5.1.30) yields

in which the *k=0* term may be extracted from the LHS,

which allows the summations on the left and right to be easily shown equal. This proves the particular solution is indeed valid if the arbitrary vector is in fact the feed term \(\overrightarrow{C} = \overrightarrow{S}\). The solution to the nonhomogeneous problem is therefore (as a series),

For the second term in Eq. (5.1.34), a new recursion relation is developed for the particular solution in the same manner as was done for the homogeneous solution,

where

For the short-lived nuclides, the secular equilibrium equations are modified to become

The Gauss-Seidel iterative method is applied to determine the solution. The complete solution to the nonhomogeneous equation in Eq. (5.1.31) is given by the sum of the homogeneous solutions described in previous sections and the particular solutions described here.

#### 5.1.3.1.2.2. CRAM

The solver kernel based on the Chebyshev Rational Approximation Method (CRAM) is described in detail in references [ORIGEN-IA11, ORIGEN-Pus11, ORIGEN-Pus13, ORIGEN-PL10]. Compared to the MATREX solver, CRAM generally has similar runtimes but is more accurate and robust on a larger range of problems. CRAM relies on the lower upper (LU) decomposition, so the SuperLU library has been used. The accuracy of CRAM is related to the order, with an order 16 solution having a truncation error less than 0.01% for all nuclides in most problems.

Unlike many methods for solving this type of system of ODEs, the len>h of a step does not significantly affect the accuracy of CRAM. However, any significant errors from CRAM will shrink rapidly over multiple steps as long as there are no large changes in reaction rates. The CRAM solver has an efficient internal substepping algorithm that can perform multiple same-sized substeps (with the same transition matrix) very efficiently by reusing the LU decomposition. When using internal substepping, 2–4 substeps are typical, with a large gain in accuracy for marginal increase in runtime.

### 5.1.3.1.3. Power Calculation

The following is formula is used to calculate power during irradiation (\(\Phi > 0\)),

where \(\kappa_{\text{fi}}\) and \(\kappa_{\text{ci}}\) are
nuclide-dependent energy released per fission and “capture,” with
*capture* defined as removal minus fission:
\(\sigma_{\text{ci}} = \sigma_{i} - \sigma_{\text{fi}}\). The
\(\sigma_{\text{fi}}\) and \(\sigma_{\text{ci}}\) terms are
extracted from the transition matrix itself, whereas the
\(\kappa_{\text{fi}}\) and \(\kappa_{\text{ci}}\) are available
from a separate ORIGEN energy resource (see ORIGEN Data Resources
chapter). If the flux \(\phi\) is specified, then the power
can be calculated at any time according to Eq. (5.1.38).
However in reactor fuel systems, it is convenient to be able to specify the power
produced by the system and internally to the depletion code, to convert
the power to an equivalent flux. Solving Eq. (5.1.38)
for the flux, however,

it is apparent that a fixed power over a time step \(n\) does not lead to a fixed flux, due to changing isotopics that produce different amounts of power per fission and capture. ORIGEN performs a flux-correction calculation to obtain an estimate of the average flux over the step. The beginning-of-step flux is first calculated for the initial compositions: Eq. (5.1.39) is evaluated as \(\Phi(t_{n - 1})\), and then Eq. (5.1.38) is solved with that flux. The flux is then recalculated at the end of step \(\Phi(t_{n})\) using the estimated end-of-step isotopics, and the step-average flux \(\Phi_{n}\) is estimated as the simple average of the beginning and end-of-step fluxes, i.e. Eq. (5.1.40),

noting that the “predicted” flux at end-of-step \(\Phi^{\text{pred}}\left( t_{n + 1} \right)\) is based on “predicted” end-of-step isotopics, based on a beginning-of-step flux level.

### 5.1.3.1.4. Decay Emission Sources Calculation

ORIGEN can calculate the emission sources (and spectra) during decay for alpha, beta, neutron, and photon particles according to

where \(w_{i,x}\left( E \right)\) is the number of particles of type \(x\) emitted per disintegration of nuclide \(\text{i\ }\)at energy \(E\), using provided energy bins defined by energy bounds \(E^{g}\) to \(E^{g - 1},\) where \(g\) is an energy index. The fundamental data resources for performing emission source calculations are described in the ORIGEN Data Resources chapter.

#### 5.1.3.1.4.1. Neutron Sources

Computed neutron sources include neutrons spontaneous fission, \(\left(\alpha,n \right)\) reactions, and delayed \(\left( \beta^-,n \right)\) neutron emission,

with components that will be described below. The method of computing the spontaneous fission and delayed neutron source is independent of the medium containing the fuel. However, \(\left(\alpha,n \right)\) production varies significantly with the composition of the medium. The homogeneous medium \(\left(\alpha,n \right)\) calculation methodology has been adopted from the Los Alamos code SOURCES 4B [ORIGEN-Sho00, ORIGEN-WPC+99, ORIGEN-WPS+83].

The total yield of spontaneous fission neutrons from decay of nuclide
*i* is

where \(\frac{\lambda_{i,SFn}}{\lambda_{i}}\) is the fraction of decays which undergo spontaneous fission. The distribution of spontaneous fission neutrons, \(w_{i,SFn}\left( E \right)\) is given by a Watt fission spectrum,

where *Ai*, *Bi,* and *Ci* are model parameters.

The \(\left(\alpha,n \right)\) neutron source is strongly dependent on the low-Z content of the medium containing the alpha-emitting nuclides and requires modeling the slowing down of the alpha particles and the probability of neutron production as the \(\alpha\) particle slows down. The calculation assumes (1) a homogeneous mixture in which the alpha-emitting nuclides are uniformly intermixed with the target nuclides and (2) that the dimensions of the target are much larger than the range of the alpha particles. Thus, all alpha particles are stopped within the mixture. The yield of a particular \(\alpha\) is given by

where \(\frac{\lambda_{i,\alpha}}{\lambda_{i}}\) is the
relative probability of \(\alpha\)–decay, and \(f^\ell_{i,\ \alpha}\)
is the fraction of those \(\alpha\)–decays producing an \(\alpha\) particle with initial
energy \(E^\ell_{i,\alpha}\), and is considered fundamental
data. The *total* neutron yield from an alpha particle
\(\ell\) emitted by nuclide *i* and interacting with target
*k* is given by the following,

where \(S\left( E_{\alpha} \right)\) is the total stopping power of
the medium, \(\sigma_{k,\left( \alpha,n \right)}(E_{\alpha})\) is
the \(\left(\alpha,n \right)\) reaction cross section for target nuclide
*k*, and \(\frac{N_{k}}{N}\) is the fraction of atoms in the medium
composed of nuclide *k*. This expression is used to calculate the neutron yield
for each target nuclide and from each discrete-energy alpha particle
emitted by all alpha-emitting nuclides in the material. The stopping
power for compounds, rather than pure elements, is approximated using
the Bragg-Kleeman additivity rule. The energy-dependent elemental
stopping cross sections are determined as parametric fits to evaluated
data. Eq. (5.1.46) is solved for the total
neutron yields from the alpha particle, as it slows down in the medium by
subdividing the maximum energy \(E^\ell_{i,\alpha}\) into a number of
discrete energy bins and evaluating stopping power and
\(\left(\alpha,n \right)\) reaction cross section at the midpoint energy of
the bin. The distribution of \(\left(\alpha,n \right)\) neutrons as required
by Eq. (5.1.41) is

with the distribution of \(\left(\alpha,n \right)\) neutrons in energy, \(X^\ell_{i,k,\left( \alpha,n \right)}\left( E \right)\), calculated using nuclear reaction kinematics, assuming that the \(\left(\alpha,n \right)\) reaction emits neutrons with an isotropic angular distribution in the center-of-mass system. The maximum and minimum permissible energies of the emitted neutron are determined by applying mass, momentum, and energy balance for each product’s nuclide energy level. The product nuclide levels, the product level branching data, the \(\left(\alpha,n \right)\) reaction Q values, the excitation energy of each product nuclide level, and the branching fraction of \(\left(\alpha,n \right)\) reactions result in the production of product levels. A more detailed discussion of the theory and derivation of the kinematic equations can be found in [ORIGEN-WPC+99].

Delayed neutrons are emitted by decay of short-lived fission products.
The observed delay is due to the decay of the precursor nuclide. The
total yield of delayed neutrons from decay of nuclide *i* is

where \(\frac{\lambda_{i,Dn}}{\lambda_{i}}\) is the fraction of
decays which emit delayed neutrons. The delayed neutrons emitted per
decay of nuclide *i* at energy *E* is given by

where the spectrum \(X_{i,Dn}\left( E \right)\) is fundamental library data. Delayed neutrons are not important in typical spent fuel applications due to the very short half-lives of the parent nuclides, dropping off significantly after ~10 seconds, but they may be of value in specialized applications where calculating time-dependent delayed neutron source spectra is important.

#### 5.1.3.1.4.2. Alpha Sources

An \(\alpha\) slowing down calculation is performed as part of the \(\left(\alpha,n \right)\) neutron calculation. However, the alpha source (i.e. without considering slowing down in the media) is also available, simply as the sum of delta functions at the discrete initial alpha particle energies \(w_{i,\alpha}\left( E \right) = \sum_{\ell \in i} {Y^{\ell}_{i,\alpha} \delta \left(E - E^{\ell}_{i,\alpha} \right)}\) with yields \(Y^{\ell}_{i,\alpha}\), as required by Eq. (5.1.41).

#### 5.1.3.1.4.3. Beta Sources

The beta source (i.e. without considering slowing down in the media) is
available as the sum of the continuous emission spectra for each
\(\beta^{-}\) decay in nuclide *i*. The total yield of beta
particles from decay of nuclide *i* is

where \(\frac{\lambda_{i,\beta}}{\lambda_{i}}\) is the fraction of
decays which emit betas. The betas emitted by nuclide *i* at energy *E*
is given by

where the spectrum \(X_{i,\beta}(E)\) is fundamental data, independent of the media. The spectrum includes betas from allowed transitions and first, second, and third forbidden transitions.

#### 5.1.3.1.4.4. Photon Sources

The total yield of photons from decay of nuclide *i* is

where \(\frac{\lambda_{i,\gamma}}{\lambda_{i}}\) is the fraction of
decays which emit photons. The photons emitted by nuclide *i* at energy
*E* is given by

where the spectrum \(X_{i,\gamma}(E)\) is fundamental data and includes both line data from x-rays, gamma-rays and continuum data from Bremsstrahlung, spontaneous fission gamma rays, and gamma rays accompanying \(\left(\alpha,n \right)\) reactions. The Bremsstrahlung component of the photon source has been tabulated for various media and no on-the-fly slowing down calculation is performed.

### 5.1.3.1.5. Library Interpolation

Accurate solution of fuel depletion with Eq. (5.1.1) requires
coupling to self-shielding and neutron transport to accurately capture the
time-dependent change in space and energy flux distribution and 1-group
cross sections with isotopic change. This is in generally a fairly
computationally intensive problem compared to stand-alone depletion. In
typical assembly design and analysis, the same basic assembly problem
must be solved repeatedly with variations in power history, different
periods of decay/burnup, different moderator density, etc. A question
naturally arises: could the isotopics from numerous well-constructed
cases be saved and interpolated to the actual system? Interpolating the
isotopics themselves is fraught with difficulty. For example, consider
two cases with the same burnup but different periods of decay between
cycles. A better approach—the ORIGEN Automated Rapid Processing
(ARP)—was developed with the key realization that one can
reconstruct very accurate isotopics from stand-alone depletion
calculations by interpolating *transition matrices* rather than
*isotopics*.

The accuracy of the interpolation methodology compared to the coupled
transport/depletion solution (e.g., with TRITON) depends on the
suitability of the interpolation parameters and the deviation of the
desired system from the systems used to create the library. For example,
for thermal systems with uranium-based fuels, it was found that
enrichment, water density, and burnup were the dominant independent
variables and thus were best suited for interpolation. An example of the
variation of removal cross sections for key actinides is shown in
Fig. 5.1.2 for a Westinghouse 17 \(\times\) 17 pressurized water reactor
(PWR) assembly type with 5% initial enrichment in ^{235}U. Each cross
section has been divided by its initial value at zero burnup to show the
variation more clearly. ^{240}Pu has been observed to have the
most variation with spectral changes, with ~60% reduction in cross
section from beginning to end of life. The variations in ^{240}Pu
with respect to enrichment and moderator density are shown in
Fig. 5.1.3, Fig. 5.1.4, Fig. 5.1.5,
and Fig. 5.1.6.

Currently there are two interpolation methods: a Lagrangian based on low-order polynomials and a cubic spline with an optional monotonicity fix-up.

#### 5.1.3.1.5.1. Lagrangian Interpolation

Lagrangian interpolation [ORIGEN-REFed68] seeks the unique *n-1* order polynomial
that will pass through *n*-points of the function and then interpolating
to the desired point by evaluating the polynomial,

where \(x_{i}\) and \(y_{i}\) are the known x- and y-values in
the neighborhood of the desired x-value *x*, with *n* the number of data
points/order of Lagrangian interpolation. Note that Lagrangian
interpolation is by definition *local*, involving only points in the
neighborhood of the desired value. Global alternatives such as Hermite
cubic splines use the entire data set to construct the interpolants.
Common interpolation methods based on polynomials can have difficulty
with data that vary quickly and have uneven *x-*spacing, as is
expected with transition data. Polynomials tend to produce unphysical
oscillations in these cases. In cases with very small *y-*values
(~10^{-10}), oscillations of the interpolant can produce negative
interpolated values.

#### 5.1.3.1.5.2. Cubic Spline Interpolation

Cubic spline interpolation has been observed to produce fewer, lower frequency oscillations. Oscillations can be effectively eliminated by enforcing monotonicity on the interpolation: that is, additional max maxima or minima are not introduced by the interpolant between known values of the function. Monotonic cubic splines [ORIGEN-WA02] have shown particularly stable behavior and have been implemented as an interpolation option.