REBOUND: An open-source multi-purpose N-body code for collisional dynamics

REBOUND is a new multi-purpose N-body code which is freely available under an open-source license. It was designed for collisional dynamics such as planetary rings but can also solve the classical N-body problem. It is highly modular and can be customized easily to work on a wide variety of different problems in astrophysics and beyond. REBOUND comes with three symplectic integrators: leap-frog, the symplectic epicycle integrator (SEI) and a Wisdom-Holman mapping (WH). It supports open, periodic and shearing-sheet boundary conditions. REBOUND can use a Barnes-Hut tree to calculate both self-gravity and collisions. These modules are fully parallelized with MPI as well as OpenMP. The former makes use of a static domain decomposition and a distributed essential tree. Two new collision detection modules based on a plane-sweep algorithm are also implemented. The performance of the plane-sweep algorithm is superior to a tree code for simulations in which one dimension is much longer than the other two and in simulations which are quasi-two dimensional with less than one million particles. In this work, we discuss the different algorithms implemented in REBOUND, the philosophy behind the code’s structure as well as implementation specific details of the different modules. We present results of accuracy and scaling tests which show that the code can run efficiently on both desktop machines and large computing clusters.

💡 Research Summary

The paper presents REBOUND, an open‑source, highly modular N‑body simulation framework designed primarily for collisional dynamics such as planetary rings, but also capable of handling classical collision‑less N‑body problems. Written entirely in ISO C99, REBOUND compiles on any POSIX‑compatible system without external dependencies, though optional OpenGL/GLUT and libpng support real‑time visualization and screenshot saving.

The code’s architecture is built around four interchangeable modules—gravity, collision, boundary, and integrator—selected via symbolic links in the Makefile. This design allows users to switch physical models or add new ones without modifying the core source, facilitating simultaneous development of multiple projects.

Time integration is performed with three second‑order symplectic schemes, all following a Drift‑Kick‑Drift (DKD) pattern:

• Leap‑frog – a simple, non‑rotating‑frame integrator suitable for generic collisional systems.
• Wisdom‑Holman (WH) – a mixed‑variable map that solves the Keplerian drift exactly; ideal for systems dominated by a central 1/r potential but slower because Kepler’s equation must be solved iteratively each step.
• Symplectic Epicycle Integrator (SEI) – tailored to Hill’s approximation and shear‑periodic (shearing‑sheet) boundaries; it retains symplecticity in a rotating frame while being as fast as a standard non‑symplectic integrator.

Gravity can be computed either by direct summation (O(N·N_active), efficient only for a few massive particles) or by a Barnes‑Hut octree (O(N log N)). The octree groups particles hierarchically; distant groups are approximated by their total mass and centre of mass, with an opening‑angle criterion θ_crit controlling the accuracy‑speed trade‑off. A compile‑time flag enables quadrupole moments for higher accuracy at additional cost.

Collision detection is offered through two algorithms. The octree‑based method uses the same spatial hierarchy as the gravity tree to limit pairwise checks to nearby particles. The plane‑sweep algorithm sorts particles along a chosen axis and sweeps a sliding slab, detecting overlaps only within the slab. This method dramatically outperforms the tree approach when the simulation domain is highly elongated in one dimension or when the system is quasi‑two‑dimensional with fewer than ~10⁶ particles.

Parallelization combines MPI and OpenMP. MPI employs static domain decomposition: each MPI rank owns one or more root boxes and exchanges only the “essential tree” (the minimal subset of the octree needed by other ranks), minimizing communication. OpenMP provides shared‑memory thread parallelism within each rank. A hybrid MPI+OpenMP mode is supported, and scaling tests demonstrate near‑linear speed‑up up to 1024 cores.

Accuracy tests cover energy and angular‑momentum conservation, collision handling, and force error versus θ_crit and quadrupole order. All integrators show the expected t⁻² energy error scaling for sufficiently small timesteps, while the octree’s force error can be tuned by adjusting θ_crit. The plane‑sweep collision detector shows 5–10× speed gains in the regimes mentioned above.

REBOUND is released under GPL‑v3 and hosted on GitHub (https://github.com/hannorein/rebound). The repository includes a rich set of example problems, a Doxygen‑generated documentation system that updates automatically with the selected modules, and a test suite for verification. By providing a flexible, well‑tested, and openly available platform, REBOUND enables reproducible research across a wide range of astrophysical and non‑astrophysical applications where collisional N‑body dynamics are essential.