REBOUNDx: Extension Library for N-body Simulations

• REBOUNDx is an extension library for REBOUND that systematically incorporates both conservative and dissipative astrophysical forces using an operator-splitting approach.
• It employs first-order and Strang symmetric splittings to reduce integration errors and supports custom physical effects via comprehensive C and Python APIs.
• The library ensures machine-independent reproducibility with its SimulationArchive format and has been validated for modeling planetary, stellar, and small-body dynamics.

REBOUNDx is an extension library for the open-source N-body integrator REBOUND, designed to enable the systematic addition of conservative and dissipative astrophysical forces to N-body simulations. It provides a plugin-based operator-splitting framework allowing researchers to incorporate complex physical processes—including tidal dissipation, general relativistic corrections, radiation forces, mass loss, and user-defined interactions—alongside state-of-the-art symplectic integrators such as WHFast and IAS15. REBOUNDx is implemented in C with a comprehensive Python interface, and supports full machine-independent reproducibility via an extended SimulationArchive binary format. Its modularity and rigor in handling velocity-dependent and non-Hamiltonian forces allow accurate modeling of planetary system evolution, stellar-planet interactions, and small-body dynamics across a wide range of astrophysical contexts (Tamayo et al., 2019).

1. Operator-Splitting Framework and Motivation

REBOUNDx addresses the limitations of conventional symplectic N-body integrators, which excel at handling conservative Hamiltonian dynamics but encounter significant errors when naive attempts are made to add dissipative or velocity-dependent perturbations. By moving to an operator-splitting formalism based on Lie derivatives, REBOUNDx enables each physical process—whether conservative (e.g., gravitational multipoles) or dissipative (e.g., tides, radiation drag)—to be coded as a separate operator acting on the system over properly ordered substeps. The first-order and Strang (second-order) symmetric splittings have error terms that are nested commutators of the operators; for perturbations that partially commute with the principal integrator (e.g., in low-eccentricity limits), REBOUNDx exploits geometric symmetries to suppress integration errors. This formalism further supports symplectic and "weak-splitting correctors," which remove the leading error terms up to arbitrary order, and work equally for dissipative systems (Tamayo et al., 2019).

2. Implemented Physical Effects

The REBOUNDx library provides a catalog of pre-implemented astrophysical effects:

Effect	Description	Key Parameters
Migration and Damping	Type I/II disk torques (Papaloizou & Larwood 2000)	$\tau_a,\,\tau_e,\,\tau_i$
General Relativity (1PN)	Both full two-body and central-mass limits	$c$ (speed of light), $\mu$
Radiation Forces	Poynting–Robertson drag, solar radiation pressure	$\beta$, $L_\star$, $R$, $\rho$
Mass Loss	Isotropic (stellar wind, mass evolution)	$\tau_m$ (mass-loss timescale)
Gravitational Multipoles	$J_2, J_4$ moments, arbitrary $\Phi(r)$	$J_2$, $J_4$, $R$
Equilibrium Tides	Hut (1981), constant time lag (CTL)	$k_2$, $\tau$, $R$, $\Omega$
Yarkovsky Effect	Diurnal/seasonal radiative accelerations	see Section 3
Parameter Interpolation	Time-variable parameters from stellar models	timeseries $M(t),R(t),L(t)$
Self-Consistent Spin and Tides	Full coupled spin/orbit evolution	$C$, $\omega_x,y,z$
User-Defined Forces	Arbitrary acceleration or torque functions	arbitrary

This architecture makes it possible to combine multiple effects, and to customize per-particle and global parameters at runtime via both the C API and the Python interface (Tamayo et al., 2019, Baronett et al., 2021, Lu et al., 2023).

3. Advanced Modules: Yarkovsky Effect and Tides

Yarkovsky Effect

The REBOUNDx Yarkovsky module introduces radiative accelerations relevant for small bodies (asteroids, debris) and is accessible in two models (Ferich et al., 2022):

• Full Version ("ye_flag = 0") implements the 3D Yarkovsky matrix

$$\frac{d\mathbf{v}}{dt} = \frac{3\,k\,L\,(1-\alpha)}{16\,\pi\,\rho\,R\,c\,r^2}\;\mathbb{Y}\,i$$

where $\mathbb{Y}$ consists of diurnal/seasonal lag rotations, $k$ is the spin constant, and $i$ describes the aberrated radiation direction. Inputs: radius, density, albedo, spin axis, thermal parameters, stellar luminosity.

• Simple Version ("ye_flag = \pm 1"): constant $\mathbb{Y}$ matrix for either prograde or retrograde drift, neglecting diurnal/seasonal separation and reducing input requirements.

Performance benchmarks indicate the Full Version increases run time by a factor $\sim$2–3, while the Simple Version adds a factor $\sim$1.3–1.8 vs. baseline N-body only (Ferich et al., 2022).

Tidal Modules and Stellar Evolution

The constant time lag (TCTL) tides module implements forces following Hut (1981), with fixed spin aligned along $z$. The self-consistent spin+tide extension (spin_tides) adds ODE evolution of spin vectors and coupled quadrupole and tidal friction torques per Eggleton, Kiseleva & Hut (1998). The operator-splitting approach underlies both, yielding robust long-term evolution of planetary obliquities, synchronous rotation, and orbital decay (Baronett et al., 2021, Lu et al., 2023).

Parameter interpolation allows researchers to load external timeseries (from MESA tracks) for stellar parameters $M_\star(t)$, $R_\star(t)$, $L_\star(t)$ via monotonic cubic spline fitting. This enables full coupling of N-body and stellar evolutionary tracks with overhead typically less than 20% (Baronett et al., 2021).

4. Library Architecture, API, and Extensibility

The REBOUNDx core provides:

• C and Python API: Effects are loaded and parameters are set at runtime via rebx_load_force/rebx_set_parameter (C) or Extras.load_force/.parameters (Python). Custom effects can be registered as arbitrary user-defined callbacks.
• SimulationArchive: An extended binary format stores all active effects, integrator configuration, and parameters per-snapshot, enabling exact reproducibility across platforms.
• Operator Scheduling: Each effect defines functions for pre-timestep, force, or post-timestep application; for symplectic integrators, dissipative steps are applied in the “kick” (interaction) stage.
• Corrector Support: For WHFast, symplectic and weak-splitting correctors (up to third order) can be activated by setting sim.ri_whfast.correctors (Python), yielding orders-of-magnitude improved energy conservation even with dissipative effects.

Sample code implementations are included for both C and Python environments, and demonstrate initialization, effect loading, parameter application, and integration loop structure (Tamayo et al., 2019).

5. Benchmarks, Performance, and Limitations

• Performance: For common planetary architectures (4 terrestrial planets, Myr–Gyr integrations), adding multiple REBOUNDx effects incurs a fractional overhead of ~18% (PI + TCTL). The Full Yarkovsky model increases runtime by a factor 2–3, while Simple Yarkovsky and standard Radiation Forces increase time by ~1.3×. Spin+tide ODE integrations are dominated by the N-body loop for $N \gtrsim 5$ (Ferich et al., 2022, Baronett et al., 2021, Lu et al., 2023).
• Timestep Guidelines: Good practice is to pick an integrator step size $\lesssim 0.1$ times the shortest dynamical period in the system; update intervals for parameter interpolation can be $\sim$10–100 yr for post-main sequence stellar evolution (much shorter for rapidly changing tracks).
• Physical Assumptions and Constraints:
  • Yarkovsky routines assume spherical, uniform-density bodies; spin evolution due to thermally-driven torques (YORP) is not included.
  • Tides may be calculated on a primary only (secondaries do not back-react), and constant time-lag implementation fixes spins if not using spin+tide ODEs.
  • For models requiring non-adiabatic stellar changes (rapid mass-loss), splitting errors can rise unless very small intervals are set (Baronett et al., 2021).
  • Code design allows N-body particles to lack subset parameters (e.g., spin, $k_2$), in which case corresponding effects are ignored for that particle.

6. Applications and Use Cases

REBOUNDx supports a broad range of astrophysical studies:

• Solar System and Exoplanet Post-MS Evolution: By combining parameter interpolation and tidal modules, REBOUNDx accurately models planetary orbital expansion or engulfment during stellar giant phases, with quantitative agreement to analytic and previous simulation results. Inner, Earth-mass planets can be tidally destroyed out to 0.8–1.4 au; outer giants expand adiabatically per expected analytic scalings (Baronett et al., 2021, Lu et al., 2023).
• Small-Body and Dust Belt Dynamics: Yarkovsky and radiation effects enable the study of asteroid belt evolution, orbital drift, and debris disk shaping under strong stellar luminosity changes (e.g., AGB phases), with robust treatment across particle size and spin regimes (Ferich et al., 2022).
• Spin-Orbit Evolution and Obliquity Dynamics: Self-consistent spin and companion tidal friction modules allow simulation of obliquity excitation in multi-resonant systems, pseudo-synchronization, and tidal migration scenarios—matching analytic predictions for spin equilibrium and dynamical histories (Lu et al., 2023).
• Custom Physics: Users can inject arbitrary velocity or position-dependent forces (e.g., gas drag, user-defined migration laws) via both high-level Python and low-level C, and mix conservative and dissipative operators as needed.

7. Distribution, Documentation, and Best Practices

REBOUNDx is freely available via GitHub and PyPI. Versioned documentation includes full API references, effect formulas, and implementation guides. Recommended practices include:

• Use corrector order ≥2 when energy conservation is critical and perturbation strengths are moderate.
• Monitor the stability of integrations, adapting timesteps especially where strong dissipative or highly time-variable effects dominate.
• For coupled N-body/stochastic evolution (parameter interpolation), choose interval steps much shorter than external evolution timescales.
• Leverage the SimulationArchive for collaborative, bitwise-reproducible results when sharing simulation outputs (Tamayo et al., 2019).

In summary, REBOUNDx has established itself as a robust, extensible, and rigorously designed library for advancing precision simulation in planetary, stellar, and debris-disk system dynamics (Tamayo et al., 2019, Ferich et al., 2022, Baronett et al., 2021, Lu et al., 2023).