REBOUND: Modular N-body Simulator

• REBOUND is a modular, open-source N-body simulator enabling high-fidelity simulations of gravitational and collisional dynamics in various astrophysical contexts.
• Its architecture features independently selectable modules for integrators, gravity solvers, and collision detection, ensuring scalability, reproducibility, and efficient code reuse.
• Advanced integrators like IAS15 and extension libraries such as REBOUNDx support simulation of complex effects including relativistic, tidal, and radiative forces for detailed system analysis.

REBOUND is a modular, open-source N-body simulation platform designed for precision gravitational and collisional dynamics across astrophysical and planetary contexts. With an extensible architecture, REBOUND couples advanced symplectic and adaptive integrators with gravity solvers, efficient collision detection, and a range of physical effect modules via core and extended libraries such as REBOUNDx. Its design enables high-fidelity simulations of planetary systems, planetesimal disks, collisional rings, exomoons, ejecta debris fields, and incorporates conservative and dissipative physics, including radiative, relativistic, and tidal forces.

1. Architecture and Core Framework

REBOUND is written in C99 with a Python API, using a plugin-style modular system in which each algorithmic component—integrator, gravity, collision, and boundary condition—is defined as an independent module selectable at runtime or compile time (Rein et al., 2011). This architecture supports high modularity, code reuse, and reproducibility. All per-particle state is managed through central data structures, with simulation workflow managed via a global reb_simulation object and integrator abstraction. The package is parallelized using both MPI (domain decomposition, distributed essential trees) and OpenMP (threaded force/collision loops), and includes bitwise-reproducible archiving for simulation checkpointing and restart.

2. Gravitational Dynamics and Integrators

REBOUND implements several gravitational solvers, from direct N² summation to an optimized Barnes–Hut octree (with monopole/quadrupole moments; O(N log N) scaling), supporting both self-gravitating and test particle configurations (Rein et al., 2011). N-body evolution is performed with a suite of integrators:

• Leap-Frog (DKD splitting): Standard second-order symplectic, optimal for non-rotating or collision-dominated regimes.
• Symplectic Epicycle Integrator (SEI): Exact epicyclic analytic drift in the Hill/sheet approximation, critical for ring and disk applications.
• Wisdom–Holman Mapping (WH/WHFast/WHFast512): Keplerian+interaction splitting for planetary systems, with machine-precision Kepler solvers and extensions for SIMD acceleration (AVX512 vectorization in WHFast512 yielding up to 4.7× speedup over traditional WHFast; 1.4 days to 5 Gyr for 8-planet system on modern CPUs) (Javaheri et al., 2023).

Advanced symplectic schemes—high-order corrector and kernel methods (WHC$_p$, WHCKL), positive- and mixed-coefficient compositions (SABA2/3/4/SABAC$_p$/SABACL$_p$), and very high-order negative-timestep methods (SABA(10,6,4))—are provided for accuracy demands up to $\Delta E/E\sim10^{-14}$ (Rein et al., 2019). Wisdom–Holman schemes are best for long-term planetary evolution, while the SABA(10,6,4) family yields the minimal energy error at fixed CPU cost in highly stable configurations (see error and runtime tables).

IAS15, REBOUND's 15th-order adaptive integrator, achieves machine-precision accuracy with robust adaptive timestep control (Pham et al., 5 Jan 2024). Recent improvements introduce a new timestep criterion based on acceleration, jerk, and snap, yielding optimal timestep selection for both pericenter passages and highly eccentric apocenter phases regardless of central offset, substantially widening practical precision limits to sub-kilometer/moonlet orbits at Kuiper-belt distances.

3. Collisional Dynamics and Specialized Modules

REBOUND supports multiple collision detection modules: brute-force direct overlap, Barnes–Hut tree-based neighbor search, and computationally optimal plane-sweep algorithms in both cartesian (x) and cylindrical (azimuthal) coordinates for disks/rings and low-dimensional geometries (Rein et al., 2011). These collision finders can be selected to match problem aspect ratio and particle number for optimal performance.

A state-of-the-art collisional fragmentation module, based on Leinhardt & Stewart (2011) and Asphaug (2010) scaling laws, resolves all possible outcome regimes (accretion, erosion, hit-and-run, catastrophic disruption). Largest remnant mass is determined via a universal law as a function of impact energy and disruption threshold, and fragments are created and tracked explicitly (Childs et al., 2022). The expansion factor $f$ (inflation of all radii) can be tuned to increase collision rates at fixed $N$; empirical dynamical time scaling follows $t'\sim t f^{1.5}$ for fragmentation runs. A separate composition-tracking code post-processes collision logs to track multi-species volatile delivery and radial mixing events throughout terrestrial planet assembly.

For debris and ejecta off small-body surfaces, the pure-Python RED extension provides initialization of ejecta cones with power-law and scaling-law distributions of size/velocity, solar radiation pressure, ellipsoidal gravity, and binary/triple compact system geometry (Larson et al., 2021). RED’s direct integration with REBOUND's API enables fast prototyping; typical scaling is O($N^2$) for gravity, with O($N$) overhead from advanced force models (SRP, ellipsoidal potential), and practical runs reach $N=10^5$ on cluster nodes.

4. Extended Physical Effects via REBOUNDx

REBOUNDx is the extension library for adding arbitrary conservative and dissipative operator steps to REBOUND simulations (Tamayo et al., 2019). It provides both force-based and operator splitting interfaces, supporting position- and velocity-dependent effects through high-accuracy stepping (RK4, implicit midpoint as needed). Effects supplied out of the box include:

• General Relativity (1PN, full or dominant-central-mass approximations)
• Tidal forces (constant time-lag, dissipation via Hut 1981 formalism)
• Stellar mass/radius/luminosity evolution (parameter interpolators): Coupling to MESA/SSE output via cubic spline interpolation, supporting operator splitting with guaranteed $\mathcal{O}(Δt^2)$ error (Baronett et al., 2021).
• Radiation pressure and Poynting–Robertson drag
• J$_2$, J$_4$ gravitational moments
• Yarkovsky effect: Two-module approach with a full model (diurnal/seasonal lag, thermal parameterization) or a reduced model (constant-rotation matrix), supporting stellar luminosity time-dependence and robust for both Solar System and post-main sequence star scenarios (Ferich et al., 2022).

All force modules are bitwise-reproducible and their integration is fully compatible with REBOUND's symplectic stepping, with optional symplectic corrector “weak splitting” extensions to further minimize secular drifts arising from external operator commutators.

5. Hierarchical and Secular Stability Analysis

REBOUND excels at hierarchical system construction and long-term stability mapping, including applications to exomoon, submoon, and multiple-satellite systems (Patel et al., 21 Jan 2025). Arbitrary hierarchy can be constructed by recursively specifying the primary for each body, with explicit coordinate rotation when defining multi-level inclinations. Stability boundaries are assessed via Hill radius crossings, unbound criteria ($e≥1$), or collision events. Quantitative mapping of stability as a function of mutual semi-major axes and inclinations employs systematic grids, with phase randomization to sample mean anomaly space. Chaos and resonance identification tools (MEGNO, secular resonance analysis via linear fits to apsidal rates, Laplace–Lagrange secular theory comparison) are included in the post-processing framework.

Parameter sweeps and color-mapped stability diagrams allow direct visualization of the intersection between stability domains and secular resonance bands, with Lyapunov times and resonance widths matched to analytical eigenfrequency combinations.

6. Variational Equations, Optimization, and Chaos Diagnostics

REBOUND's variational module supports first and second-order variational equation integration alongside the primary N-body trajectory using IAS15 (Rein et al., 2016). This enables derivative-based chaos quantification (MLE, MEGNO), direct calculation of Jacobians and Hessians for optimization (e.g., Newton's method for orbit fitting, trajectory optimization, or Riemann-manifold MCMC), and machine-precision evaluation of phase-space divergence at vanishing computational error compared to finite-difference approaches. Analytical initialization of partial derivatives eliminates “step-tuning” and guarantees formal correctness for all common orbital parameters and mass.

Best practices include always using analytic “vary” routines for variational initialization, restricting interval lengths in strongly chaotic orbits, and leveraging high-level Python hooks in the REBOUND API for flexible adoption in fitting pipelines or global search→local refinement workflows.

7. Performance, Scalability, and Best Practice Recommendations

REBOUND benchmarks across direct, tree, and sweep algorithms exhibit near-linear strong and weak scaling for $N\gtrsim2,000$ per MPI rank, with bottlenecks only in essential tree communication at very large processor counts (Rein et al., 2011). For small-N, OpenMP outperforms MPI; for massive runs, pure MPI dominates. Plane-sweep collision detection outperforms all tree-based approaches for low-dimensional or highly elongated geometries, and the octree’s quadrupole expansion yields $1$–$2$ orders of magnitude accuracy improvement in force calculation at modest opening angles.

Choosing an integrator and timestep should be driven by accuracy requirements, dynamical regime, and physics (e.g., SABA(10,6,4) or IAS15 for strongly interacting, high-eccentricity systems; WHCKL for long-term, moderately accurate planetary systems). Always run short IAS15 reference integrations to calibrate secular cycles and resonance features. For operator splitting, document all non-gravitational effect module choices, their parameters, and updating cadences to enable full reproducibility of dynamical boundaries.

For the latest feature-specific modules and usage examples (including Python scripts for hierarchically nested bodies, complex ejecta initialization, collisional fragmentation, or REBOUNDx operator chaining), users should refer to the respective GitHub repositories and published test suites as cited in the primary literature.