Fusion Neutron Sources

D-D/D-T fusion neutron spectra, the Miller tokamak source model, and MCNP/Serpent card emission.

Transport codes need to know where neutrons are born and with what energy. This page pins the physics behind nucleide.plasma_source, the tokamak fusion-neutron source capability: temperature-broadened D-D (2.45 MeV) and D-T (14.1 MeV) spectra, a parametric Miller-geometry plasma driven by caller-supplied L/H/A-mode profiles, seeded sampling to particle vectors, and emission of MCNP source-definition cards (SDEF) and Serpent source cards (src) with an emission report that quantifies what the cards preserve and what they lose.

The spectrum of thermonuclear neutrons is Gaussian to good approximation, with a width proportional to $\sqrt{T_i}$^brysk-1973. Nucleide implements the closed-form coefficient fits of Ballabio et al.^{ballabio-1998}, which refine the Brysk scaling with a small mean-energy shift and a weakly temperature-dependent width; at $T_i = 0$ the spectrum collapses to the monoenergetic nominal line. The spatial side is the public ITER/EU-DEMO parametrization of Fausser et al. — a Miller flux-surface map with L/H/A-mode density and temperature profiles^fausser-2012 — and the neutron emission is reactivity-weighted with the Bosch–Hale fits in the Atzeni–Meyer-ter-Vehn parametrization^{bosch-hale-1992}. All equations and coefficients are transcribed as published facts and pinned by independently recomputed reference values; no simulator source was read to produce them.

Overview

nucleide.plasma_source builds fusion neutron sources at three levels of fidelity. The simple tier is a point source or an axisymmetric ring on the full torus (a tokamak midplane circle): positions are exact, and the only physics is the energy spectrum. The parametric tier models the plasma cross-section with Miller flux surfaces and caller-supplied L/H/A-mode profiles of ion density and temperature; birth positions weight the local neutron production, so the source follows the profiles and the geometry rather than a hand-drawn shape. The lattice tier drops axisymmetry itself for non-axisymmetric (stellarator-class) plasmas: the caller supplies the full 3D birth distribution as an explicit point cloud, and the same fuel mixtures, spectra, and card emission apply per node. All tiers sample seeded, deterministic particle vectors (position, direction, energy, weight) and emit SDEF and src cards with an emission report.

Units follow the transport-code card convention: lengths in centimetres, energies in MeV, ion temperature in keV. Profile densities are caller data in m⁻³ (the Fausser convention) and converted internally; reactivity is reported in m³/s.

Mathematical details

S1 — Temperature-broadened neutron spectrum

The spectrum of thermonuclear neutrons is Gaussian to good approximation, with a width proportional to $\sqrt{T_i}$^brysk-1973. Nucleide uses the closed-form fits of Ballabio et al., Table III^{ballabio-1998}, which refine the Brysk scaling with a small mean-energy shift and a weakly temperature-dependent width:

Here $T_i$ is in keV and the mean shift is computed in keV and added to $E_0$ in MeV; $(a_1, \dots, a_4)$ and $(w_0, b_1, \dots, b_4)$ are the published Table III coefficients for each reaction. $E_0$ is the nominal cold line: 14.021 MeV for D-T and 2.4495 MeV for D-D — the “14.1 MeV” and “2.45 MeV” lines of the literature. At $T_i = 0$ both corrections vanish and the spectrum is the monoenergetic nominal line.

S2 — Thermonuclear reactivity

Emission weights the profiles by the thermonuclear reactivity $\langle\sigma v\rangle$ of Bosch & Hale^{bosch-hale-1992}, evaluated in the compact parametrization reproduced by Atzeni & Meyer-ter-Vehn:

$T$ is in keV and num, den are per-reaction rational correction polynomials:

The reactivity vanishes as $T \to 0$ (the exponential beats the power-law growth), which the parametric sampler relies on at cold separatrix temperatures; the crate returns exactly 0 at $T = 0$. The D-D row is the published D-D fit — the neutron-relevant D(d,n)³He branch shares it to the precision of the coefficients. Where the fit leaves its validity domain (the D-D $\eta$ goes non-positive on roughly 965–2720 keV, far above the published range) the crate raises a clear out-of-domain error instead of returning NaN.

S3 — Miller flux-surface map and volume element

The parametric plasma cross-section is the Miller map of Fausser et al.^fausser-2012, in the spelling common to that paper and the MIT openmc-plasma-source package:

Major radius $R_0$, minor-radius label $a \in [0, a_\text{minor}]$, elongation $\kappa$, triangularity $\delta$, and Shafranov factor $\mathrm{esh}$ are the shape parameters (lengths in cm). Neutron production lives on flux surfaces, so the volume element matters:

— the poloidal Jacobian times the toroidal radius $R$. The $R$ factor is mandatory: without it the core is over-represented. Closed forms pinned by the crate: in the torus limit ($\kappa = 1$, $\delta = 0$, $\mathrm{esh} = 0$) $J = a$ exactly; at $\delta = 0$ (any $\kappa$, $\mathrm{esh}$)

the cross-section area is $\pi\kappa a^2$, and the torus volume is $2\pi^2\kappa a^2 R_0$ for any $\mathrm{esh}$ — the Shafranov terms integrate out, so the shift moves surfaces outboard without changing the volume. (The $|J| = J$ step needs $J > 0$ everywhere, which the crate’s $|\mathrm{esh}| < a_\text{minor}/2$ validation guarantees.) The map is up-down symmetric and the vertical extent is $\max|Z| = \kappa\,a_\text{minor}$.

S4 — L/H/A-mode plasma profiles

Profiles over the minor radius are the L/H/A-mode parametrizations of Fausser et al.^fausser-2012. L mode is a single power-law falloff:

H and A mode (same algebra; the modes differ by parameter choice, not formula) add a pedestal at $r_\text{ped}$ and a linear edge:

the temperature core instead carries a $\beta_T$ exponent, $(1 - (r/r_\text{ped})^{\beta_T})^{\alpha_T}$. Profiles are inputs: the caller supplies the centre, peaking, pedestal, and separatrix values and the mode label; nothing here computes or fits profiles, and no equilibrium is solved. Density is caller data in m⁻³ and temperature in keV.

S5 — Reactivity-weighted emission strength

Neutron production at minor radius $r$ follows the local ion density and temperature through the reaction rate:

$f_\text{fuel} = 1/4$ for equimolar D-T ($n_D = n_T = n_i/2$) and $1/2$ for pure D-D — the ½ avoids double-counting identical reactant pairs and follows the Fausser/openmc-plasma-source strength convention. A D/T fuel mixture with fractions $f_D + f_T = 1$ instead weights both branches of the rate:

the standard $n_a n_b \langle\sigma v\rangle$ rate with the $1/(1+\delta_{ab})$ same-species guard (the DRESS rate construction of Eriksson et al., and the form upstream openmc-plasma-source implements). At $f_D = f_T = 1/2$ the first term reproduces the equimolar factor $1/4$ exactly, and at $f_D = 1$ the second term reproduces the pure D-D factor $1/2$; birth energies draw from the local mixture of the two Ballabio lines, branch by branch. Birth positions weight $S \cdot R \cdot |J|$, the volume element of (S3-V); the total relative strength is its toroidal integral, $S_\text{tot} = 2\pi \int_0^a\!\int_0^{2\pi} S\,R\,|J|\,d\theta\,dr$ — arbitrary global scale, but ratios between configurations are meaningful. Because $\langle\sigma v\rangle \to 0$ as $T_i \to 0$, a cold separatrix makes no neutrons.

S5-P — Per-species temperatures without a fuel mixture

A single-fuel config carrying a species pair $(T_D, T_T)$ reacts at the landed mixture pair temperatures — single-fuel D-T at the mass-weighted $T_{DT} = T_D + \tfrac{2}{5}(T_T - T_D)$, single-fuel D-D at $T_D$ — with the Ballabio lines following per branch (D-T at $T_{DT}$, D-D at $T_D$). No new temperature concept enters; the pair composes with tail and sector exactly as in mixtures^{eriksson-2016}, and the profile ion temperature is then unused for rate and spectrum (still validated). An equal pair $(T, T)$ reproduces the shared-temperature kernel bit-for-bit. Hand vector: $(T_D, T_T) = (20, 30)$ keV gives $T_{DT} = 24$ keV with $\langle\sigma v\rangle_{DT} = 5.414667327922193\times 10^{-22}$ m³/s and the D-T line at mean $14.077934372230146$ MeV, sigma $0.3700390508551786$ MeV.

S5-T — Deuterium hot-tail fraction

Auxiliary-heated plasmas carry a non-thermal tail that a pure Maxwellian mis-models. The one pinned non-Maxwellian shape^{eriksson-2016} splits the deuterium population into a bulk $(1 - \eta)$ at $T_D$ and a hot tail $\eta$ at $T_\mathrm{tail}$ (uniform scalars, flat in minor radius), while tritium stays Maxwellian at $T_T$. Each sub-pair reacts at its own mass-weighted relative temperature:

and the strength rule expands into the bulk/tail sub-rates:

The D-T arm carries two spectrum lines (at $T_{DT}$ and $T_{DTt}$) and the D-D arm three (at $T_D$, $T_\mathrm{mix}$, $T_\mathrm{tail}$), weighted by the same sub-rate rule, so strengths, sampled energies, card marginals, and the axis summary stay consistent. Each line is the Ballabio Gaussian at its effective temperature — an approximation in the same spirit as the species-pair lines, since the published fits assumed single-temperature Maxwellians. The pinned hand vector (flat profile, $n = 10^{20}$ m⁻³, $f_D = 0.7$, $f_T = 0.3$, $T_D = 20$ keV, $T_T = 30$ keV, $\eta = 0.05$, $T_\mathrm{tail} = 60$ keV) gates the rate at $10^{-12}$ relative, and $\eta = 0$ reproduces the no-tail mixture kernel exactly.

S6 — Arbitrary-3D birth-rate lattice

The lattice tier drops the last tokamak assumption, axisymmetry. The caller supplies the full 3D birth distribution as a point list — positions, relative birth-rate weights, and per-point ion temperatures — and the same fuel mixtures, Ballabio lines, and card emission apply per node. An optional field-period symmetry declares the points to be one base-sector cloud replicated $n$ times around the machine axis (toroidal angles identified modulo $2\pi/n$), with the symmetric total scaled by $n$; $n = 1$ is the identity. Anchors pinned by tests: a two-point hand vector (rates 2 and 1, total 3, mean birth position $(100, 0, 25)$, mono 14.021 MeV line), single-point reproduction of the point source, ring-cloud convergence to the analytic ring moments, and bit-for-bit reproduction of the expanded full-cloud stream by the folded base-sector stream. Cards carry the cloud’s cylindrical-$R$/vertical/energy marginals with a lattice row in the emission report, exactly like the parametric product-form cards.

S7 — ECRH accessibility

Heating access is a separate closed-form kernel in the same crate: given a wave frequency plus the caller’s magnetic-field and density profiles tabulated along a beam path, it reports where the wave meets the cyclotron resonance and where it meets a cut-off — no ray tracing, no launcher design. The cold electron gyrofrequency (Hutchinson, Principles of Plasma Diagnostics, 2nd ed. (2002), §7.2) is

and fundamental (harmonic $n$) absorption sits where $f = n\,f_{ce}$, i.e.

the cold resonance condition pinned by the ECRH review of Bornatici et al.^{bornatici-1983}. At finite electron temperature the resonance carries the Lorentz factor, $f = n\,f_{ce}/\gamma$, so the layer moves out to

the weakly-relativistic Maxwell–Jüttner downshift (Bornatici et al.^{bornatici-1983}). The O-mode cut-off ($P = 0$ of the cold Appleton–Hartree dispersion; Stix, Waves in Plasmas (1992), Ch. 1) sits at the plasma frequency,

and the X-mode right-hand cut-off inverts to

(Stix 1992, Ch. 1 — the low-field-side launch cut-off). Where $f \le f_{ce}$ the X1 wave is evanescent and the kernel returns no cut-off, never a number. Beamline locators take caller profiles at strictly increasing positions $s$ [m] and return crossing positions [m] by linear interpolation between bracketing nodes; exact nodal touches are reported. Scalar inputs are GHz / T / m⁻³ / keV; the per-port absorbed power itself lives downstream (blanket bookkeeping consumes these positions as plain floats).

Pinned hand vectors (170 GHz, ITER-class, $10^{-9}$ relative): cold resonance $B_{res} = 6.07306$ T at $n = 1$ ($3.03653$ T at $n = 2$); $\gamma(20\ \mathrm{keV}) = 1.03914$ shifts it to $6.31075$ T; O1 cut-off $3.58488\times 10^{20}$ m⁻³; X1 cut-off at 5 T $6.33418\times 10^{19}$ m⁻³ — while 100 GHz at 5 T ($f_{ce} = 139.96$ GHz) is evanescent. Beamline vectors: a 7→5 T fall over $s = 0..2$ m crosses the cold resonance at $0.92694$ m; a $10^{20}$→$5\times 10^{20}$ m⁻³ rise crosses the 140 GHz O1 cut-off at $0.71563$ m.

Algorithm

Ring and point sampling. A point source sits at its configured position; a ring source draws a uniform azimuth on the midplane circle, so radius and height are exact for every particle and only the azimuth carries sampling error. Energies are the monoenergetic line or a Box–Muller Gaussian draw ($\bar{E} + \sigma z$); directions are isotropic, with $\cos\theta$ uniform in $[-1, 1]$ and azimuth uniform in $[0, 2\pi)$.

Parametric sampling. At construction the sampler tabulates the strength-weighted radial distribution on a fixed 256-cell grid over $[0, a_\text{minor}]$; each cell precomputes its poloidal cumulative distribution on 256 points (trapezoid quadrature of the periodic $R \cdot |J|$ integrand), so drawing a particle needs no per-particle quadrature: $r$ from the radial table, $\theta$ from the owning cell’s poloidal table, $\varphi$ uniform, then $(R, Z)$ from the Miller map. The birth energy comes from the Ballabio Gaussian at the local ion temperature $T_i(r)$. Both samplers run on a small self-contained RNG (SplitMix64 seed expansion feeding xoshiro256**), so a pinned seed reproduces the stream on a given platform.

Tail sampling. With a deuterium hot tail set, the per-particle branch choice runs over the five bulk/tail sub-lines weighted by the sub-rate rule (S5-T) — one extra draw per particle on the mixture path, none when $\eta = 0$, so the no-tail stream is preserved exactly.

Lattice sampling. A lattice particle picks its birth node from the rate-weighted discrete distribution, then (with symmetry set) its field-period copy from the within-bin fraction of the same single draw — no extra draws on any path, so the folded stream reproduces the expanded stream exactly. Energies follow the node’s Ballabio line (or its mixture branch choice); directions are isotropic as above.

Source cards and the emission report. A monoenergetic spectrum renders as a single line with probability 1. A Gaussian renders as discrete lines at bin centers over a symmetric $\bar{E} \pm 4\sigma$ window (21 lines by default), with probabilities from the exact normal CDF (erf via the Abramowitz–Stegun 7.1.26 rational approximation, $|\epsilon| \le 1.5\mathrm{e}{-7}$, far below the six-significant-digit card format); the probability mass outside the window is dropped and reported. Parametric cards can only carry the marginals — radial and vertical birth histograms plus the global marginal energy spectrum — so the emission report also carries the joint-correlation distance $\tfrac{1}{2}\sum_{ij}|p_\text{true}(i,j) - p_r(i)\,p_z(j)|$, half the L1 distance between the true $(r, z)$ birth joint and the product of the marginals (0 for independent). Every report row records the fraction of the source strength the card accounts for and the fraction lost (rel_drift); reparsed marks rows machine-verified by feeding the emitted text back through that code’s own reader — true for SDEF, which must round-trip byte-identically through nucleide-mcnp-io (a failed round trip raises a clear error rather than shipping unverified text), and false by design for Serpent src, where the workspace has no source reader and the rows are checked analytically against the public input-manual spelling.

Worked example: ITER-like D-T ring

An ITER-like D-T ring at $R = 300$ cm, $z = 25$ cm with $T_i = 20$ keV. First the spectrum and reactivity closed forms through spectrum_moments (the Ballabio moments of a reaction at an ion temperature) and reactivity (the Bosch–Hale $\langle\sigma v\rangle$):

import nucleide.plasma_source as ps

ps.spectrum_moments("dt", 20.0)

#  'mean_mev': 14.072874253353246, 'sigma_mev': 0.33772176292161377}
ps.spectrum_moments("dd", 0.0)["sigma_mev"]  # 0.0 — cold spectrum is monoenergetic
ps.reactivity("dt", 20.0)                    # 4.330220397689176e-22 m^3/s
ps.reactivity("dd", 20.0)                    # 2.6025829587215236e-24 m^3/s

Evaluating (S1a)–(S1b) at $T_i = 20$ keV with the D-T Table III coefficients gives mean 14.0729 MeV and $\sigma$ 0.33772 MeV against the 14.021 MeV nominal line; at $T_i = 0$ the width vanishes and both reactions collapse to their monoenergetic lines (14.021 MeV and 2.4495 MeV). The reactivity ratio at 20 keV, $\langle\sigma v\rangle_{DT} / \langle\sigma v\rangle_{DD} \approx 166$, is why D-T dominates neutron production at fusion temperatures.

Then sampling and card emission through particles (a seeded particle vector) and emit_source_cards (both cards plus the emission report):

import numpy as np
import nucleide.plasma_source as ps

spec = {"kind": "ring", "radius": 300.0, "height": 25.0,
        "reaction": "dt", "ion_temperature_kev": 20.0}
out = ps.particles(spec, 4, seed=2026)
np.hypot(out["x"], out["y"])  # array([300., 300., 300., 300.]) — exactly R
out["z"]                      # array([25., 25., 25., 25.])
out["energy"]                 # array([14.17846035, 13.81800956, 14.3161297 , 14.02004292])

cards = ps.emit_source_cards(spec)
cards["sdef"]["card"].splitlines()[:4]
# ['SDEF POS=0 0 25', '     AXS=0 0 1', '     RAD=D1', '     ERG=D2']
cards["sdef"]["drift"][0]["rel_drift"]  # 6.337207213114038e-05

import nucleide.mcnp as mcnp
mcnp.parse_sdef(cards["sdef"]["card"])["card"] == cards["sdef"]["card"]  # True

Hand checks: every birth sits exactly on the ring (radius 300 cm to machine precision, height exactly 25 cm), so the azimuth is the only sampled spatial coordinate. The first particle’s energy, 14.1785 MeV, sits $0.31\sigma$ above the 14.0729 MeV mean. Re-running with the same seed replays the identical stream — the vectors above match the crate’s pinned Rust reference stream. The card keeps the ring exactly (RAD=D1 with SI1 L 300 300) and tabulates the Gaussian at 21 bin centers spanning $\bar{E} \pm 4\sigma$ (12.7863–15.3594 MeV). The emission report’s first row accounts for 0.999937 of the probability mass; the lost $6.34\mathrm{e}{-5}$ is the Gaussian tail beyond $\pm 4\sigma$, $1 - \mathrm{erf}(4/\sqrt{2}) = 6.334\mathrm{e}{-5}$. The SDEF text re-parses byte-identically through nucleide.mcnp.parse_sdef (the typed SDEF reader).

Assumptions and limitations

• Toroidal sectors. Partial-torus sources take start_angle/rotation_angle (birth angles uniform over the sector; strength scales by rotation/2π); full rotation reproduces the full-torus kernel exactly. Out-of-range angles raise a clear error — never a silent guess.
• Maxwellian fuel, several choices — plus one pinned tail. Reactants are Maxwellian: equimolar D-T (strength factor $1/4$), pure D-D ($1/2$), a D/T mixture {"D": f_D, "T": f_T} weighted by the two-branch rate (S5-M), or a mixture with per-species ion temperatures {"D": T_D, "T": T_T} (the D-T branch reacts at the mass-weighted effective temperature) — the equimolar, pure, and shared-temperature limits recover the simpler kernels exactly. A deuterium_tail fraction {"fraction": η, "temperature_kev": T_tail} adds the one pinned non-Maxwellian shape (S5-T: bulk plus hot tail with per-sub-pair lines); η = 0 recovers the no-tail kernel exactly, and the tail needs a fuel mixture. Fuel fractions that are negative, non-finite, or do not sum to 1 raise a clear error, as do negative temperatures and out-of-range tail parameters. General non-Maxwellian distributions beyond this shape stay out.
• Neutron branches plus proton bookkeeping. The sampler and the emitted cards are neutron-only (D-T + D-D lines; no T-T neutron branch — a pure-T mixture is a loud error). The D(d,p)T proton branch is reported as bookkeeping next to the source (the proton density is the D-D neutron-branch density, pinned 50/50); proton transport is out of scope.
• Profiles are caller inputs. Nothing computes or fits profiles and no equilibrium is solved. Malformed geometry or profile parameters (non-positive minor radius or elongation, triangularity outside $(-1, 1)$, pedestal radius outside $(0, a)$, non-positive centre density or temperature) raise a clear, named error, and a configuration whose total emission strength is zero raises a clear error instead of sampling nothing.
• Lattice clouds are caller inputs. Nothing computes the 3D birth distribution and no equilibrium is read: node positions, rates, and ion temperatures arrive as data (a point list, not a grid). Empty clouds, negative rates, zero totals, and non-symmetric fold inputs raise clear, named errors.
• Published fit domains. The Ballabio and Bosch–Hale fits are used within their published ranges. Negative or non-finite ion temperatures raise a clear error, and the D-D reactivity’s correction factor goes non-positive on roughly 965–2720 keV, where the crate raises a clear out-of-domain error rather than returning NaN.
• Cards carry marginals, not correlations. A product-form card cannot represent the correlated $(r, z)$ birth joint or the position–energy correlation; the emission report quantifies both losses (the truncated energy tail and the joint-correlation distance).
• Determinism per platform. A pinned seed reproduces the particle stream on a given platform; platform libm differences can shift the last ulp of individual streams, with no effect on moments or card text at any meaningful tolerance.
• No file writing in the crate. Particle vectors and card strings are returned to the caller; writing MCPL or MCNP files stays caller-side (nucleide.mcpl, nucleide.mcnp).
• ECRH reports positions, never designs launchers. Resonance and cut-off crossings are plain floats for downstream bookkeeping; ray tracing, beam optics, and absorbed-power computation stay out. An evanescent X1 reports no cut-off rather than a number, and empty, ragged, or non-monotonic beamlines raise a clear error.

Validation

• Rust unit tests pin the transcriptions: spectrum moments and reactivities at 1/10/20 keV recomputed independently from the published coefficients, the Miller closed forms (torus limit, the $\delta = 0$ area and volume identities, finite-difference agreement of the Jacobian at general $\delta$), profile definitions and pedestal continuity, sampler determinism, and a pinned four-particle stream for the RNG.
• validation/plasma_source_vs_openmc.py rebuilds every quantity from a second, independent in-script transcription of the published formulas. Always-on checks cover ring geometry and azimuth moments, Ballabio moments (< 1e-12 relative), pinned-seed determinism, the SDEF round trip, the tabulation tail drift against the $\pm 4\sigma$ tail mass, sampled flat-profile moments against the $\delta = 0$ closed forms ($\langle r^2\rangle = a^2/2$, $\langle z\rangle = 0$, and the outboard volume bias $\langle R\rangle = R_0 + a^2/(4R_0)$), the reactivity transcription (< 1e-12), and H-mode sampled birth moments against independent fine quadrature of $S\,R\,|J|$, plus the P11 species gates (pair-temperature branch fire, single-fuel pair arms at the landed pair temperatures, equal-T stream/card recovery) and the P12 lattice gates (two-point hand vectors, single-point recovery, ring-limit closure, symmetry-fold replication, card round trip, mixture branch fire) and the P13 hot-tail gates (sub-rate-rule moments and branch share, the drift tail-neutron share, zero-fraction recovery, loud errors).
• Container-only cross-checks against the MIT-licensed openmc-plasma-source package and NeSST (Ballabio helpers, sampled D-D and ring moments, the Miller map, L/H profiles, reactivities to ~5e-9 — the published coefficients’ own precision — end-to-end parametric moments, the lattice ring-limit probe, and the hot-tail blend probe). When the packages are unavailable these checks are reported as SKIP with the reason, never silently.
• ECRH hand vectors in ecrh.rs unit tests and tests/test_ecrh_access.py pin the closed forms at $10^{-9}$ relative (the 170 GHz resonance/shift/cut-off vectors above plus beamline crossings and loud-error vectors). ECRH has no oracle leg: closed forms need no cross-code leg, only the two independent transcriptions.
• tests/test_plasma_source.py replays the user-facing Python surface.

Related work beyond Nucleide

Adjacent capabilities Nucleide does not implement: the Eriksson et al. generalization beyond the one pinned single-tail-temperature shape — arbitrary non-Maxwellian reactant distributions, which covers beam-target spectra beyond the hot-tail approximation (the documented “not yet supported” boundary above)^{eriksson-2016}; adaptive Monte Carlo source iteration in the spirit of KDSource (GPL-licensed, not integrated); OpenMC’s arbitrary mesh- and CAD-based spatial sources — OpenMC’s fusion-specific sources themselves ship in the separate openmc-plasma-source package, which the validation harness cross-checks against at runtime; and synthetic diagnostic responses beyond the source description itself, the territory of NeSST.