EMPYREAN

Two axes of validation: the distribution channels agree with each other (every channel returns bit-identical numbers from one shared engine), and the engine agrees with external sources (JPL Horizons, ASSIST). The stack under test:

Every external comparison in this report is driven by a runner script in empyrean-validation. Each script is the full configuration we used to invoke that tool — settings, force model, weighting, rejection. If you want to reproduce a number in any panel, the script tells you exactly how we got there. (The tools themselves are external dependencies and are not linked into empyrean — they're run independently and their JSON outputs are merged into the report at the comparison step.)

How closely empyrean's propagated positions match an independent external reference. The scorecard is the at-a-glance agreement; the heatmap breaks it out per object and propagation offset (color scale tagged to encounter-distance thresholds). Choose the reference below.

Position error vs JPL Horizons over time, log-y. Median curve plus IQR band per population reduces the 44 individual lines to 13 population curves + named outliers. Toggle between rust-only (default) and 5-channel overlay; in All channels mode, c/cli/python lie exactly on rust (bit-exact, the lines overlap at chart resolution) and core is bit-identical to rust on the chaotic-divergence rows and on all but two objects, differing from rust by at most ~0.3 m on Didymos and 2026 FQ12 across all dt — a benign round-off-order difference in the core binary, not chaotic divergence. On this median-of-all-rows log-y plot all five channels overlap.

Log-y |empyrean − ASSIST| vs |ASSIST − Horizons| over time. ASSIST (Holman et al. 2023, REBOUND/IAS15/DE440) is the independent N-body reference. Triangular markers = comets and ISOs (g(r) implementation differs from ASSIST). Circles = self-perturbers (ASSIST treats as test particles). Threshold lines at 1 km / 100 km / 1 AU.

Per-row timing (log-y) for empyrean vs ASSIST, broken out by propagation_uncertainty mode so the comparison is apples-to-apples: empyrean's f64 path against single-particle ASSIST (f64), and empyrean's Jet1 + 6×6 covariance path against ASSIST with 6 first-order variational particles. The merge keys on (object, dt, propagation_uncertainty) so each empyrean row pairs with the matching ASSIST mode. Caveat: both ASSIST and empyrean f64 rows are fresh per-call integrations from the initial epoch to the target (the ASSIST runner constructs a new rebound.Simulation() and re-integrates from t₀ per dt), so the f64 panel is a head-to-head comparison. Where ASSIST's f64 column reads sub-microsecond, the cause is REBOUND/ASSIST's cached ephemeris-interpolation state and IAS15 step reuse across the closely spaced output epochs, not a different timing methodology.

empyrean's per-call cost across the STM-bearing surface in one frame. First-order pairs empyrean's Jet1 (6 partials) against ASSIST's 6 first-order variational particles — the one mode where a like-for-like external comparison exists, since ASSIST is first-order only. Second-order (empyrean's Jet2, 6 + 21 partials) and Auto (empyrean's adaptive mode) are shown as empyrean per-call cost on their own.

Sky-plane (RA/Dec) agreement between empyrean's predicted positions and an independent external reference, by population. Median and 95th-percentile separation, in milliarcseconds (1 mas ≈ the Gaia astrometric noise floor). The chart below shows how that separation grows with propagation offset.

Predicted-vs-Horizons RA/Dec separation per object and propagation offset (milliarcseconds); color scale anchored to the Gaia 1 mas floor.

RA·cos(δ) vs Dec residuals (mas) against the JPL Horizons reference. Clipped to ±100 mas to show the main cluster; rows beyond ±100 mas are listed in the outlier sidecar table below (these are typically pathological cases like 2020 CD3 propagated 10 yr before its 2020 mini-moon capture). 1σ (solid) and 3σ (dotted) ellipses overlaid from the cleaned cloud, assuming a zero-mean iid Gaussian model — for context, Gaia DR3 single-frame astrometric precision is ≲ 1 mas, ground-based survey typical residuals are ~100–300 mas. Rust channel.

dRA·cos(δ) (top) and dDec (bottom) residuals vs JPL Horizons over the test-epoch time window, symlog y so chaotic outliers and machine-epsilon-grade rows both render. Solid lines: median per population. Shaded band: 25-75 percentile (IQR) per population. Two stacked panels eliminate the solid/dashed-overlay readability problem. Rust channel.

For every OD test case, four diagnostics drive the "did the differential corrector land in the same minimum?" question across channels. Post-fit RMS, reduced χ² = χ²/ν, iteration count, and fitted-state drift to the core channel are all signals that an OD pipeline is not solving the same problem on the same data. Cross-channel OD agreement in this run is bit-identical in fitted state for rust / c / cli / python (see the per-object scatter below — every point sits on the 10⁻⁹ km chart floor). The remaining story is Empyrean vs find_orb on the same observations: find_orb fit 31 of the objects this run (the other 8 hit its 10 min/fixture cap and read “not run” below). Of the 27 that also have a converged Empyrean fit to pair against, the 23 with full-coverage find_orb solutions agree to within a factor of ~1.6 in post-fit RMS (median 1.09×, all within 2×) — no large regressions this run. The four exceptions (2008 TC3, 2023 CX1, 2024 BX1, 67P) are not clean comparisons: find_orb falls back to a degenerate low-coverage placeholder on these short / sparse arcs (≤ 2% of observations used, fitting only its 2 anchor points), so their Δ is flagged “rejection-mismatch” in grey rather than counted as a fit-quality difference. Empyrean fits the full arc on those same objects (e.g. 2008 TC3: 809 obs at 0.84″ RMS).

No orbit-determination rows in this report. Run the OD subset to populate this section.

For every object whose JPL SBDB reference carries a non-zero non-gravitational signal (the Yarkovsky-detected NEOs Apophis and Bennu, plus the comets), the runner re-fits the same optical arc with solve_for = StateAndNonGrav and emits a second OD row (test_type = non_grav_recovery) carrying the fitted Marsden A1/A2/A3 and their 1σ. The σ are read straight from the diagonal of the fitted 9×9 state+A1/A2/A3 covariance: σ_A1 = √C[6][6], σ_A2 = √C[7][7], σ_A3 = √C[8][8]. The pass criterion is a per-coefficient σ-consistency check against JPL: z = (od_a· − ic_a·) / od_a·_sigma, PASS when |z| ≤ 3 on every coefficient whose JPL reference is non-zero. A fitted value that comes back None is a loud FAIL, not a blank: it means the 9×9 covariance was absent (the fit silently fell back to a 6-parameter state-only solution) or the fitted coefficient was non-finite — exactly the regression this section exists to catch. All four distribution channels (rust / c / cli / python) are shown side-by-side so a per-channel FFI marshaling drop of the fitted non-grav block is immediately visible.

No non-gravitational-recovery rows in this report. Run the OD subset against objects with a known SBDB non-grav signal to populate this section.

Cross-channel agreement broken out by test type. The five channels run end-to-end through libempyrean (or empyrean-core directly for the core channel). Bit-identical results are the expected outcome for propagation and ephemeris; OD is bit-identical in fitted state (see §9 fitted-state-drift chart).

On the denominators: rust's 2588 / 2760 / 35 totals cover the full uncertainty grid (Auto + first_order_with_cov + second_order_with_cov + f64_no_cov). Non-rust channels run only the modes their public API exposes (first_order_with_cov + f64_no_cov for c / cli / python; the same two for core, replayed twice through different entry points), so their compared denominators are ≈ rust ÷ 2 by design — not a coverage gap. The fidelity claim is per-row bit-exactness against the matching core row; the denominator difference is a sampling artefact of which modes are end-to-end production-shipped.

cli·propagation 77/77  ·  cli·ephemeris 60/60  ·  python·propagation 77/77  ·  python·ephemeris 60/60  ·  rust·propagation 154/154  ·  rust·ephemeris 120/120

Empyrean is uncertainty-first by design. When you call empyrean.propagate(orbit, t) with an orbit that carries a covariance, the propagator dispatches automatically to first-order Jet1 / STM integration — no flag, no opt-in. The covariance is propagated through the State Transition Matrix and a 6×6 covariance falls out of the same call. When the input orbit has no covariance, the propagator dispatches to plain f64 state-only integration and you pay nothing for what you don't use. The validation suite exercises both modes side-by-side, tagging each propagation row with propagation_uncertainty = "first_order_with_cov" or "f64_no_cov". This section shows the cost of carrying uncertainty per population — i.e., the price of the production hot path users hit every day.

No paired Jet1/f64 rows in this report. Re-run rust + core channels to populate this section.

For every OD test case where an SBDB solution (or find_orb fit) is available, Empyrean's fit is compared to that reference in Keplerian element space, bidirectionally — once at Empyrean's epoch (reference propagated via STM) and once at the reference's epoch (Empyrean propagated). For linear propagation, Mahalanobis is direction-invariant so the bidirectional pair acts as a sanity check on the STM chain. Three pathology fingerprints:

• CONSISTENT (σ_equiv < 1, green): the two ellipsoids agree within their joint uncertainty.
• CORRELATION (joint d² >> marginal d²): off-diagonal correlation in the joint covariance is driving the discrepancy, not any single element. Canonical example: short-arc fit with a/e degeneracy.
• ROTATION (principal-axis rotation angle > 30° with similar eigenvalue spectra): the two ellipsoids have different orientations. Often artifact of propagation across a non-linear region.
• (unlabeled high σ_equiv): large σ_equiv that triggers neither CORRELATION nor ROTATION typically means a scale-only disagreement (missing perturber, weight scheme mismatch) or a true difference in the fit. Canonical example: 67P (σ_equiv ≈ 67.7 in both directions, near-zero rotation and joint ≈ marginal). Asymmetric bidirectional pairs (e.g., 2020 CD3 σ_equiv = 106.82 one way / 89.71 the other, both with principal-axis rotation flagging it ROTATION) reflect non-linear regions like chaotic-capture geometries — the STM Mahalanobis is only direction-invariant in the linear limit.

The Mahalanobis decomposition column (d²_combined / d²_marginal) flags whether the joint metric is or is not marginally explained.

No orbit-comparison rows in this report (no SBDB or find_orb sidecar found).

Provenance for this run — the frame, ephemeris, force model, and external references behind the numbers in this report.

References
· Holman, M. et al. 2023, "ASSIST: An ephemeris-quality test-particle integrator", PSJ 4(4), 69 (DOI 10.3847/PSJ/acc9a9).
· Vereš, P. et al. 2017, "Statistical analysis of astrometric errors for the most productive asteroid surveys", Icarus 296, 139.
· Eggl, S., Farnocchia, D., Chamberlin, A. B., Chesley, S. R. 2020, "Star catalog position and proper motion corrections in asteroid astrometry II: the Gaia era", Icarus 339, 113596.
· Park, R. S. et al. 2021, "The JPL Planetary and Lunar Ephemerides DE440 and DE441", AJ 161, 105.
· Marsden, B. G., Sekanina, Z., Yeomans, D. K. 1973, "Comets and Nongravitational Forces. V", AJ 78, 211.
· Vokrouhlický, D. et al. 2015, "The Yarkovsky and YORP Effects", in Asteroids IV, p. 509.
· Everhart, E. 1985, in "Dynamics of Comets" (Reidel), p. 185.