3.1. Hermite Scheme

The Hermite scheme and hierarchically blocked time-steps (Hut and McMillan Reference Hut and McMillan1986; Makino and Aarseth Reference Makino and Aarseth1992) improved the performance on vector and parallel computers and turned out to be efficient for all of recent parallel and innovative hardware (general and special purpose parallel computers, GRAPE and GPU, see Sect. 3.5 below). The Hermite scheme uses the analytic Nbody expression for the Newtonian gravitational acceleration $\vec{a}_0$ and its time derivative $\vec{{\dot a}}_0$ at some time $t_0$ , to extrapolate low-order predicted positions $\vec{x}_p$ and velocities $\vec{v}_p$ at the next time step $t\gt t_0$ . Then again gravitational acceleration $\vec{a}_1$ and its time derivative $\vec{{\dot a}}_1$ at t are computed, using only predicted data. These four values of $\vec{a}$ and $\vec{{\dot a}}$ on two time points $t_0,t$ allow the the Hermite extrapolation to compute the $2^{\mathrm{nd}}$ and $3^{\mathrm{rd}}$ derivatives of acceleration. These are used finally to obtain the high order corrected new positions and velocities from

(1) \begin{align} \vec{x}(t) &= \vec{x}_p(t) + \frac{1}{24}(t-t_0)^4 \vec{a}_0^{(2)} +\frac{1}{120}(t-t_0)^5 \vec{a}^{(3)} + \mathcal{O}(\Delta t^6).\end{align}

(2) \begin{align}\vec{v}(t) &= \vec{v}_p(t) + \frac{1}{6}(t-t_0)^3 \vec{a}_0^{(2)} +\frac{1}{24}(t-t_0)^4 \vec{a}_0^{(3)} + \mathcal{O}(\Delta t^5).\end{align}

The effective error in the methods is of order $\mathcal{O} (\Delta t^5)$ , as was shown in test simulations (Makino Reference Makino1991), which is not obvious due to the intrinsic low order predicted values used in the Hermite step.

3.2. Time-step choice

Aarseth (Reference Aarseth, Brackbill and Cohen1985) provides an empirical time-step criterion

(3) \begin{equation}\Delta t = \sqrt{\eta \frac{ \vert\vec{a}\vert \vert\vec{a}^{(2)}\vert + \vert\vec{{\dot a}}\vert^2 }{ \vert\vec{{\dot a}}\vert \vert\vec{a}^{(3)}\vert + \vert\vec{a}^{(2)}\vert^2 }} \ .\end{equation}

The error is governed by the choice of $\eta$ , which in most practical applications is taken to be $\eta = 0.01 - 0.04$ . In Makino (Reference Makino1991) a different time step criterion has been suggested, which appears simpler and more straightforwardly defined, and couples the time-step to the difference between predicted and corrected coordinates. The standard Aarseth time-step criterion from Eq. (3) has been used in most N-body simulations so far, because it seems to achieve an optimal step better than (on average) the mathematically more sound Makino step (see the time-step related discussion in Sweatman Reference Sweatman1994). For the purpose of vectorization and parallelization the time steps of particles are quantized, only a variation by a factor of 2 is permitted, and the time t should be an integer multiple of all time steps (Hut and McMillan Reference Hut and McMillan1986).

3.3. Ahmad-Cohen neighbour scheme

For every particle a neighbour radius is defined, accelerations and their time derivatives according to the Hermite scheme are separately computed for particles inside and outside the neighbour radius. Consequently two different time steps are defined, the smaller irregular one related to the more fluctuating neighbour particle’s forces, and another more steady, regular one related to the distant particles. The Ahmad-Cohen (AC) neighbour scheme goes back to Aarseth and Zare (Reference Aarseth and Zare1974) and typically yields regular timesteps, which are 5–10 times larger than the irregular steps. During irregular steps the regular forces are extrapolated using low order Taylor series. This reduces the amount of full force calculations significantly (see Huang et al. Reference Huang, Spurzem and Berczik2016, and also an earlier analysis by Makino and Hut Reference Makino and Hut1988). The AC scheme is currently rarely used elsewhere, with the exception of a recent simulation of star clusters in a cosmological context using Enzo-N (Jo et al. Reference Jo, Kim, Kim and Bryan2024). Note that the AC scheme is algorithmically similar to smoothed particle hydrodynamics (see for a more historical quote Spurzem et al. Reference Spurzem, Berczik, Berentzen, Merritt, Nakasato, Adorf, Brüsemeister, Schwekendiek, Steinacker, Wambsganß, Martinez, Lienhart, Kugel, Männer, Burkert, Naab, Vasquez and Wetzstein2007). Note also that the concept of a neighbour radius is not strictly implemented, there are precautions to search for fast incoming particles in a larger environment already in Nbody5, and nowadays a mass ratio far away from unity (in both directions) triggers larger search radii for neighbours in the current AC scheme (Wang et al. Reference Wang, Spurzem, Aarseth, Nitadori and Berczik2015; Kamlah et al. 2022a; Spurzem and Kamlah Reference Spurzem and Kamlah2023).

3.4. Regularizations

As the relative distance r of the two bodies becomes small, their time-steps are reduced to prohibitively small values, and truncation errors grow due to the singularity in the gravitational potential (von Hoerner Reference von Hoerner1960, Reference von Hoerner1963, Reference von Hoerner2001, see Sect. 2). The idea is to take both stars out of the main integration cycle, replace them by their centre of mass and advance the usual integration with this composite particle instead of resolving the two components. The internal motion of the two members of the regularized pair (KS pair, for Kustaanheimo and Stiefel Reference Kustaanheimo and Stiefel1965) is done in a new 4D coordinate system, which is similar to a quaternion space, including perturbing forces of other nearby stars. The theoretical base of this can be found e.g. in Neutsch and Scherer (Reference Neutsch and Scherer1992); the algorithm including the time transformation is in Mikkola (Reference Mikkola1997).

Close encounters between single particles and binary stars are a central feature of cluster dynamics. The chain regularization (Mikkola and Aarseth Reference Mikkola and Aarseth1998) is invoked if a KS pair has a close encounter with another single star or another pair (see Fig. 3). The chain regularization is effectively a chain of two body KS regularization, with the other particles perturbing the two body forces. Especially, if systems start with a large number of primordial (initial) binaries, such encounters may lead to stable (or quasi-stable) hierarchical triples, quadruples, and higher multiples. They are treated by using special stability criteria (Mardling and Aarseth Reference Mardling and Aarseth2001). Algorithmic regularization is a relatively recent different method based on a time transformed leap-frog method (Mikkola and Merritt Reference Mikkola and Merritt2008); it does not employ the KS transformation. See for its use and application the next subsection and Banerjee (this volume, on Nbody7).

Figure 3.

Sverre Aarseth and Seppo Mikkola, discussing chain regularization, in 2007.

3.5. Nbody with Parallel and GPU computing

The importance of vector and parallel computing was discussed in a pioneering conference (Hut and McMillan Reference Hut and McMillan1986). The GRAPE chip (Sugimoto et al. Reference Sugimoto, Chikada, Makino, Ito, Ebisuzaki and Umemura1990; Makino et al. Reference Makino, Kokubo and Taiji1993) is an application specific integrated circuit (ASIC), which computes in parallel gravitational forces between particles. It also computed the time derivative of the force, to be directly applicable to the Hermite scheme of Nbody6 (Makino et al. Reference Makino, Fukushige, Koga and Namura2003). A first implementation of the Nbody-code on GPU rather than GRAPE was reported by Nitadori and Aarseth (Reference Nitadori and Aarseth2012) using CUDA kernels. Many of these kernels written by Keigo Nitadori are still in current use, even for the massively parallel programs such as Nbody6++GPU.

At the same time another development started, parallelization of Nbody6 with the (at that time) new standard MPI (message passing interface), yielding Nbody6++ (Spurzem Reference Spurzem1999). Nbody6++ parallelizes both force loops with MPI, for the regular and the neighbour force in the Ahmad–Cohen scheme.

With the advent of clusters, where nodes would be running MPI, and each node having a GPU accelerator, Nbody6++GPU was created – on the top level MPI parallelization is done for the force loops (coarse grained parallelization) and at the bottom level each MPI process calls its own GPU to accelerate the regular force calculation (Berczik et al. Reference Berczik, Spurzem, Zhong, Wang, Nitadori, Hamada and Veles2013), and an AVX/SSE implementation on the CPUs accelerates prediction and neighbour (irregular) forces (Wang et al. Reference Wang, Spurzem, Aarseth, Nitadori and Berczik2015). With this code the Dragon-I simulation was done (Wang et al. Reference Wang, Spurzem, Aarseth, Nitadori and Berczik2015, Reference Wang, Spurzem, Aarseth, Giersz and Askar2016), at time of publication the award winning largest direct N-body simulation of a globular cluster with all required astrophysics (single and binary stellar evolution, stellar collisions, tidal field), simulated over 12 Gyrs.

In recent years, inspired also by LIGO/Virgo/KAGRA gravitational-wave detections (Abbott et al. Reference Abbott, Abbott, Abbott, Abernathy and Acernese2016, Reference Abbott, Abe, Acernese, Ackley, Adhikari, Adhikari, Adkins, Adya, Affeldt and Agarwal2024), numerous current updates have been made with regard to stellar evolution of massive stars in singles and binaries (summarized in Kamlah et al. Reference Kamlah, Spurzem, Berczik, Arca Sedda, Flammini Dotti, Neumayer, Pang, Shu, Tanikawa and Giersz2022b; Spurzem and Kamlah Reference Spurzem and Kamlah2023), and with regard to collisional build up of stars (mass loss at stellar collisions allowed) and intermediate-mass black holes (Rizzuto et al. Reference Rizzuto, Naab, Spurzem, Giersz and Ostriker2021; Rizzuto et al. Reference Rizzuto, Naab, Spurzem, Arca-Sedda and Giersz2022; Arca-Sedda et al. Reference Arca-Sedda, Rizzuto, Naab, Ostriker, Giersz and Spurzem2021). The Dragon-II simulations using Nbody6++GPU show nicely that the ensemble of binary black hole mergers observed by LIGO/Virgo/KAGRA (including the most massive ones) can be well reproduced from star cluster simulations (Arca Sedda et al. Reference Arca Sedda, Kamlah, Spurzem, Rizzuto, Naab, Giersz and Berczik2023, 2024a,b). Dragon-III simulations are in progress for globular and nuclear star clusters (Wu et al. Reference Wu, Cho, Spurzem, Wang, Flammini Dotti and Amiri2025, Cho et al., this volume). Recent work has brought the largest published number of particles to $1.8\cdot 10^6$ (Barber and Antonini Reference Barber and Antonini2025; Rantala et al. Reference Rantala, Lahén, Naab, Escobar and Iorio2025), using PeTar and BiFrost codes, respectively (see below Sect. 4). Current construction sites regarding Nbody6++GPU are a better model for neutron stars (Song et al. this volume), and for tidal disruption events near massive black holes (Cho et al. this volume, and Lee et al. Reference Lee, Lee, Kim, Spurzem, Hong and Chung2025). Finally a variant called Nbody6++GPU-massless has been developed which gives the possibility to integrate large numbers of very small mass particles (free floating planets, asteroids, or comets) inside a star cluster (Flammini Dotti et al. Reference Flammini Dotti, Kouwenhoven, Berczik, Shu and Spurzem2025a, b). The current standard code Nbody6++GPU is available via GitHubFootnote ^‡.

Table 1.

Features of Nbody1-7, adopted with some updates from Spurzem and Kamlah (Reference Spurzem and Kamlah2023).

Table 1 Long description

A table comparing features of Nbody1-7. The table has 8 rows and 10 columns. The columns are labeled ITS, ACS, KS, HITS, PN, AR, CC, MPI, and GPU. The rows are labeled Nbody1 to Nbody7. Each cell contains a checkmark, a special symbol, or is empty, indicating the presence or absence of specific features. Notable features include the use of Individual time-steps (ITS), Ahmad-Cohen neighbour scheme (ACS), KS-regularization of few-body subsystems (KS), Hermite scheme and block steps (HITS), Post-Newtonian (PN), Algorithmic regularization (AR), Classical chain regularization (CC), Message Passing Interface (MPI), and GPU acceleration. Special symbols indicate inclusion in special versions of the code or restrictions on certain features.