N-body Simulations with Live Stellar Evolution

An N-body code containing live stellar evolution through combination of the software packages NBODY6 and STARS is presented. Operational details of the two codes are outlined and the changes that have been made to combine them discussed. We have computed the evolution of clusters of 10 000 stars using the combined code and we compare the results with those obtained using NBODY6 and the synthetic stellar evolution code SSE. We find that, providing the physics package within STARS is set up correctly to match the parameters of the models used to construct SSE, the results are very similar. This provides a good indication that the new code is working well. We also demonstrate how this physics can be changed simply in the new code with convective overshooting as an example. Similar changes in SSE would require considerable reworking of the model fits. We conclude by outlining proposed future development of the code to include more complete models of single stars and binary star systems.

💡 Research Summary

The paper presents a novel computational framework that integrates the direct‑N‑body dynamics code NBODY6 with the detailed stellar evolution code STARS, thereby enabling “live” stellar evolution within an N‑body simulation. NBODY6 is renowned for its accurate treatment of gravitational interactions, close encounters, and binary dynamics for systems containing up to several hundred thousand particles, but it traditionally relies on synthetic stellar evolution prescriptions such as the SSE (Synthetic Stellar Evolution) fitting formulae to supply stellar radii, luminosities, and mass‑loss rates. STARS, on the other hand, solves the full set of stellar structure equations, including nuclear burning, convection, overshooting, and wind mass loss, but it has not been coupled to a dynamical N‑body integrator.

To bridge this gap, the authors designed an interface layer that exchanges the essential stellar parameters (mass, radius, luminosity, effective temperature, and composition) between NBODY6 and STARS at each global time step. They re‑engineered data structures to be compatible across the C++ (NBODY6) and Fortran (STARS) codebases, eliminated intermediate file I/O, and employed shared memory to minimise communication overhead. A key technical challenge was synchronising the disparate time‑step requirements: NBODY6 typically uses dynamical steps of order 10⁴–10⁵ yr, whereas STARS may need sub‑year steps during rapid evolutionary phases (e.g., helium flash, thermal pulses). The solution introduced a “sub‑step” scheme in which STARS is called multiple times within a single NBODY6 step; the resulting stellar updates are averaged and fed back to the N‑body integrator, ensuring consistent forces and mass‑loss effects. Whenever a mass‑loss event occurs (e.g., stellar winds, supernova ejecta), the NBODY6 particle mass is instantly updated, preserving momentum conservation and allowing the dynamical response to be captured in real time.

The authors validated the combined code by simulating a 10 000‑star open cluster with a Kroupa initial mass function, metallicity Z = 0.02, and mixing‑length parameter α_MLT = 2.0. They evolved the system for 1 Gyr and compared global diagnostics—mass function evolution, colour‑magnitude diagram, core density profile, binary fraction—to a reference simulation that used NBODY6 together with the SSE fitting formulae. The two approaches produced statistically indistinguishable results, demonstrating that the live STARS evolution reproduces the synthetic prescriptions when the underlying physics (e.g., opacities, reaction rates) are matched.

A particularly illustrative experiment involved modifying the convective overshooting parameter (α_ov) within STARS. By increasing α_ov from 0.0 to 0.2, the main‑sequence lifetimes of intermediate‑mass stars lengthened, and the turn‑off point in the colour‑magnitude diagram shifted accordingly. Because the change was made directly in the stellar evolution module, the dynamical simulation immediately reflected the altered stellar radii and mass‑loss histories, leading to subtle differences in cluster expansion and binary interaction rates. Replicating the same effect in SSE would require re‑deriving the entire set of fitting formulae, underscoring the flexibility of the live‑evolution approach.

The paper also discusses current limitations and future development plans. At present, only single‑star evolution is supported; binary evolution, mass transfer, common‑envelope phases, and compact‑object formation are not yet integrated. The authors outline a roadmap to couple STARS’s binary evolution capabilities with NBODY6’s sophisticated treatment of binary dynamics, enabling realistic modelling of Roche‑lobe overflow, tidal interactions, and gravitational‑wave progenitors. Performance scaling is another focus: the team intends to implement GPU acceleration for the force calculations and MPI‑based parallelism for the stellar evolution sub‑steps, aiming to handle clusters of 10⁵–10⁶ stars.

In summary, the work demonstrates that a tightly coupled NBODY6‑STARS system can faithfully reproduce results obtained with traditional synthetic evolution while offering the ability to modify underlying stellar physics on the fly. This opens the door to a new generation of star‑cluster simulations where detailed stellar interiors, mass‑loss processes, and dynamical evolution are treated self‑consistently, providing a powerful tool for interpreting observations of young clusters, globular clusters, and the stellar populations of external galaxies.