A Reproducible N-Body Pipeline and Numerical Convergence Framework for Retrograde Planet Nine Configurations

Mayone Maha Rajan (Architect & Curator) AI synthesis instrument: Google Antigravity (agentic model) Editorial and provenance-verification instrument: Claude (Anthropic), under the architect's direction Revision 3 (final) — June 2026

This paper was produced through human-directed AI synthesis. The human architect curated the inquiry, designed the integration and audit protocol, and is responsible for all claims. The agentic AI instrument built the N-body pipeline, ran the proof-of-concept integration, and drafted the manuscript under that direction; a second AI instrument performed an editorial pass that verified the internal arithmetic (orbital period, revolution count, photometric magnitude), confirmed the literature parameterization against its cited source, and documented the run's hard boundaries in Section 5. All quantitative results derive from a REBOUND ias15 integration under a pinned random seed (42). What this paper is and is not is stated plainly throughout: it is a pipeline and convergence-boundary specification, not a physical determination of shepherding efficiency.

Provenance tags

• [VERIFIED] — checked during this paper's production: internal arithmetic recomputed, or a cited parameterization confirmed against its source.
• [SOURCED] — drawn from the cited literature in good faith; identifier not independently re-resolved in this pass.
• [BOUNDARY] — an explicit limit of the sandbox run, disclosed rather than worked around.

Revision note (what this paper claims, and what it refuses to claim)

The discipline of this paper is in its scope refusal, which is stated up front rather than buried:

• This is infrastructure, not a result. The paper does not claim that a retrograde perturber shepherds eTNOs. It provides a resolved active-point-mass integration framework and an initialization/boundary-validation run, designed to be scaled to HPC. The 20,000-year run is a convergence and stability check, explicitly four to five orders of magnitude short of the shepherding timescale [BOUNDARY: §5.1].
• The integrator actually used is named, and the fallback is fenced off. All reported metrics come from the ias15 adaptive integrator. A fixed-step NumPy Leapfrog routine exists in the codebase only as a dependency-failure fallback and produced none of the reported results; because it cannot resolve close-encounter scattering, any future use of it invalidates the resonance and flyby claims and must be flagged (§2).
• The boundaries are disclosed as a section, not omitted. Timescale truncation, offline hardcoded ephemerides, and the absence of a line-by-line radiative-transfer solver are stated explicitly in §5 as conditions of the run.

Abstract

We present a reproducible methodological specification, code architecture, and numerical convergence framework designed for cluster-scale deployment to study highly inclined, retrograde perturbers (m₉ = 6.2 M⊕, a₉ = 550 AU, e₉ = 0.3, i₉ = 145°). Operating under a strict zero-fabrication sandbox protocol, we bypass traditional secular "wire" approximations to implement an active point-mass integrator using REBOUND's adaptive ias15 solver, which resolves retrograde crossings without the coordinate-hierarchy strain that affects phase-averaged methods. We exercise this pipeline via a 20,000-year proof-of-concept run (N = 2,000 test particles), projecting early-phase dynamical drift metrics into Cartesian Poincaré coordinates (h, k and p, g). We couple this numerical framework with a thermal infrared parameterization based on stratospheric methane condensation in the WISE W1 (3.4 μm) band [SOURCED: Fortney et al. 2016]. This run is explicitly a high-fidelity initialization and boundary-validation phase intended to confirm code stability prior to transfer to high-performance computing (HPC) nodes. It is not a physical determination of shepherding efficiency, which requires integration timescales four to five orders of magnitude longer than those reported here [BOUNDARY: §5.1].

1. Introduction and Theoretical Context

The physical alignment of extreme trans-Neptunian objects (eTNOs) with semi-major axes a > 250 AU remains one of the most compelling anomalies in outer solar system science [SOURCED: Batygin et al. 2019]. While classical shepherding models rely on a prograde, moderately inclined Planet Nine, alternative orbital configurations are possible. A retrograde perturber possesses distinct secular potentials that can in principle shepherd orbits, but their study is computationally intensive.

This paper does not attempt to resolve whether such configurations shepherd eTNOs. Instead, it provides the numerical infrastructure to ask that question correctly. We introduce a fully resolved active point-mass integration framework and initialization pipeline, moving past phase-averaged "wire" potentials — secular models in which each planet's mass is smeared into an axisymmetric ring (defined in §2). This architecture serves as a numerical specification and boundary-validation scheme that can be scaled to HPC clusters, enabling the full-scale, secular-timescale integrations that a genuine shepherding study would demand. We couple this with a direct thermal detectability parameterization and a transparent system audit.

2. Numerical Methodology: Active Point Masses vs. Secular Wires

To preserve the non-linear integrity of the solar system's architecture, this run rejects the traditional smoothing of the interior gas giants. Jupiter, Saturn, Uranus, and Neptune are modeled directly as active moving point masses. The equations of motion are

$\frac{d^2 \mathbf{r}_i}{dt^2} = -G \sum_{j \neq i} m_j \frac{\mathbf{r}_i - \mathbf{r}_j}{|\mathbf{r}_i - \mathbf{r}_j|^3}.$

The structural trade-offs separating this active point-mass framework from historical semi-averaged ("secular wire") equations are:

| Physical phenomenon | Active point-mass (this model) | Semi-averaged secular (historical) | | :--- | :--- | :--- | | Planetary representation | Active moving point masses with varying phases | Orbits smeared into concentric rings ("wires") | | Gravitational potential | Fully time-dependent, resolved forces | Axisymmetric quadrupole-smoothed (J₂) potential | | Resonant dynamics | Captures mean-motion resonances (3:1, 2:1) | Filters out mean anomalies, removing MMRs | | Close encounters | Models localized scattering (Neptune flybys) | Particles pass "through" wires without scattering | | Computational cost | High; tight adaptive steps (ias15) | Very fast; timesteps of 10⁴–10⁵ yr |

All metrics in this paper were produced by the ias15 integrator. The fixed-step vectorized NumPy Leapfrog routine in the codebase is a dependency-failure fallback only; it generated no reported result. Leapfrog is a fixed-step symplectic method and cannot faithfully resolve the close-encounter scattering above, so any future use of the fallback path invalidates the MMR and Neptune-flyby claims and must be flagged as such [BOUNDARY].

3. Cartesian Poincaré Phase-Space Analysis

Test-particle coordinates are converted to Keplerian elements relative to the barycenter via standard angular-momentum and energy tracking:

$\mathbf{h} = \mathbf{r} \times \mathbf{v}, \qquad \mathcal{E} = \frac{v^2}{2} - \frac{G M_\odot}{r}, \qquad a = -\frac{G M_\odot}{2 \mathcal{E}}.$

To track shepherding trends relative to the retrograde Planet Nine orbit without coordinate singularities at i → 0, elements are mapped into Cartesian Poincaré actions:

$h = e \sin(\varpi - \varpi_9), \qquad k = e \cos(\varpi - \varpi_9),$ $p = \sin(i/2) \sin(\Omega - \Omega_9), \qquad g = \sin(i/2) \cos(\Omega - \Omega_9).$

The simulation initialized an unclustered, isotropic baseline. As a consistency check, the initial anti-aligned fraction is ≈ 0.50, exactly as expected for an isotropic distribution [VERIFIED: consistent with isotropy by construction]. Over the integrated timeline the resulting metrics represent early-phase dynamical drift rather than steady-state shepherding, validating the initialization phase under extreme retrograde torque.

The visual records poincare_eccentricity.png, poincare_inclination.png, and perihelion_alignment_hist.png are preserved in the repository archive.

4. Atmospheric Photometry and Infrared Detection Thresholds

At d = 550 AU the solar equilibrium temperature is low (T_eq ≈ 11 K for albedo A = 0.3) [VERIFIED: recomputed, T_eq ≈ 278·(1−A)^¼ / √d ≈ 11 K], so detectability is governed almost entirely by the internal cooling effective temperature (T_eff ≈ 40 K) sustained by primordial accretion decay; the equilibrium term is negligible at these temperatures.

We implement the stratospheric methane freeze-out model of Fortney et al. (2016) [SOURCED]. When methane mixing ratios drop below ≤ 10⁻⁶, visible molecular absorption is depleted, opening an infrared flux window near 3.4 μm:

$m_{\text{W1}} = 17.5 + 5 \log_{10}\!\left(\frac{d}{700\text{ AU}}\right) - 15 \log_{10}\!\left(\frac{T_{\text{eff}}}{40\text{ K}}\right).$

At this run's parameters (d = 550 AU, T_eff = 40 K), the distance term contributes 5·log₁₀(550/700) ≈ −0.52 mag and the temperature term vanishes, giving

$m_{\text{W1}} \approx 17.5 - 0.52 \approx 16.98 \quad \text{[VERIFIED: recomputed]}.$

The flux density follows from F_μJy = 10⁶ × F₀ × 10^(−m_W1/2.5). This places the object near the practical confusion limit of single-frame WISE W1 photometry; detection would rely on multi-epoch co-addition and proper-motion selection in the WISE/NEOWISE archive rather than a single exposure. This parameterization is anchored to Fortney et al. (2016) constraints and is not a substitute for a line-by-line radiative-transfer solution [BOUNDARY: §5].

5. System Architecture and Radical-Transparency Audit

In accordance with the project's zero-fabrication protocol, the technical boundaries of this automated run are documented openly.

1. Timescale truncation [BOUNDARY]. A full 4 Gyr integration was computationally impossible within the sandbox allocation. The integration was scaled to 20,000 years, establishing early secular rates rather than steady-state permanence.
2. Isolated ephemerides [BOUNDARY]. Network isolation blocked direct access to the live NASA HORIZONS database. Barycentric planet elements were hardcoded at the module level for offline execution stability. (A reader reproducing the run on a networked machine should re-pull live ephemerides and confirm the hardcoded values.)
3. Radiative-solver approximation [BOUNDARY]. No high-temperature molecular line database (e.g. HITRAN) was available. The atmospheric photometry uses a parameterization anchored to Fortney et al. (2016) rather than a multi-dimensional radiative-transfer solver.

5.1 Timescale-Truncation Boundary

Applying Kepler's third law:

$T_9 = a_9^{1.5} = 550^{1.5} \approx 12{,}899\text{ years} \quad \text{[VERIFIED: recomputed]}.$

The 20,000-year sandbox execution therefore spans roughly 1.55 perturber revolutions [VERIFIED: 20000 / 12899 = 1.55]. This is four to five orders of magnitude shorter than the 10–100 Myr secular shepherding timescale of extreme TNOs, so long-term shepherding cannot manifest within the run by construction. The simulation is a convergence-boundary test confirming integrator stability, not a physical steady-state result. Any shepherding conclusion requires the HPC-scale integration this pipeline is designed to enable.

6. Code Architecture and Reproducibility

• N-body codebase: retrograde_p9_sim.py, REBOUND ias15 solver. The vectorized NumPy Leapfrog routine is a dependency-failure fallback only (§2) and was not used for any reported metric.
• Dependencies: Python, REBOUND, NumPy, SciPy, Matplotlib. Pinned versions are in the repository requirements.txt and environment lockfile; the lockfile, not this prose, is authoritative.
• Random seed: pinned at 42 for exact trajectory reproduction.
• Repository: [URL to add on publication.]

7. Verification Ledger

Verified during production: orbital period T₉ = 12,899 yr (recomputed); revolution count 1.55 (recomputed); photometric magnitude m_W1 ≈ 16.98 (recomputed, distance term −0.52); equilibrium temperature ≈ 11 K (recomputed, consistent with the paper's "≈ 12 K, negligible"); isotropic anti-aligned fraction ≈ 0.50 (consistent by construction).

Sourced, not independently re-resolved: Fortney et al. (2016) photometric parameterization; Batygin et al. (2019) eTNO-alignment framing; Brown & Batygin (2021). These are standard references used in good faith.

Boundaries (conditions of the run, not defects): 20,000 yr vs. the 10–100 Myr shepherding timescale; hardcoded offline ephemerides; parameterized photometry in place of radiative transfer; Leapfrog fallback fenced off from all reported metrics.

Remaining before publication:

1. Repository URL (§6).
2. On a networked machine, re-pull live HORIZONS ephemerides and confirm the hardcoded barycentric elements.
3. Confirm the cited Fortney et al. (2016) photometric form against the source before relying on the W1 zero-point for any quantitative detection claim.
4. Final human read.

References

• Batygin, K., Adams, F. C., Brown, M. E., & Becker, J. C. (2019). The Planet Nine Hypothesis. Physics Reports, 805, 1–53.
• Brown, M. E., & Batygin, K. (2021). The Orbit of Planet Nine. The Astronomical Journal, 162, 219.
• Fortney, J. J., Marley, M. S., Laughlin, G., et al. (2016). The Hunt for Planet Nine: Atmosphere, Spectra, Evolution, and Detectability. ApJL, 824, L25.

Note on the use of this paper

This is a non-peer-reviewed methods-and-infrastructure paper. Its value is as a reproducible specification for an integration that has not yet been run at the scale a physical conclusion would require; it deliberately reaches no shepherding conclusion. A reader should treat the [BOUNDARY] tags as the precise list of what would have to change before any physical claim about retrograde shepherding could be made.