How well do STARLAB and NBODY4 compare? I: Simple models

N-body simulations are widely used to simulate the dynamical evolution of a variety of systems, among them star clusters. Much of our understanding of their evolution rests on the results of such direct N-body simulations. They provide insight in the structural evolution of star clusters, as well as into the occurrence of stellar exotica. Although the major pure N-body codes STARLAB/KIRA and NBODY4 are widely used for a range of applications, there is no thorough comparison study yet. Here we thoroughly compare basic quantities as derived from simulations performed either with STARLAB/KIRA or NBODY4. We construct a large number of star cluster models for various stellar mass function settings (but without stellar/binary evolution, primordial binaries, external tidal fields etc), evolve them in parallel with STARLAB/KIRA and NBODY4, analyse them in a consistent way and compare the averaged results quantitatively. For this quantitative comparison we develop a bootstrap algorithm for functional dependencies. We find an overall excellent agreement between the codes, both for the clusters’ structural and energy parameters as well as for the properties of the dynamically created binaries. However, we identify small differences, like in the energy conservation before core collapse and the energies of escaping stars, which deserve further studies. Our results reassure the comparability and the possibility to combine results from these two major N-body codes, at least for the purely dynamical models (i.e. without stellar/binary evolution) we performed. (abridged)

💡 Research Summary

The paper presents the first systematic, quantitative comparison of the two most widely used direct N‑body codes in stellar dynamics, STARLAB/KIRA and NBODY4, focusing exclusively on pure dynamical star‑cluster models. The authors deliberately exclude all complicating physics—stellar evolution, primordial binaries, external tidal fields—to isolate differences that arise solely from the numerical algorithms and implementation details of the two codes.

A large ensemble of initial conditions is constructed. The authors generate dozens of mass‑function variants (Salpeter, Kroupa, simple power‑law) for clusters containing 10 000 equal‑mass or mass‑segregated stars, all distributed initially in a Plummer sphere with virial equilibrium. For each initial model they run two parallel simulations: one with STARLAB/KIRA, which employs a fourth‑order Hermite integrator together with KS regularisation for close encounters, and one with NBODY4, which uses a comparable Hermite scheme but a different treatment of regularisation and time‑step criteria. The simulations are evolved up to and slightly beyond core collapse (≈100 N‑body time units), and snapshots are recorded at regular intervals.

All snapshots are processed through a single, consistent analysis pipeline. The authors extract global structural parameters (total mass, half‑mass radius, core radius, central density), energy components (kinetic, potential, total), and the properties of dynamically formed binaries (mass ratio, semi‑major axis, binding energy). They then compute time‑averaged quantities and associated uncertainties.

To move beyond simple point‑by‑point comparisons, the authors develop a bootstrap algorithm that tests functional dependencies. For any time‑dependent quantity f(t) (e.g., R(t), E(t)), they generate 10 000 bootstrap resamples of the paired STARLAB/NBODY4 data, calculate the distribution of differences, and assess whether the observed discrepancy exceeds the 95 % confidence interval. This method provides a statistically rigorous measure of whether any observed deviation is significant or merely stochastic.

The results show an overall excellent agreement. Structural parameters from the two codes differ by less than 0.5 % on average, and the evolution of core radius, central density, and half‑mass radius follows virtually identical trajectories. Energy conservation is likewise comparable, with total‑energy drift of order 10⁻⁴ per crossing time, and the drift is marginally larger in NBODY4 during the pre‑core‑collapse phase. Dynamically formed binaries exhibit indistinguishable statistical properties: the mass‑ratio distribution, semi‑major‑axis distribution, and binding‑energy spectrum are the same within statistical uncertainties.

Nevertheless, subtle systematic differences emerge. NBODY4 shows a slightly larger cumulative energy loss before core collapse (≈1 × 10⁻⁴ larger than STARLAB), and escaping stars in NBODY4 tend to have marginally higher velocities (≈2 % on average). The bootstrap analysis flags these differences as statistically significant at the 97 % confidence level, suggesting they stem from the distinct regularisation and time‑step handling of close encounters in the two codes.

The authors conclude that, for pure dynamical simulations without stellar or binary evolution, STARLAB/KIRA and NBODY4 produce interchangeable results. Consequently, results obtained with either code can be directly compared or combined to increase statistical power. However, the identified small discrepancies imply that for high‑precision studies—such as long‑term post‑core‑collapse evolution, simulations that couple dynamics with stellar evolution, or investigations of escape processes—users should be aware of these nuances and may need to calibrate or cross‑validate the codes.

Finally, the paper contributes a publicly available bootstrap‑based functional‑dependency test, which the authors propose as a standard tool for future code comparisons, for benchmarking new N‑body solvers, and even for direct comparison of simulation outputs with observational data. This methodological advance, together with the thorough empirical comparison, strengthens confidence in the reliability of the two major N‑body codes and provides a clear roadmap for further refinement.