REBOUND's IAS15 Integrator

• REBOUND's IAS15 integrator is a high-precision 15th-order N-body method that employs Gauss–Radau quadrature and a predictor–corrector loop to achieve near machine-level accuracy.
• It features adaptive timestep selection based on fractional error control, efficiently handling close encounters, high-eccentricity orbits, and non-conservative forces.
• Implemented in C99 with a Python interface, IAS15 is versatile for simulating diverse gravitational dynamics, including planetary systems and resonant interactions.

REBOUND’s IAS15 integrator is a high-precision, 15th-order N-body integrator based on Gauß–Radau quadrature, designed for gravitational dynamics in both conservative and non-conservative systems. Developed for integration accuracy near the machine precision limit across up to $10^9$ orbits, IAS15 automatically adapts its timestep to system dynamics and is implemented within the open-source REBOUND framework in C99, with a Python interface. It efficiently handles close encounters, high-eccentricity orbits, and non-conservative forces such as radiation pressure and damping. Applications include long-term integrations of planetary systems, Kozai–Lidov oscillations, close encounters, and scenarios requiring both extremely high accuracy and adaptability.

1. Algorithmic Principles and Predictor–Corrector Architecture

IAS15 employs a 15th-order integration scheme built on the Gauß–Radau quadrature, which leverages seven internal substeps (nodes $h_1, ..., h_7$ plus the starting time $h=0$) per time step. The method expands the acceleration $\mathbf{y}''(t)$ for each particle in a truncated series: $$\mathbf{y}''[t] \approx \mathbf{y}''_0 + a_0 t + a_1 t^2 + \ldots + a_6 t^7$$ Translating to normalized time $h = t/dt$ and using projected coefficients $b_k = a_k dt^{k+1}$, this becomes: $$\mathbf{y}''[h] \approx \mathbf{y}''_0 + b_0 h + b_1 h^2 + \ldots + b_6 h^7$$ An alternate, operationally convenient expansion involves computing $g_k$ coefficients directly from force evaluations at the Gauß–Radau nodes, which correspond to algebraic combinations of acceleration differences divided by polynomial terms in $h_j$. The quadrature weights and nodes are computed in extended precision to suppress scheme error $E_{\text{scheme}}$ below floating-point limits.

Each integration step proceeds via a predictor–corrector loop:

• Prediction: Assume a provisional state (e.g., $b_k=0$) to estimate substep positions and velocities.
• Force Evaluation: Calculate force at each substep via $F[\mathbf{y},\dot{\mathbf{y}},t]$.
• Correction: Update $g_k$ and thus $b_k$ by comparing predicted and calculated accelerations.
• Iteration: Repeat prediction and correction until the highest-order $b_6$ coefficient change is below the threshold (typically $10^{-16}$). Most steps converge in two iterations.

2. Adaptive Timestep Selection and Robust Error Control

IAS15’s timestep control implements a physics-driven, dimensionless error estimator. The core step-size selection equation is: $$dt_{\text{required}} = dt_{\text{trial}} \cdot \left( \frac{\epsilon_b}{\tilde{b}_6} \right)^{1/7}$$ where $\epsilon_b$ is a user-defined tolerance (default $\sim 10^{-9}$), and

$$\tilde{b}_6 = \frac{\max_i |b_{6,i}|}{\max_i |y''_i|}$$

measures the fractional error in the top (7th-degree) coefficient. The $1/7$ exponent derives from $b_6$'s scaling as $dt^7$.

This step-size selection is invariant under unit rescalings, as no absolute scale enters the criterion. The approach maintains system-agnostic performance and ensures smooth interpolation of acceleration, allowing the timestep to vary flexibly and reliably capture dynamical changes such as close encounters or high-eccentricity excursions.

3. Floating Point and Systematic Error Minimization

IAS15 employs compensated summation to suppress floating-point accumulation (both random and bias components) when updating positions and velocities. All arithmetic on step coefficients, constants, and quadrature nodes is performed in extended (or, in some implementations, arbitrary) precision. As a result, systematic integration error is driven well below the floating-point floor ($E_{\text{scheme}} < E_{\text{floor}}$), so that long-term energy error growth obeys Brouwer’s Law: $\text{Error}_{\text{energy}} \propto \sqrt{t}$ There is no detectable systematic drift in energy error even after $10^9$ orbits; the random-walk nature ensures that error growth resembles an unbiased accumulation of round-off, rather than a secular drift. This optimal error behavior is a critical distinction in long-term n-body simulations.

4. Performance, Accuracy, and Comparison with Symplectic Methods

In typical problems, IAS15 achieves near machine-precision errors ($\sim 10^{-16}$) using only $\sim 100$ time steps per orbit. When comparing against mixed-variable symplectic (MVS) integrators such as Wisdom–Holman, or higher-order SABA-type schemes (Rein et al., 2019), IAS15 is found to be superior in both speed and accuracy when required precision is $10^{-8}$ or better. Unlike symplectic integrators, which preserve Hamiltonian structure exactly for fixed step sizes and conservative systems, IAS15 is agnostic to non-conservative or velocity-dependent forces and therefore suitable for a broader class of gravitational and astrophysical simulations.

IAS15's robust timestep adaptation allows efficient exploration of close encounters without step-size preselection issues faced by symplectic or fixed-step integrators. On the other hand, for quasi-Keplerian, long-term integrations in well-separated planetary systems, recent advances in high-order symplectic methods may afford even higher speed at moderate accuracy, though with much less flexibility for handling non-Hamiltonian physics (Rein et al., 2019).

5. Mathematical Formulation and Error Estimates

IAS15 directly solves the second-order ODE: $\mathbf{y}'' = F[\mathbf{y}, \mathbf{y}', t]$ with $F$ supporting arbitrary dependence on position, velocity, and time. Gaussian quadrature theory provides an upper bound on local integration error: $E_y \approx \mathcal{C}_8 \, y^{(16)}[\xi] dt^{16}$ where $\mathcal{C}_8 \approx 1.15 \times 10^{-21}$ and $y^{(16)}$ is the 16th derivative of the integrated variable. For harmonic oscillators: $$\frac{E_y}{y_0} \approx \mathcal{C}_8 \, \omega^{16} dt^{16}$$ This formalism quantifies the suppression of scheme error inherent in the high quadrature order.

The global error over many timesteps accumulates as: $$E_\text{tot} = E_\text{floor} + E_\text{rand} + E_\text{bias} + E_\text{scheme}$$ with $E_\text{floor}$ (machine precision floor) and $E_\text{rand}$ (random-walk error) ultimately dominating. $E_\text{bias}$ and $E_\text{scheme}$ are eliminated, or rendered negligible, by the integrator construction.

6. Implementation, Extensibility, and Use Cases

IAS15 is implemented as part of the REBOUND codebase, available open-source and written in C99. Its Python API enables rapid scripting and integration in large simulation pipelines. The package includes extendable examples spanning resonant planetary chains, Kozai–Lidov cycles, close-encounter scattering, and scenarios involving generic damping (e.g., planet–disc migration), radiation pressure, and quadrupole moments.

IAS15 is particularly designed to handle:

• Astrophysical systems specializing in close approaches or requiring rapid, localized error reduction,
• Dynamical studies of planetary architectures with non-conservative processes (e.g., migration, tides),
• Mixed gravitational and non-gravitational forces including PR drag, radiation, and custom damping,
• Problems where accurate phase conservation and robustness against floating-point drift are essential for billions of orbits.

Its generality, robust error control, and adaptive features make IAS15 viable for any N-body simulation (regardless of force complexity or orbital architecture) where precision is the primary requirement.

7. Application Domains and Role within the REBOUND Ecosystem

IAS15 is central to the REBOUND package’s offering of N-body integrators. While for some applications—such as long-term, fixed-step integrations in Keplerian-dominated, low-perturbation regimes—alternatives like high-order symplectic integrators, Encke-based methods (Hernandez et al., 2020), or operator-splitting hybrids (Lu et al., 6 May 2024) may yield higher efficiency, IAS15 is uniquely positioned as the integrator of choice for:

• Simulations involving arbitrary force models (including arbitrary user-defined external or dissipative forces),
• Scenarios where automatic control of step size is crucial (for example, during rapid dynamical transitions), and
• Comprehensive, high-accuracy, and extensible orbital dynamics studies where non-conservative physics or close encounters must be robustly accommodated.

IAS15 thus underpins a wide range of research in planetary dynamics, minor planet evolution, tidal system modeling, and general N-body applications requiring maximal accuracy and minimal systematic error growth over arbitrarily long timescales.