Depletion

Burnup matrices, the Bateman equation, and the CRAM solver used by Nucleide.

Nucleide solves nuclide transmutation over a single time step using a burnup matrix formulation and the Chebyshev Rational Approximation Method (CRAM) for the matrix exponential.

Overview

For a vector of nuclide atom densities $\mathbf{n}(t)$, the linearized depletion equation is

where $\mathbf{A}$ is the burnup matrix. Each column of $\mathbf{A}$ describes how a parent nuclide is lost and what it produces; each row describes how a daughter nuclide is gained.

Building the burnup matrix

Nucleide builds $\mathbf{A}$ from a Chain and a set of one-group reaction rates.

Decay

For nuclide $i$ with half-life $t_{1/2,i}$, the decay constant is

The diagonal entry loses atoms at rate $\lambda_i$:

Each decay mode with branching ratio $b_{i \to j}$ adds production to the off-diagonal:

A decay entry may omit target (or set target="nothing", the same convention as target-less reactions): the daughter falls outside the modeled chain, so the diagonal loss still applies but no gain term is built — pure loss. An explicit target naming a nuclide absent from the chain stays a loud error, never silent loss.

Reactions

For a one-group reaction rate $\sigma\phi$ of type $r$ on nuclide $i$:

• Loss: $A_{ii} \mathrel{-}= \sigma\phi$ once per reaction type, even if the chain lists multiple branched entries for that type.
• Production: each reaction entry with branching ratio $b$ and target $j$ contributes $$A_{ji} \;\mathrel{+}=\; \sigma\phi \, b.$$
• Light-particle secondaries: reactions such as $(n,\alpha)$, $(n,p)$, and $(n,d)$ emit light nuclides that are tracked when they appear in the chain (e.g. $(n,\alpha) \to \mathrm{He\text{-}4}$).
• Fission: if fission yields $y_{i \to k}$ are present, production is distributed over all fission products: $$A_{ki} \;\mathrel{+}=\; \sigma\phi \, y_{i \to k}.$$

The final matrix $\mathbf{A}$ has units of $\mathrm{s}^{-1}$ and is sparse, with one diagonal entry per chain nuclide even if it is stable.

CRAM: matrix exponential by partial fractions

CRAM evaluates the action of the matrix exponential $e^{\mathbf{A}t}$ on $\mathbf{n}_0$ in incomplete partial fraction (IPF) product form:

where $m = 8$ for order-16 and $m = 24$ for order-48. The coefficients $\alpha_0$, $\alpha_k$, and $\theta_k$ are pre-computed complex numbers published by Pusa. Each pole is applied sequentially to the running solution vector; these IPF coefficients are not interchangeable with the residues of the classic partial-fraction sum form.

Coefficients

Nucleide ships the CRAM-16 and CRAM-48 coefficient sets verbatim. CRAM-16 is cheaper; CRAM-48 is the production default.

The CRAM-16 scaling constant is $\alpha_0 = 2.124853710495224 \times 10^{-16}$.

Analytic Bateman fast path

For decay-only chains the burnup matrix is lower-triangular $\mathbf{\Lambda}$ and the depletion equation admits a closed form (Bateman, Proc. Cambridge Philos. Soc. 15 (1910), 423–427)^bateman-1910:

with $\lambda_j = 0$ for stable nuclides. $\mathbf{\Lambda}$ is assembled from the same decay data as $\mathbf{A}$ (spontaneous-fission daughters skipped, $\mathrm{He\text{-}4}$/$\mathrm{H\text{-}1}$ secondaries from alpha/proton decay included).

When the diagonal entries are distinct and the chain is acyclic, the eigendecomposition $\mathbf{\Lambda} = \mathbf{C}\,\mathbf{D}\,\mathbf{C}^{-1}$ with $\mathbf{D} = \mathrm{diag}(-\lambda)$ has the unit-diagonal recursion (Amaku–Pascholati–Vanin)^amaku-2010

Activity and cumulative-decay observables follow (E8) of crates/depletion/src/bateman.rs: $A_i = \lambda_i N_i$ (implemented by activity_vec in crates/depletion/src/integrate.rs) with the diagonal cumulative integral $\int_0^t e^{-\lambda s}\,ds = (1 - e^{-\lambda t})/\lambda$ ($\lambda = 0$ limit $\to t$; cumulative_decays in crates/depletion/src/inventory.rs). Labels (E6)–(E7) are intentionally unassigned in both the code and this page so the two stay consistent; a code-side close-the-gap renumber is left to the depletion owner (no formula changes either way).

$\mathbf{C}$/$\mathbf{C}^{-1}$ are time-independent, so Nucleide builds them once per chain (BatemanCache in crates/depletion/src/bateman.rs, generated at runtime from the crate’s own chain data) and reuses them across timesteps; each solve is two triangular matvecs. The single-chain Bateman sum is the special case of this form for a linear chain, and the two/three-nuclide $A \to B \to C$ formulas in the worked example above drop out directly.

Two precision variants share the cache: bateman (plain accumulation) and bateman_hp (f64-careful summation only — per-row terms sorted by magnitude and summed with Neumaier compensation, plus exp_m1 forms for small $\lambda t$; no extended-precision types, no new dependencies). On a synthetic 10-order-spread chain (year-scale parent, microsecond-scale progeny) the HP variant beats-or-matches standard summation while both conserve atoms to $10^{-8}$.

Algorithm flow

Worked example: sequential decay chain

Consider a three-nuclide chain:

with $\lambda_A = 10^{-6}\,\mathrm{s}^{-1}$, $\lambda_B = 10^{-5}\,\mathrm{s}^{-1}$, and $C$ stable. The burnup matrix is

Starting from $n_A(0) = 10^{15}$ atoms and $n_B(0) = n_C(0) = 0$, the analytic solution after $t = 10^5\,\mathrm{s}$ is

1. $n_A(t) = n_A(0)\,e^{-\lambda_A t} = 10^{15}\,e^{-0.1} \approx 9.0484 \times 10^{14}$.

2. $n_B(t) = n_A(0)\,\frac{\lambda_A}{\lambda_B - \lambda_A}\, (e^{-\lambda_A t} - e^{-\lambda_B t}) \approx 5.9662 \times 10^{13}$.

3. $n_C(t) = n_A(0) - n_A(t) - n_B(t) \approx 3.5501 \times 10^{13}$.

Nucleide’s CRAM-48 result matches this analytic solution to better than one part in $10^7$ for $n_B$ and one part in $10^9$ for $n_C$. This exact chain is exercised in crates/depletion/src/lib.rs as analytic_two_step_decay_cram48.

Time-series integrators

For schedules with per-step rates, crates/depletion/src/integrate.rs chains single-step CRAM solves. Each Step carries a length dt in seconds plus prescribed one-group rates in the same layout as DepletionSystem::build; the matrix is rebuilt from each step’s rates, so time-varying rates are handled exactly per step. Integrator selects the method: Predictor (one CRAM solve per step), Cecm (predictor plus a corrector re-solve from the step start with the averaged matrix), or Cf4 (fourth-order commutator-free quasi-Magnus, two half-step CRAM solves). Each step’s exponential additionally honors the Method kernel: Cram runs the CRAM order named above, while Bateman/BatemanHp evaluate the analytic closed form per step (steps with live rates fall back to CRAM-48; CF4 stages scale the timestep rather than the matrix so per-stage systems stay in decay form). Because rates here are prescribed rather than recomputed from the predicted composition, the corrector matrix equals the predictor matrix and Cecm reproduces Predictor up to roundoff; the two-solve structure is kept so flux-coupled extensions can plug a recomputed matrix into the same code path. Under the same piecewise-constant reconstruction each Cf4 stage is a half-step exponential, so its accuracy is at least as good as Predictor on the analytic test chain.

Results come back as TimeSeries with one row per time node in chain-nuclide order (times[0] == 0.0 echoes the initial state):

• times: cumulative time nodes in seconds.
• atoms: atom counts per node.
• activity: $A_i = \lambda_i N_i$ in Bq, clamped at zero.
• decay_heat: $H_i = \lambda_i N_i E_i$ in W, clamped at zero.

Per-nuclide decay energies resolve chain-first: the chain’s decay_energy attribute (eV converted to MeV) when positive and finite, else the ENDF/B-VII.1 prompt-decay-energy table (nucleide_nuclei::data::decay_energy_mev_by_name), else 0.0 (stable and synthetic nuclides) — one evaluated table, no parallel hand-kept fallback.

Year-length convention

Production code uses the Julian year: 1 yr = 365.25 d = 31,557,600 s (SECONDS_PER_YEAR in crates/depletion/src/inventory.rs; the y unit in crates/alara-io/src/schedule.rs and crates/alara-io/src/deck.rs). Rationale: the Julian year is the standard year-scale interval for time-unit conversion, and both the depletion solver and the ALARA schedule/deck layers convert through it.

Two deliberate neighbors differ and must not be “unified” by accident:

• crates/fispact-io/src/inventory.rs uses a 365-day year to match the ALARA reference tool’s file-format convention — format fidelity, not physics.
• Some test anchors divide by the Gregorian mean year 365.2425 d ($3.1556952 \times 10^7$ s/yr), e.g. the Cs-137 30.08 yr anchor in the nuclei half-life tests. The gap between conventions there is $\approx 6 \times 10^{-4}$ yr against a $0.01$ yr tolerance, so no validation gate moves either way.

Decision (no gate movement): keep the Julian production convention, document the split here, and migrate test anchors to 31,557,600 s/yr only where the current tolerance already hides the drift — tolerances themselves stay untouched. The alara-io and fispact-io year definitions are referenced from here rather than restated.

Assumptions and limitations

• One-group reaction rates: CRAM in Nucleide uses a single scalar rate per reaction channel. Multigroup fluxes must be collapsed externally.
• Fixed cross sections over the step: $\mathbf{A}$ is treated as constant during $t$; multi-step schedules rebuild it from each step’s prescribed rates, with no flux renormalization between steps.
• Fission yields are single-energy: the yield set at the lowest incident neutron energy is used, matching OpenMC’s default.
• Decay branching ratios: chains read from XML keep the file values verbatim; only programmatically built chains are renormalized (the largest branch is adjusted so the modes sum to exactly 1.0).
• Stable nuclides have zero diagonal: they act as sinks in the chain. The Bateman fast path handles them through limit forms ($e^0 = 1$).
• CRAM order: order-16 is adequate for short, simple chains; order-48 is recommended for production burnup matrices. Pure-decay chains may use the bateman/bateman_hp closed form instead (same method= selector; automatic CRAM-48 fallback outside decay-only triangular topology).
• Bateman conditioning: the closed form sums alternating-sign terms, so long chains (tens of members) or extreme half-life spreads can suffer catastrophic cancellation in float64 — bateman_hp (sorted Neumaier summation) extends the usable range, but ill-conditioned chains should use CRAM-48, which remains the default. The automatic fallback covers non-triangular topology, not ill-conditioning: pick the method deliberately for long chains.

Related work beyond Nucleide

Adjacent capabilities found in PyNE, OpenMC, and the literature that Nucleide does not implement:

• Multi-step time integration: OpenMC ships predictor-corrector integrators (CECM, CELI, CF4, EPCRK4, and their SI variants) that recompute reaction rates between steps; Nucleide’s Integrator::{Predictor, Cecm, Cf4} uses prescribed per-step rates, so Cecm currently degenerates to Predictor.
• Substepping: OpenMC’s CRAM solver can subdivide a step while reusing the LU factorizations; Nucleide has no substep parameter.
• Transport-coupled depletion: OpenMC normalizes tally reaction rates to absolute rates from a power or source normalization; Nucleide assumes rates are already given in $\mathrm{s}^{-1}$.
• Energy-dependent fission yields: OpenMC selects between constant, thermal, and thermal+fast yield modes; Nucleide uses the lowest-energy yield set.
• Feed and removal terms: OpenMC adds transfer rates and external source vectors to the burnup matrix for MSR-style continuous processing.
• Chain truncation: PyNE’s transmutagen-generated solver truncates the transmutation tree by density tolerance; Nucleide uses a fixed chain.
• Krylov subspace methods: an alternative to CRAM for the Bateman equations; see Josey et al. (2017)^josey-2017.
• Precomputed decay datasets: the radioactivedecay package ships precomputed decay matrices for fast lookup; Nucleide’s BatemanCache follows the same time-independent-cache idea but generates the eigendecomposition at runtime from its own chain data (clean-room textbook implementation, no vendored data).