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:

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}$.

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.

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$.
• 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.
• CRAM order: order-16 is adequate for short, simple chains; order-48 is recommended for production burnup matrices.

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 exposes a single constant-rate step.
• 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.