ROLLIN': Rotating globular cluster simulations. I. The kinematic evolution of realistic direct N-body models

Abstract

Internal rotation has emerged as a fundamental feature of globular clusters (GCs), yet its origin and long-term evolution remain poorly understood. We explore the evolution of rotating GCs over a Hubble time under the combined influence of two-body relaxation, tidal field, and stellar evolution. We introduce the ROLLIN' simulations, a suite of 25 N-body models characterized by a realistic number of stars from 250k to 1.5M, ran with the direct N-body code NBODY6++GPU and evolved for 14 Gyr. With present-day masses of 5 x 10^4 - 5x10^5 M_sun, the models cover the parameter space of low-density MW GCs. Our analysis reveals that rapidly rotating GCs experience earlier and more pronounced core collapse, efficiently segregating massive objects and remnants in their centers within the first few 100 Myr. In the long-term, internal rotation declines and a correlation emerges between rotation and GC mass, in agreement with observations. The primary driver of this evolution is mass loss, capturing both internal (stellar evolution, evaporation) and external processes (tidal stripping). The velocity anisotropy also evolves in response to mass loss: GCs initially near isotropy develop radial anisotropy, peaking around 40% mass loss, before progressing toward isotropy or tangentiality. The GC orbital history also plays a role, as retrograde rotators retain rotation more effectively than prograde rotators. Finally, we quantify the long-term changes of GCs after 12 Gyr: (1) The surface density decreases by up to 2 orders of magnitude. (2) The half-mass radius increases by a factor of 3-5. (3) The rotation decreases by a factor >5 for GCs that have lost >50% of their mass. The ROLLIN' simulations demonstrate that angular momentum is crucial to understand the origin, evolution, and survival of GCs. These models provide a benchmark for interpreting GC observations in the local and high-z Universe.

Figures (15)

• Figure 1: Time evolution of the mass; half-light radius, $r_{50\%}$; and filling factor, $r_{50\%}/r_j$, of our set of simulations. Orange crosses are the initial conditions, blue dots are the snapshots at 12 Gyr, and the gray lines indicate the time evolution. The comparison with MW GCs (gray dots, Baumgardt GCs Database) highlights that our simulated GCs cover the current parameter space of real GCs in the low-density regime.
• Figure 2: Visualization of the 1.5M-A-R4-10 model as if it were observed by the JWST NIRcam with the F070W, F115W, and F356W filters. The three images correspond to snapshots at 0, 1, and 12 Gyr observed at a distance of d=200 kpc and with an inclination angle between the line of sight and the rotation axis of $i=45^\circ$ and a field of view of 132x132 arcsec$^2$. A video illustrating the time evolution of this simulation is https://www.youtube.com/watch?v=o_C2nwJq560.
• Figure 4: Top row: Evolution of 50%, 10%, and 1% Lagrangian radii (orange, gray, and black lines) for simulations 500k-A-R2-10, 500k-A-R4-10, and 500k-C-R4-10. The three simulations are representative of different densities and rotation strengths. The vertical lines indicate the time of core collapse. Higher density GCs show an earlier and deeper collapse, whereas low rotating GCs show a weaker signature of collapse. Bottom row: Evolution of the mean mass within the same Lagrangian radii of the top panel. The 1% Lagrangian radii mainly contain massive objects, corresponding to BHs (or their progenitors). High-density and strong rotating GCs mass segregate more efficiently, in timescales of a few hundred million years.
• Figure 5: Evolution of the kinematic profiles of the 1.5M-A-R4-10 simulation from time t=0 Gyr (orange lines) to t=12 Gyr (blue lines). The intermediate time snapshots are color coded in gray scale. The panels report the rotation profile $V_\phi$(r), the $V_\phi/\sigma$(r) profile indicating the ratio between ordered and random motion, the z-component of the angular momentum per unit mass $l_z(r)$, and the anisotropy profile $\beta(r)$. The time resolution between the plotted lines is $\approx64$ Myr.
• Figure 6: Kinematic profiles, as in Fig. \ref{['1.5Mprofiles']}, for the entire set of simulations at 12 Gyr. The profiles have been normalized by the respective half-mass radii $r_{50\%}$ (vertical dashed lines) and are color coded by the value of their initial tidal filling factor ($r_{50\%}/r_j$) tracing the tidal influence on a GC. In the rotation profile, $V_\phi$(r), we indicate the rotation peaks with red dots and their average radial position with a vertical red line.