Nii-body: Bayesian Inference of Multiplanet Dynamics via N-body Simulations

Abstract: Many exoplanetary systems are multiplanet configurations whose long-term dynamics are governed by N-body gravitational interactions. Consequently, their detection signatures cannot be adequately described by Keplerian orbits. Accurately interpreting the observational data of these systems -- including radial velocity (RV), astrometry, and transit timing variations (TTVs) -- requires N-body integration. Motivated by this need, we developed a Bayesian fitting framework that couples N-body integration with Markov chain Monte Carlo (MCMC) to retrieve the system parameters of multiplanet systems. The code, named \texttt{Nii-body}, integrates an adaptive Runge--Kutta--Fehlberg 7(8) (RKF78) solver with an automated parallel tempering MCMC algorithm. Using simplified synthetic astrometric observations, we evaluated the efficiency and robustness of \texttt{Nii-body}'s N-body orbit retrieval on an idealized two-planet model, demonstrating its potential for future application to real observational data. The N-body fitting workflow can be readily extended to RV, TTVs, or combined datasets, providing a versatile engine for high-precision orbital inference in multiplanet systems.

• The paper introduces a framework that integrates high-accuracy N-body simulations with Bayesian MCMC estimation to model exoplanet dynamics.
• It employs an adaptive RKF78 integrator for machine-precision propagation and synthetic signal generation, validated against established benchmarks.
• The methodology outperforms traditional Keplerian models by reliably recovering orbital parameters and addressing multimodal ambiguities in planetary systems.

Bayesian Inference of Multiplanet Dynamics with Nii-body

Motivation and Context

Multiplanet systems are ubiquitous in exoplanetary discoveries, often exhibiting complex architectures and strong dynamical coupling between components. Accurate inference of their orbital parameters and dynamical stability from observational data (RV, astrometry, TTVs) demands that forward models encapsulate $N$-body gravitational interactions self-consistently. Traditional Keplerian superposition schemes, commonly used in RV and astrometry orbit fitting, are fundamentally limited in scenarios exhibiting significant planet–planet perturbations, especially near mean motion resonances, leading to biased or entirely incorrect parameter inference.

"Nii-body: Bayesian Inference of Multiplanet Dynamics via N-body Simulations" (2604.09383) introduces a computational framework for addressing these challenges by tightly integrating high-accuracy $N$-body orbit integration with modern Bayesian parameter estimation through Markov chain Monte Carlo (MCMC) sampling. The framework is tailored for exoplanet inference, with extensions possible for any dynamical gravitational system.

Numerical Implementation and Integration Scheme

Nii-body employs an embedded Runge–Kutta–Fehlberg 7/8 (RKF78) integrator, implemented in C, to propagate system state vectors. The solver features adaptive stepsize and local error control, offering both high performance and machine-precision solutions necessary for robust sampling over high-dimensional, stiff dynamical models. The RKF78 implementation was rigorously benchmarked against the widely-adopted REBOUND/IAS15 integrator, demonstrating negligible discrepancies ($\mathrm{RMSE} \sim 2 \times 10^{-14}$ AU over $10^4$ yr for the Sun–Earth–Jupiter system), at a competitive computational cost.

Figure 1: Discrepancies in $x$, $y$, and $z$ coordinates between Nii-body/RKF78 and REBOUND/IAS15 trajectories, demonstrating numerical equivalence over 10,000 years.

State vector propagation is performed in barycentric Cartesian coordinates, with seamless conversion utilities to and from Keplerian elements $(P, e, i, \omega, \Omega, M_0)$ to interface with observation generation modules. The solver enables native production of synthetic astrometric and RV signals, fully capturing the system geometry and orientation.

Figure 2: Geometry of synthetic astrometric and RV signals, with observer location and projection conventions.

Synthetic Astrometry and Limitations of Keplerian Models

The importance of direct $N$-body integration is elucidated via comparison to the Keplerian superposition approach in a two-planet model of the Kepler-9 system, which contains a b–c pair in a 2:1 mean motion resonance. The Keplerian model fails to reproduce the secular drift and amplitude modulations present in the N-body solutions, as highlighted by the residuals between the two prescriptions. The differential signal is already substantial over five years and will amplify over longer baselines, further confirming the need for self-consistent $N$-body forward models in resonant and tightly packed multiplanet configurations.

Figure 3: Astrometric wobbles of the Kepler-9 host star—comparison of N-body (red) and Keplerian superposition (blue), with their residuals.

Bayesian Fitting Pipeline and MCMC Convergence

Bayesian inference is performed by coupling the RKF78 integrator to Nii-C, an automated parallel tempering MCMC engine supporting high-dimensional, multimodal posterior exploration. The sampler optimizes the proposal distributions during burn-in and robustly identifies the global posterior mode through temperature-stratified chains. Posterior convergence is stringently monitored over multiple independent runs using the Gelman–Rubin criterion.

Astrometric likelihoods are computed assuming independent, isotropic Gaussian errors; future extensions will incorporate covariance structures for application to real datasets. For the synthetic Kepler-9 system, the retrieval task involves 15-parameter posteriors (two planet masses, $N$0 orbital elements, $N$1 observational noise scale). Notably, the sampler achieves convergence from wide, uninformed priors—i.e., it is not reliant on precise initialization.

Figure 4: Marginal posteriors for all model parameters from one MCMC chain, illustrating single-mode recovery for $N$2 and $N$3.

Empirical Validation: Kepler-9 Astrometry

Astrometric time series were synthesized from the N-body model and injected with Gaussian observational noise ($N$4as). Nii-body successfully recovered the full set of input parameters to sub-percent accuracy in 6 out of 8 parallel MCMC runs (75% recovery rate), with the remaining two either falling into local degeneracies or displaying insufficient mixing due to multimodal structure inherent to the posterior. In contrast, none of the Keplerian superposition-based retrievals converged to the correct modes; their posteriors were broad and systematically biased.

Figure 5: Comparison of N-body ground-truth wobbles (black), noisy synthetic observations (dots), and MCMC-retrieved best fit (red), demonstrating recovery fidelity.

Nii-body offers quantitative posterior credible regions—the only meaningful approach to inference in such nonlinear, degenerate parameter spaces. The sampler also identifies and flags multimodal ambiguities, such as the degeneracy between $N$5 and $N$6, which are intrinsic to astrometry and require, e.g., RV or TTV data fusion to resolve.

Practical Performance

Fitting the two-planet synthetic Kepler-9 dataset in 15 parameters achieves reliable convergence within $N$7 hours on a commodity 12th-gen Intel laptop (with $N$8 MCMC steps). Scaling tests indicate that computational cost grows with the number of planets and epochs, but remains feasible for systems up to $N$9–$\mathrm{RMSE} \sim 2 \times 10^{-14}$0) planets and several hundred epochs, making the method applicable to ongoing and future high-precision astrometric missions such as Theia and CHES.

Implications and Future Work

Nii-body unambiguously demonstrates that Keplerian superposition models are insufficient for interpreting high-precision observations of dynamically active multiplanet systems; N-body orbital integration is required to avoid biased inference. This is especially pertinent for exoplanet population statistics, dynamical stability assessment, and precise mass/inclination determinations. The architecture admits straightforward inclusion of RV and TTVs for joint dynamical modeling.

Future development will address:

• Incorporation of parallax, proper motion, and more realistic astrometric error models to enable analysis of real datasets.
• Dynamical stability constraints and long-term integration of posterior samples to rule out unphysical configurations.
• Scalability enhancements, including hierarchical modeling over ensembles and acceleration via GPU-parallelization or symplectic mapping for cold systems.

Practically, Nii-body serves not only the exoplanet characterization community but can be adapted to any hierarchical gravitational inference problem (stellar clusters, galactic cores, etc.) where N-body effects are non-negligible.

Conclusion

Nii-body provides a rigorously validated, open source pipeline for Bayesian inference of multiplanet dynamics with self-consistent N-body integration, robust MCMC-based posterior estimation, and native synthetic data generation for astrometry and RV. Direct comparison with Keplerian approaches shows that only N-body forward modeling captures the signatures of dynamical perturbations, enabling reliable recovery of true system parameters in regimes relevant to current and next-generation exoplanetary observations. Its systematic application will be crucial in the dynamical mapping of planetary systems where gravitational interactions are dynamically significant and must not be neglected.