© The Authors 2026

1 Introduction

Globular clusters (GCs) are among the oldest and most compact stellar systems in the Universe, comprising roughly 10⁶ stars, and they are widely regarded as prominent fossils of the early cosmic epochs. Their ubiquitous presence across galaxies, together with the wealth of detailed observations available, makes GCs privileged probes of substructure formation at high redshift and of the subsequent tumultuous phases of galaxy assembly.

Modeling and interpreting the physical properties of these stellar systems remains a challenging task, as multiphysics processes–such as stellar and dynamical evolution–must be followed simultaneously over a Hubble time, from the small scales of individual stars to the large scales of the galactic environment. A particular difficulty arises from the fact that GCs are “collisional” systems in which frequent close encounters between stars play a major role in shaping their global properties, driving their evolution, and ultimately leading to their dissolution.

Direct summation, through direct integration of Newtonian gravitational forces, has proven to be a key method for investigating GC evolution (see review by Spurzem & Kamlah 2023). For example, the NBODY6 code (Aarseth 2003) has enabled studies of the mechanisms driving GC mass loss (Baumgardt & Makino 2003; Baumgardt & Sollima 2017), the formation of tidal tails (e.g., Küpper et al. 2010; Webb & Bovy 2019), the influence of time-varying external fields in galactic environments (Bianchini et al. 2015; Miholics et al. 2016; Bianchini et al. 2017; Moreno-Hilario et al. 2024; Webb et al. 2024), and processes such as binary star formation, the formation of stellar remnants, and the production of gravitational wave sources (e.g., Sippel & Hurley 2013; Arca Sedda et al. 2024).

Simulations of GCs with a realistic number of stars (≈10⁶) are now possible thanks to high-performance computing (HPC) facilities, which enable simultaneous modeling of dynamical evolution and stellar evolution in the presence of the external environment. Using the code NBODY6++GPU (Wang et al. 2015), one-to-one simulations of GCs with N = 10⁵–10⁶ stars have been developed. Notable examples include the simulation presented in Sippel & Hurley (2013), the simulation of M4 (Heggie 2014), the suite of simulations DRAGON-I (Wang et al. 2016), DRAGON-II (Arca Sedda et al. 2024), and those from Kamlah et al. 2022a. These remarkable works have shown that the direct study of the long-term evolution of GCs remains challenging, with a realistic simulation still requiring computational times of more than ∼1 year even on parallel GPU architectures. However, alternative numerical methods relaxing the assumption of direct summation have proven to be extremely competitive in the task of simulating GCs with 10⁶ stars, particularly the code PeTar (Wang et al. 2020) and Monte Carlo method codes (the MOCCA code, see, e.g., Hypki & Giersz 2013; Giersz et al. 2013, and the CMC code, Rodriguez et al. 2022).

At the same time, observational facilities have advanced dramatically over the past ∼10 years, with high-quality data from both space- and ground-based telescopes providing precise kinematic measurements. These observations have revealed the unambiguous presence of internal rotation and velocity anisotropy in GCs, offering new constraints for theoretical models. Rotation was first measured in a few selected GCs thanks to line-of-sight observations (Lane et al. 2011; Bellazzini et al. 2012; Bianchini et al. 2013; Kacharov et al. 2014; Fabricius et al. 2014; Lardo et al. 2015; Cordero et al. 2017; Ferraro et al. 2018; Lanzoni et al. 2018) and proper motion data (van Leeuwen et al. 2000; Anderson & King 2003; Bellini et al. 2017; Häberle et al. 2024). These types of measurements became feasible for an extended number of GCs thanks to large spectroscopic surveys (Kamann et al. 2018; Szigeti et al. 2021; Martens et al. 2023; Petralia et al. 2024 and all-sky astrometric data from the Gaia mission (Bianchini et al. 2018a; Sollima et al. 2019; Vasiliev 2019; Vasiliev & Baumgardt 2021). To date, the combination of these extensive datasets gives us access to the 3D rotational information for ∼30 GCs (e.g., Leitinger et al. 2025 for a recent overview). The remarkable precision of these data has permitted measurements of rotation signals below the kilometer per second regime, of rotation differences between multiple stellar populations (Cordoni et al. 2020; Martens et al. 2023; Dalessandro et al. 2024; Leitinger et al. 2025), and between stars with different masses (Scalco et al. 2023). Velocity anisotropy has also been characterized on the plane of the sky for many GCs, showing that a variety of flavors are possible (from isotropy to mild radial and tangential anisotropy), using HST proper motions (Watkins et al. 2015; Libralato et al. 2022) and Gaia data (e.g., Jindal et al. 2019; Vasiliev & Baumgardt 2021; Cordoni et al. 2025). The newly released high-precision astrometry and photometry from the Euclid mission seem particularly promising for carrying out further studies of this kind, particularly when combined with Gaia data (see, e.g., Libralato et al. 2024).

On the modeling side, early efforts to understand rotation in GCs employed both direct N-body simulations and Fokker-Planck methods. Pioneering N-body studies explored the dynamical effects of rotation in self-consistent cluster models (Akiyama & Sugimoto 1989; Ernst et al. 2007; Kim et al. 2008; Fiestas & Spurzem 2010; Hong et al. 2013). At the same time, Fokker-Planck approaches (Einsel & Spurzem 1999; Kim et al. 2002; Fiestas et al. 2006; Tep et al. 2024) provided a complementary view of the long-term evolution, allowing researchers to track the influence of rotation on cluster structure and kinematics over extended timescales. Together, these early efforts set the stage for modern simulations that aim to capture the full complexity of rotating GCs.

Recent N-body simulations of rotating GCs have focused on the exploration of an increasing number of important physical ingredients, namely, the interplay between rotation and the tidal field (Tiongco et al. 2018, 2022), the effect of a stellar mass function (Tiongco et al. 2021; Livernois et al. 2022), the impact of stellar evolution (Kamlah et al. 2022b), the presence of multiple stellar populations (Berczik et al. 2025; White et al. 2026), and rotation in low-N clusters (Bissekenov et al. 2025). These simulations are limited because either they have a relatively low number of particles, <10⁵, or they have only been evolved for a fraction of the Hubble time. Today, only one large rotating N-body model with 10⁶ stars, stellar evolution, and an external tidal field has been produced (Wang et al. 2016, DRAGON-I).

In this work, we explicitly aim to fill this gap by simulating realistic GCs with internal rotation for a Hubble time and for a wide set of initial conditions. The goal is to provide a description of present-day Milky Way (MW) GCs and to connect their current properties to their properties in the early Universe, which is now accessible via JWST high-redshift observations (Vanzella et al. 2023; Adamo et al. 2024; Mowla et al. 2024). In this paper, we present a suite of 25 rotating GC simulations with N=250k–1.5M stars run with the code NBODY6++GPU, which we name the ROLLIN’ simulations (ROtating globular cLusters Long-term evolutioN simulations). In particular, we focus on the characterization of their long-term kinematic evolution as an important step for interpreting and building a modeling benchmark for real GCs.

The paper is structured as follows. In Sect. 2, we describe the simulation setup, the parameter space explored, and the initial and final GC global properties. Section 3 is devoted to the early and long-term evolution of the clusters, with a particular emphasis on kinematic properties. In Sect. 4, we connect the present-day GC properties to their primordial properties, and in Sect. 5 we perform a comparison with available kinematic observations. Finally, in Sect. 6 we present our summary and conclusions.

2 N-body simulations

We performed the simulations of rotating GCs using the code NBODY6++GPU, a GPU-accelerated direct N-body code (Wang et al. 2015; Kamlah et al. 2022a) based on the family of N-body codes originally developed by Sverre Aarseth (Aarseth 1985; Spurzem 1999; Aarseth 1999, 2003; Nitadori & Aarseth 2012). The code allows the Newtonian equation of motion for an N-body system to be directly integrated and simultaneously includes a series of astrophysical ingredients, such as stellar and binary evolution and the interaction with an external tidal field.

2.1 Initial conditions

As initial conditions, we adopted the family of distribution function–based models developed by Varri & Bertin (2012) and constructed to describe quasi-relaxed stellar systems with realistic rotation. These models feature axisymmetry, differential rotation, and a velocity anisotropy profile that is isotropic at the center, radially anisotropic in the intermediate regions, and mildly tangentially anisotropic in the outskirts.

The models are expressed as functions of the integral of motion I = I(E, J_z), with E and J_z as the energy and z-component of the angular momentum, defined as $$I=E-\frac{\omega J_z}{1+b{J}_z^{2c}}.$$(1)

The models are characterized by four dimensionless parameters:

• the concentration parameter, W₀, defined as the central depth of the dimensionless potential well;

• the rotation strength parameter, $\widehat{\omega }$, controlling the amount of rotation present in the cluster;

• the parameters b and c defining the shape of the rotation profile.

For a full description of each parameter, we refer to Varri & Bertin 2012¹.

To construct initial conditions representative of real GCs, we based our choice of the four dimensionless parameters on the model–observation comparisons presented in Bianchini et al. (2013) and Kacharov et al. (2014). In these works, the typical structural parameters derived from the best fit models are W₀ = 7, b = 0.01, and c = 1². While these parameters are kept fixed, we further explore three different values of the rotation parameter, $\widehat{\omega }=0.3,0.2,0.1$, which correspond to progressively decreasing rotation strengths. Models with higher $\widehat{\omega }$ are characterized by a stronger rotational support, quantifiable as a higher V/σ parameter, as indicated in Table 1.

The initial conditions are characterized by three values of initial number of stars, N = (2.5, 5, 15) × 10⁵, and two values for the initial half-mass radius r_50% = 1, 2 and 4 pc². The initial mass function (IMF) is assigned following a Kroupa (2001) mass function with minimum and maximum stellar mass of 0.08 and 100 M_⊙ and with no primordial binaries. Previous simulations of rotating GCs often consider an initial mass function truncated at lower stellar masses (e.g., Livernois et al. 2022 use a Kroupa mass function with a mass range of 0.1–1 M_⊙). To allow for an exploration of the effects due to the presence or absence of high mass stars, we set an additional simulation with maximum stellar mass of 50 M_⊙. The stellar metallicity is set to [Fe/H]=−1.3 (i.e., Z=0.001). All clusters are evolved in circular orbits with distances of d = 5, 10, 25, 50 kpc, around an idealized point-mass galaxy with mass M = 10¹¹ M_⊙. The rotation axes of the GCs are initialized along the z-axis, which is perpendicular to the orbital plane x-y around the point-mass galaxy. The majority of simulations are set with prograde internal rotation with respect to the orbital motion of the GC in the galaxy (i.e., the internal angular momentum vector is aligned with the orbital angular momentum vector), while three simulations are set with retrograde rotation (i.e., rotation axis along the negative direction of the z-axis). These orbital properties set all clusters in the tidal underfilling regime, with filling factors (i.e., the ratio between the half-mass ratio and the Jacobi radius, r_j³) of r_50%/r_j ∼ 0.01–0.05.

The implementation of stellar evolution for single and binary stars follows the prescriptions defined in Wang et al. (2016), based on the original SSE and BSE prescriptions of Hurley et al. (2000, 2002), commonly referred in NBODY6++GPU as level A (we refer to Kamlah et al. 2022a for details on the implementation of stellar evolution). In this context, for most of the simulations, the natal kick experienced by neutron stars at formation is set as a Maxwellian distribution with a dispersion of σ_nk = 265 km/s (see Hobbs et al. 2005; but also see Disberg & Mandel 2025). In order to explore different stellar remnant retentions, we employ for three simulations a natal kick dispersion of σ_nk = 30 km/s (low natal kick), and for two simulations σ_nk = 5 km/s (very low natal kick). For all simulations the natal kick of black holes is modeled including a fallback mechanism (see Belczynski et al. 2002). To allow for a comparison with the most updated version of the stellar evolution code, we ran one further simulation with 500k stars using the level C stellar evolution implementation⁴. The effects of the natal kick and stellar evolution prescriptions, and in particular their connection to the number of stellar remnants retained by the GCs, will be examined in a follow-up paper.

Finally, we run two control simulations with a different rotating model as initial condition. Instead of using the distribution function of Varri & Bertin (2012), we set a spherical, isotropic and nonrotating Wilson (1975) model with W₀ = 6, in which solid body rotation is introduced using the “Lynden-Bell demon” (Lynden-Bell 1962). This method corresponds to imposing positive values to the v_ϕ velocity component for a fraction of stars that otherwise had negative values (see, e.g., recent applications Rozier et al. 2019; Breen et al. 2021). In our case, we randomly select 50% of the stars in our sampled distribution function and impose the absolute value of v_ϕ and consider one prograde and one retrograde initial condition.

The naming convention used to identify the simulations, as reported in Table 1, explicitly describes the initial conditions. The format is the following: number of stars N, model type (indicated as -A, -B, or -C for the Varri & Bertin 2012 models with rotation strengths of $\widehat{\omega }=0.3,0.2,0.1$ respectively, and -W6 for the rotating Wilson models), half-mass radius in parsec, size of the circular orbit in kpc, plus complementary information if a property differs from the standard one. This additional information includes the use of level C stellar evolution (indicated as -lC), low or very low natal kicks of stellar remnants (indicated as -lk or -vlk), truncation of the IMF to 50 M_⊙ (indicated as -imf50), and retrograde internal rotation with respect to the galactic orbital motion (-retr). As an example, 1.5M-A-R4-10 designates the simulation with 1.5M stars, rotating Varri & Bertin (2012) model with $\widehat{\omega }=0.3$, half-mass radius r_50% = 4 pc and circular orbit with d = 10 kpc.

The suite of simulations consists of a total of 25 GC models evolved up to 14 Gyr, with the exception of the simulation 1.5M-A-R4-10 which is run until 12 Gyr. The simulations were run on the supercomputer Jean-Zay (GENCI–IDRIS) on NVIDIA V100 GPUs over a ∼3 yr period, under the time allocations Grand Challenge-101470 (08/2020), A10-A0100412451 (05/2021), and A13-A0130412451 (11/2022). The entire set of simulations was run for a total of ∼350 000 GPU hours. A typical simulation with 250k stars took approximately 20–50 days with 4–16 GPUs, while the largest 1.5M star simulation took ∼400 days with 12–48 GPUs. The total CO₂ equivalent consumption⁵ for our simulations is ∼9 tons of CO_2,eq.

Fig. 1 Time evolution of the mass; half-light radius, r50%; and filling factor, r50%/rj, 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.

2.2 General properties

The goal of this work is to investigate the long-term evolution of the kinematic and structural properties of rotating GCs, and their interplay with tidal effects and initial GCs’ properties. Our grid of simulations was specifically selected to explore these effects and these dependencies. Figure 1 shows the evolution of the mass, half-mass radius, and filling factor of the ROLLIN’ simulations compared to the observed values of MW GCs (from the Baumgardt Galactic GCs Database⁶). The snapshots at 12 Gyr of our simulations cover the parameter space of real MW GCs in the low-density regime. In essence, they allowed us to explore the following:

• a variety of initial mass surface densities ranging from Σ_∗ ≈ 10³ to Σ_∗ ≈ 10⁴ M_⊙/pc²;

• different strengths of the tidal field, including the effects of prograde versus retrograde cluster rotation with respect to their orbits in the galactic potential;

• different stellar evolution prescriptions (level A and level C) and the effects of different natal kicks.

In Tables 1 and 2, we report the summary of the initial and 12 Gyr properties, respectively. In particular, we include the estimate of internal rotation quantified by the peak of the rotation profile V_peak, the V_peak/σ₀ parameter with σ₀ the central value of the 1D velocity dispersion, and the β anisotropy parameter defined as $\beta =1-\left({\sigma }_{\theta }^2+{\sigma }_{\varphi }^2\right)/2{\sigma }_r^2$, with (r, θ, ϕ) representing the components in spherical coordinates of the velocity dispersion. Note that, not reported in Table 1, the initial values of β within the half-light radius for all the simulations are β ≈ 0, except for 250k-W6-R4-25 and 250k-W6-R4-25-retr which have β ≈ 0.08. The relaxation time at the half-mass radius t_rh is calculated using Eq. (2) in Bianchini et al. (2016).

In Figs. 2 and 3, we show a visualization of the 1.5M-A-R4-10 simulation as if it were observed with a JWST NIRcam setup. The color images were produced with three channels using the F070W, F115W, and F356W photometric filters, and using a filter-dependent PSF produced with the Python package WebbPSF (Perrin et al. 2012)⁷. The stellar magnitude of each star was calculated from the stellar parameters (luminosity, mass, stellar radius; contained in the outputs of the simulations) using the software FSPS (Conroy et al. 2009; Conroy & Gunn 2010)⁸ and employing BaSel spectra library and Padova isochrones library. Finally, the field of view was defined to have 132x132 arcsec⁹ size, with a pixel scale of 0.062 arcsec. Figure 2 displays three snapshots, at 0, 1, and 12 Gyr, observed at a distance of 200 kpc and with an inclination angle between the line-of-sight at the rotation axis of i=45^◦. These panels show the effect of dynamical evolution (e.g., expansion of the clusters) and of stellar evolution (e.g., change of color of the stellar population). Figure 3 shows instead the 12 Gyr snapshots observed face on (along the rotation axis, perpendicular to the orbital plane of the cluster) at different distances, underlying the amount of star-by-star details available in the dense central region as well as the large-scale effects of the tidal field (e.g., formation of tidal features). A video illustrating the time evolution of this simulation is also available at this link.

Fig. 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◦ and a field of view of 132 × 132 arcsec2. A video illustrating the time evolution of this simulation is available at this link.

Fig. 3 Visualization of the 1.5M-A-R4-10 model at 12 Gyr as if it were observed by the JWST NIRcam (with the F070W, F115W, and F356W filters) with a field of view of 132 × 132 arcsec2 at two different distances in an extragalactic and galactic context (d=800 kpc in the left image, and d=20 kpc in the right image). The GC is observed face on, along the direction of the rotation axis, which is perpendicular to the orbital plane.

3 Evolution of GC properties

In this section, we present the evolution of our simulated GCs. We first focus on their early evolution connected to the gravogyro-gravothermal collapse and then on the long-term evolution of their kinematic radial profiles. Finally, we discuss the evolution of global quantities as a function of relevant structural properties.

3.1 Early evolution

The evolution of GCs is strongly shaped by collisional dynamics causing energy exchanges between stars which, due to the negative heat capacity of self-gravitating stellar systems, lead to core collapse, a rapid increase of the central density of a cluster, which can be halted by the injection of energy due to binary formation (e.g., Hénon 1961; Lynden-Bell & Wood 1968; Meylan & Heggie 1997), in particular black hole binaries (Breen & Heggie 2013). Core-collapsed clusters are usually identified observationally via the presence of a central density cusp in their surface brightness profiles (e.g., Djorgovski & King 1986); but, more recently, identification methods based on kinematic data have also been put forward (e.g., Bianchini et al. 2018b; Libralato et al. 2018). However, in the early stages of GC evolution an increase of the central density can also take place due to the quick segregation of massive objects, including BHs (see, e.g., Kamlah et al. 2022b); this process is referred to as the initial collapse of the core. From a general point of view, the collapse of the core in a GC can be thought of as the outcome of the gravothermal instability that collisional stellar systems experience. It has been extensively studied in GCs under various conditions (e.g., Breen et al. 2017; Pavlík et al. 2024) and it has been noted that the presence of internal rotation adds further features to the collapse of the core. This behavior has often been interpreted as due to an additional catastrophic process, known as gravogyro instability, which is caused by the transport and redistribution of angular momentum within a cluster (see Akiyama & Sugimoto 1989; Ernst et al. 2007; Kamlah et al. 2022b and references therein).

In order to study the onset of core collapse in our simulations, we analyze the evolution of Lagrangian radii, in particular the radii containing 50%, 10%, and 1% of the total mass, r_50%, r_10% and r_1%, respectively. The analysis carried out here assumes spherical symmetry, an assumption that is not strictly valid in these early phases of evolution where deviations from sphericity can be severe due to the onset of dynamical instabilities. This point will be discussed in the next paper of the series.

In the first row of Fig. 4, we show the evolution of the Lagrangian radii for three selected simulations with 500k stars, representative of different initial densities and rotation strengths. The first 2 models, 500k-A-R2-10 and 500k-A-R4-10, are characterized by the same rotational support (V_peak/σ₀ = 1.20) and different densities (11.5 and 2.9 M_⊙ pc⁻²). The third model, 500k-C-R4-10, is characterized by weaker rotation (V_peak/σ₀ = 0.71) and low density (2.8 M_⊙ pc⁻², comparable to 500k-A-R4-10 model). The figures for all simulations are reported in Appendix A.

Figure 4 shows the expansion of the intermediate and outer layers (10% and 50% Lagrangian radii) and the contraction of the central region within r_1%. The minimum value of r_1%, corresponding to the maximum value of the central density, is reached by all simulations within the first 500 Myr of evolution. Higher density clusters show an earlier and deeper collapse, whereas low rotating clusters show a weaker signature of collapse. These results are in agreement with what was shown by Kamlah et al. (2022b) for rotating star clusters with 10⁵ stars. After the time of core collapse, the r_1% expands, as a result of heating mechanisms, such as binary formation or stellar evolution mass loss.

The second row of Fig. 4 shows the evolution of the mean stellar mass contained in each Lagrangian radii. At time zero, the mean mass is the same throughout the radial extent of the clusters, indicating homogeneous non mass-segregated initial conditions. While time evolves, massive stars segregate in the center, and the mean stellar mass contained within r_1% reaches values of 4–6 M_⊙. These high values indicate the presence of massive objects, mostly BHs (and their progenitors) formed in the first few million years through SNe. Interestingly, the simulations with stronger rotations reach the maximum of mean mass around core collapse, while the simulations with lower rotation strengths reach a peak later on, beyond 1 Gyr of evolution. This trend indicates that rotation has a key role in establishing the level and the timescale of mass segregation that a cluster can reach dynamically, which can then influence the timescale in which BHs segregate in a GC center. Finally, we note that after peaking, the mean-mass declines, indicating that massive objects which sank in the center are dynamically heating and ejecting each other out of the cluster through repeated 2- and 3-body interactions (see also Breen & Heggie 2013). A much milder decline is instead observable for the mean mass within r_50%, mostly reflecting stellar mass loss due to stellar evolution on a larger time-scale.

Fig. 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.

Fig. 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ϕ(r), the Vϕ/σ(r) profile indicating the ratio between ordered and random motion, the z-component of the angular momentum per unit mass lz(r), and the anisotropy profile β(r). The time resolution between the plotted lines is ≈64 Myr.

3.2 Long-term evolution of kinematic profiles

Using a spherical coordinate system (r, θ ϕ), we computed the radial profiles of four kinematic quantities:

1. the rotation profile V_ϕ(r), calculated as the mean of the ϕ-component of the velocity vector;

2. the V_ϕ/σ(r) profile, with σ as the 1D velocity dispersion profile tracing the ratio between ordered and random motions;

3. the z-component of the angular momentum per unit mass l_z(r);

4. the anisotropy profile $\beta (r)=1-\left({\sigma }_{\theta }^2+{\sigma }_{\varphi }^2\right)/2{\sigma }_r^2$, with β > 0 indicating radial anisotropy.

All profiles were calculated by dividing a GC in spherical shells such that each bin contains an equal number of stars. Figure 5 reports the time evolution of such profiles for the 1.5M-A-R4-10 simulation. Each curve is color-coded with a gray scale from 0 to 12 Gyr (from back to white), and the initial profiles are indicated with orange lines and the final with blue lines.

By analyzing the differences between the initial and final profiles of Fig. 5, we can appreciate the effects of the long-term evolution of GCs. In particular, we note an expansion of the cluster and a loss of angular momentum, clearly visible in the V_ϕ(r), V_ϕ/σ(r), and l_z(r) profiles as a reduction of the respective peaks and a shift toward larger radii. The anisotropy profile β(r) shows that the cluster starts with isotropy in the center, mild radial anisotropy in the intermediate regions, and mild tangential anisotropy in the outskirts. The long-term evolution imprints a similar trend to the β(r) profile at 12 Gyr, but with an enhanced radial anisotropy that builds up with time (consequence of the radial expansion of the cluster). Moreover, we note that the evolution of the kinematic profiles is more prominent in the early times compared to the much smoother evolution occurring on a longer timescale. This is due to the changes induced by early mass loss driven by stellar evolution, most effective in the early epochs of evolution (see also Sect. 3.4 and Fig. 9).

These general features described for the 1.5M-A-R4-10 simulation also apply to the entire set of simulations. However, the large parameter space covered by the initial conditions does imprint a remarkable variety of properties. In Fig. 6 we report the kinematic profiles of all simulations at 12 Gyr, normalized by their respective half-mass radius. Since we are interested in unveiling the imprint of the initial conditions on the final state of a GC, we color code each profile according to their initial tidal filling factor, defined as the ratio between the half-mass radius and the Jacobi radius, r_50%/r_j. This quantity provides both an estimate of the initial density of a cluster and a measure of the intensity of the tidal effects that a cluster experiences due to the presence of the external environment (a higher filling factor means that a GC is more vulnerable to tidal effects).

The first panel of Fig. 6 shows the variety of rotation profiles V_ϕ(r) produced by the set of simulations. The three negative rotation profiles correspond to clusters with retrograde internal rotation with respect to their orbital motion around the galaxy (we refer the reader to Sect. 3.3 for a discussion of the differences between prograde and retrograde clusters). For each profile, we identified the rotation peak by performing a fit of the profiles with the following relation routinely used to model rotation profiles in GCs (see, e.g., Mackey et al. 2013; Kacharov et al. 2014; Bianchini et al. 2018a): $$V(r)=\frac{2V_{\text{peak}}}{R_{\text{peak}}}\frac{r}{1+{\left(r/R_{\text{peak}}\right)}^2},$$(2)

where V_peak and R_peak represent the peak of the rotation curve and its radial position.

In Fig. 6, we indicate the rotation peaks with red points and the average of their radial position with a red vertical line. The profiles are characterized by a wide range of rotation amplitudes (from 0.2 to 3 km/s) and, noticeably, all rotation peaks are located around the half-mass radius r_50% (as demonstrated by the red vertical line almost coinciding with r_50%). Interestingly, the clusters with retrograde rotation have flatter rotation profiles, with peaks beyond the half-mass radius (see also Sect. 3.3). The shapes and amplitudes resulting from the long-term evolution of our simulations are fully consistent with observations of today’s MW GCs (e.g., Bianchini et al. 2018a).

As in the case of the rotation profiles, the V_ϕ/σ(r) and l_z(r) profiles also show an analogous variety of shapes. We note here that the three types of kinematic profiles that trace rotation in GCs do not show a sharp dependence on the initial filling factor. However, GCs with higher initial filling factors tend to show lower rotation signals, consistent with a stronger loss of angular momentum driven by the tidal field. We discuss evidence of this trend further in Sect. 3.4.

Of all kinematic profiles in Fig. 6, the anisotropy profiles β(r) show the widest diversity of shapes. In particular, despite the fact that all simulations start with a similar and globally weak radial anisotropy, different anisotropy flavors build up over time. All profiles are characterized by isotropy in the central region, but in the intermediate and outer parts they can exhibit: (i) strong radial anisotropy, (ii) mild radial anisotropy with quasi-isotropy in the outskirts, (iii) mild radial anisotropy with tangentiality in the outskirts, or (iv) fully tangential anisotropy. All these different flavors of anisotropy show a clear correlation with the initial filling factor, as is evident from the color coding: GCs with initially higher filling factors develop more tangential anisotropic profiles, while GCs with more underfilling initial conditions develop stronger radial anisotropy.

From a theoretical perspective, dynamical evolution, combined with mass-loss driven by stellar evolution, can naturally imprint radial anisotropy, due to the ejection of stars from the core to the halo of a cluster, a phenomenon that is particularly efficient during the post-core collapse phase (Spitzer & Shapiro 1972; Giersz & Heggie 1996). In contrast, the presence of a tidal field can quench radial anisotropy, bringing the outer part of the system toward isotropy and tangential anisotropy (Giersz & Heggie 1997), due to the preferential stripping of stars in radial orbits (e.g., Takahashi & Lee 2000; Baumgardt & Makino 2003). These effects, due to the coupling of the tidal field, internal dynamics, and stellar evolution, have been explored in previous simulations clearly showing that GCs can develop tangential anisotropy if they experience strong tidal fields (Tiongco et al. 2016b; Zocchi et al. 2016; Bianchini et al. 2017). Our simulations confirm these results and extend their validity to the mass regime of real GCs with >10⁵ stars, while also including realistic stellar evolution prescriptions.

Before concluding this section, a particular case worth mentioning is the simulation 250k-A-R2-5 characterized at 12 Gyr by an extremely low rotation signal. This model starts as a compact and fast rotating GC (r_50% = 2 pc and V_peak = 8.31 km/s) and experiences a particularly strong mass loss of ∆M/M = 0.92 (see properties in Tables 1 and 2). This is the result of the interplay between the initial high concentration of the cluster, the strong tidal field that the model experiences (due to the proximity of its orbit to the center of the galaxy, d=5 kpc), and the escape of stellar remnants due to their natal kick. These three ingredients trigger a strong mass loss that can lead to an abrupt cluster dissolution (a behavior often referred as “jumping”; see, e.g., Contenta et al. 2015; Giersz et al. 2019). The cluster quickly expands, becomes strongly tidally filling and it loses substantial mass, such that its r_50% starts to decrease. This strong mass loss leads to a significant loss of angular momentum, erasing any original rotation features, which explains the almost constant rotation profile with V_ϕ ≈ 0 km/s within the half-mass radius and the slightly negative rotation in the outer parts. The rotation profile outside the half-mass radius has a solid-body pattern with angular velocity (ω = 28.03 km/s/pc) approximately equal to half of the angular velocity of the cluster orbital motion around the host galaxy (Ω = 58.66 km/s/pc). This is fully consistent with the results of Tiongco et al. (2016a), predicting that a tidal field can imprint an internal rotation signature that is partially synchronized with the orbital angular momentum $(\Omega \approx -\frac{1}{2}\omega )$. Finally, the same model shows a strong signature of tangential anisotropy already within the half-mass radius. This specific case is rather instructive, as it shows how the observation of a lack of rotation in a cluster (or of a mild solid body rotation), combined with a tangential anisotropy signature, can indicate that the stellar system has undergone significant mass loss and is close to disruption.

Fig. 6 Kinematic profiles, as in Fig. 5, for the entire set of simulations at 12 Gyr. The profiles have been normalized by the respective half-mass radii r50% (vertical dashed lines) and are color coded by the value of their initial tidal filling factor (r50%/rj) tracing the tidal influence on a GC. In the rotation profile, Vϕ(r), we indicate the rotation peaks with red dots and their average radial position with a vertical red line.

Fig. 7 Kinematic profiles at 12 Gyr of three simulations with prograde rotation (with respect to the orbital angular momentum) compared to their retrograde counterparts. For the retrograde profiles, we plot the absolute values. Each profile has been normalized by its corresponding Jacobi radius, rj. Retrograde GCs are indicated as solid lines, while prograde ones are shown as dotted lines. Although prograde and retrograde simulations pairs have identical initial structures and rotation strengths, they significantly differ after 12 Gyr of evolution due to the complex coupling of internal rotation and the tidal field. Retrograde GCs retain a higher amount of angular momentum and show significant differences in rotation both in their outskirts and in the intermediate regions (∼0.5 rj).

3.3 The case of retrograde rotation

Our suite of models contains three pairs of GCs with identical structural initial conditions but different orientations of the rotation axis. Specifically, three simulations have prograde internal rotation with respect to the orbital angular momentum (250k-A-R4-25, 250k-A-R4-10, and 250k-W6-R4-25) and three simulations have retrograde rotation (250k-A-R4-25-retr, 250k-A-R4-10-retr, and 250k-W6-R4-25-retr). As highlighted in previous works, the tidal field will act differently in the prograde and retrograde cases: stars in prograde orbits are energetically less stable and more prone to be stripped by the tidal field; this leaves an excess of stars in retrograde orbits in a given cluster (see, e.g., Keenan & Innanen 1975; Fukushige & Heggie 2000; Ernst et al. 2008; Tiongco et al. 2016a, 2017; Claydon et al. 2019). As a consequence, the amount of angular momentum lost by a cluster will depend on the interplay between internal rotation and tidal field, and it will ultimately shape the rotation curves of present-day GCs.

Figure 7 shows the kinematic profiles of these simulations, with dashed lines indicating prograde GCs and solid lines retrograde ones. To facilitate the comparison, we plot the absolute values for the retrograde profiles. We normalized the profiles by the Jacobi radius r_j (the limiting radius of the GCs) to highlight whether the differences affect only the outer regions, where the tidal field is most effective, or also the inner and intermediate parts. We stress that for each pair all the kinematic profiles at t=0 Gyr are identical, so the differences in the plots can only be attributed to the prograde/retrograde configurations. Four main significant differences are noticeable:

1. Retrograde GCs are able to retain a higher amount of rotation. This is clearly visible by the overall higher values of rotation and the corresponding higher peaks in the profiles;

2. The V/σ and l_z profiles of retrograde GCs are characterized by a flat (or slightly increasing) shape in the intermediate/outer regions, whereas prograde GCs show a peak (around 0.2–0.4 r_j) followed by a decrease;

3. Prograde GCs show counter rotation signatures (with respect to the inner and intermediate regions) for radii >0.5–0.7 r_j, contrary to retrograde GCs where the rotation pattern keeps the same direction along the entire radial extent. These counter-rotating signatures appear in all prograde models, as clearly visible in Fig. 6 and are consistent with tidal field induced internal rotation, as described in Sect. 3.2 for the 250k-A-R2-5 model;

4. Differences in anisotropy between retrograde and prograde GCs are small, except in the outer regions, where retrograde GCs are more isotropic and can display tangentiality, specifically in the case of stronger tidal fields.

These results highlight the importance of the coupling between tidal field and internal rotation in shaping rotation signatures and indicate that the properties of kinematic profiles provide a powerful diagnostic to distinguish prograde from retrograde GCs. In particular, our work shows that the differences between prograde and retrograde GCs are already detectable in the intermediate regions of the cluster, at distances ∼0.5 r_j, rather than exclusively in the outermost parts, where observational analyses are complicated by low number statistics and high contamination from background and foreground stars. We note, however, that real clusters may be subject to more complex configurations between their orbital and internal rotation inclinations.

Fig. 8 Time evolution of the global properties of all our simulations: absolute value of the peak of the rotation profile Vpeak, |Vpeak|/σ0 parameter with σ0 (the central velocity dispersion), absolute value of the z-component of the angular momentum per unit mass calculated within the half-mass radius, and anisotropy parameter β within the half-mass radius. The initial values are indicated with orange crosses, and the curves are color coded according to their initial tidal filling factor, as in Fig. 6.

3.4 Long-term evolution of global kinematic quantities

Next, we focus on characterizing the evolution of global kinematic properties, namely:

1. the absolute value of the peak of the rotation profile |V_peak|;

2. the |V_peak|/σ₀ parameter, with σ₀ as the central value of the 1D velocity dispersion profile;

3. the z-component of the specific angular momentum within the half-mass radius r_50%;

4. the velocity anisotropy within the half-mass radius r_50%, β_<50%.

Figure 8 presents the evolution of these four global properties as a function of age. Each simulation is color-coded according to its initial filling factor, as in Fig. 6. As already noticed for the kinematic profiles, models with a higher initial filling factors show a more rapid change of their kinematic properties in the long term, highlighting again the importance of the tidal field in shaping the evolution of these quantities. This is particularly evident for the velocity anisotropy evolution, as it was also the case for Fig. 6. However, the most rapid and significant changes occur at earlier times (≲1 Gyr), when stellar evolution processes are the most efficient and induce a mass loss of up to 30% to all GCs, as shown in Fig. 9. The mass loss rate due to stellar evolution alone is dominant up to an average mass loss of 0.42 (see horizontal line of Fig. 9): when the clusters reach this total mass loss value, the main mechanism driving mass loss is the tidal escape of stars¹⁰. The time at which tidal effects become dominant is different for each model, depending on the tidal field strength they experience, and it ranges from 800 Myr to ∼9.5 Gyr.

The difference in mass loss due to the initial tidal conditions of a cluster is noticeable in the right panel of Fig. 9, where mass loss is plotted against the relaxation state of a cluster n_rel, defined as the number of relaxation times a cluster has experienced, n_rel=t_age/t_rh (with t_rh representing the relaxation time of a given time snapshot). After a common initial evolution, GCs with similar n_rel (i.e., similar internal dynamical conditions) lose mass (and stars) more efficiently if they were born with a higher filling factor, specifically r_50%/r_t ≳ 0.035.

In the case of the evolution of the specific angular momentum (Fig. 8), we note that most of the simulations reach a minimum around the time of core collapse (see Figs. A.1 and A.2). This reflects the behavior of the spatial distribution of a cluster that reaches a maximum concentration (i.e., minimum radial extent of the core region) at core collapse. Finally, we note that the different natal kicks prescriptions used in our suite of models (see Sect. 2.1) have only a minor effect on the kinematic quantities analyzed here. Table 2 shows that simulations with different natal kicks experience similar mass loss at 12 Gyr, and the models characterized by lower kicks display marginally lower rotational support, marginally higher radial anisotropy and are more spatially extended. These differences can be attributed to the different number of stellar remnants retained by a cluster depending on the natal kick prescription. We will analyze this topic in more depth in a follow up paper.

In the last part of this section, we focus on two main kinematic properties, |V_peak| and β_<50%, to identify the physical ingredients that primarily shape their evolution. In Figs. 10 and 11, we plot the rotation and anisotropy parameters against the half-mass relaxation time t_rh, n_rel, the total mass of the system, the mass loss, the filling factor, and the cumulative tidal strength experienced by a cluster¹¹. All of these quantities are measured for every snapshot. In each plot, the evolutionary tracks are shown as gray lines and the initial and 12 Gyr snapshots with orange and blue symbols, respectively. For the data points at 12 Gyr we perform a Pearson and Spearman correlation test, and indicate the corresponding coefficients (r_p and r_s) in each panel.

In the case of rotation (Fig. 10), the strongest correlation was obtained for |V_peak| versus log(M/M_⊙), with r_p = 0.78 and r_s = 0.5. The corresponding linear fit of the data is plotted as a dotted line and gives the following relation: $$\frac{\left|V_{peak}\right|}{\left[\text{km}{\text{s}}^{-1}\right]}=1.63\log M/M_{\odot }-7.06$$(3)

The rotation-mass correlation was previously observed for MW GCs (e.g., Bellazzini et al. 2012; Kamann et al. 2018; Bianchini et al. 2018a; Leitinger et al. 2025), implying that more massive clusters rotate faster. This is often taken as the result of the dissipation of angular momentum (e.g., Einsel & Spurzem 1999; Tiongco et al. 2017), as more massive systems would lose angular momentum less efficiently due to their longer relaxation timescales. However, we notice that our simulations do not show a strong correlation between |V_peak| and n_rel (or t_rh), suggesting that the observed rotation-mass relation could be due to both the formation mechanism of GCs and evolutionary processes. This would imply that more massive clusters already rotate faster at formation, as is the case in our suite of simulations¹². Interestingly, by normalizing V_peak by the initial peak rotation value (i.e., factoring out the impact of the initial rotation) we obtain a milder V_peak versus M correlation (r_p = 0.49 and r_s = 0.35).

This supports the idea that the initial conditions have a key role in setting the rotation at later ages.

In Fig. 10, other milder (anti)correlations are present, in particular for |V_peak| versus filling factor (r_p = −0.49 and r_s = −0.62) and for |V_peak| versus mass loss ∆M/M (r_p = −0.48 and r_s = −0.37). They indicate that also the tidal field has an impact in shaping rotation (a closer look at this is presented in Sect. 4.1).

In the case of velocity anisotropy (Fig. 11), the strongest correlation is found between anisotropy and mass loss (r_p = −0.85 and r_s = −0.71), but also other comparable correlations are found between anisotropy and relaxation conditions (t_rh and n_rel), the tidal field (filling factor and tidal strengths), and, to a lesser extent, the total mass. This suggests that both internal processes driven by two-body relaxation, and external processes, driven by the interaction with the host galaxy, contribute to shaping the velocity anisotropy of a cluster. The strong correlation obtained for the mass loss ∆M/M, indicates that this parameter efficiently captures these complex physical processes, as it traces the stellar escapes due to the interplay between two-body relaxation and tidal-field stripping¹³. A linear fit gives the following relation: $${\beta }_{<50\text{\%}}=-0.56\text{Δ}M/M+\mathrm{0.42.}$$(4)

Interestingly, all simulations reach a maximum value of radial anisotropy around a mass loss of 40%.

The general validity of the anisotropy and mass-loss relation is also supported by previous studies. In particular, Bianchini et al. (2017) found a similar relation for a set of nonrotating simulations, with a smaller number of stars (N=50k), evolved in a large variety of time-dependent tidal fields. Their work indicates that the transition between radially and tangentially anisotropic GCs occurs at mass loss ≈60% (see their Fig. 4). In our case, only two simulations have tangential anisotropy (β_50% < 0), but the linear fit suggests that this transition happens in a similar mass-loss regime, around 60–80%. Finally, the observational study of Watkins et al. (2015) detected a correlation between anisotropy and relaxation time, with a transition between radial and tangential anisotropic clusters at t_rh = 10⁹ yr. This is fully consistent with the relation found in the first panel of Fig. 11 and with Fig. 2 of Bianchini et al. (2017).

In summary, our simulations confirm that the evolution of the kinematic properties, in particular the velocity anisotropy, is driven by mass loss, which in turn is due to the combination of internal and external dynamical processes. For internal rotation, the results indicate that the values observed at present are also dependent on the amount of rotation at formation. We discuss this point further in Sect. 4.1, where we directly link the present-day rotation to the primordial rotation.

Fig. 9 Left panel: globular cluster mass loss ∆M/M vs. age. Right panel: globular cluster mass loss ∆M/M vs. number of relaxation times experienced by the GC, nrel = age/trh. Each GC curve is color coded according to the initial tidal filling factor. The red lines indicate the mass loss due to stellar evolution alone, common for all simulations and dominating in the first gigayear of evolution. The horizontal dashed line indicates the point where the mass loss starts to be dominated by tidal effects (stellar escapers) rather than stellar evolution processes. The plot shows the importance of the initial filling condition in defining the amount of mass loss cumulatively experienced by a GC after 12 Gyr of evolution. The nrel plot clearly shows that GCs formed with initial filling factors higher than 0.035 lose significantly more mass given the same relaxation condition.

Fig. 10 Correlation between the peak of the velocity profiles |Vpeak| of all simulations versus various quantities: the relaxation time trh at a given snapshot; the number of relaxation times nrel; the total mass M; the mass loss ∆M/M; the tidal filling factor r50%/rt; and the cumulative tidal strength experienced by a cluster (defined in Sect. 3.4). Orange crosses indicate the initial conditions, blue circles the snapshots at 12 Gyr, and the gray lines the evolutionary tracks. We highlight in black the track corresponding to the simulation 1.5M-A-R4-10. For each panel we performed a Pearson and a Spearman correlation test for the data points at 12 Gyr, and we indicate the corresponding coefficients (rp and rs), and plot the corresponding linear fit (dotted blue line) for the panel showing the strongest correlation, |Vpeak| vs. M.

Fig. 11 Correlation between the anisotropy β<50% calculated within the half-mass radius and the same quantities as in Fig. 10. The strongest correlation at 12 Gyr was found for β<50% vs. ∆M/M.

Fig. 12 Evolution of the absolute value of the rotation peak Vpeak normalized by the initial value at time t=0, Vpeak,t=0 as a function of mass loss ∆M/M. All simulations follow a similar track with a small variation, indicating that mass loss is the key parameter driving the decrease of rotation strength. The orange line is a fit of a simple function able to reproduce the shape of this evolutionary track. By measuring Vpeak and mass loss at the present time, one can estimate the initial rotation peak (and the rotation peak at any given time). The gray shaded area indicates the phases where mass loss is dominated by stellar evolution or by tidal stripping (see also Fig. 9).

4 Connecting present-day and primordial GC properties

In order to connect GCs present-day properties with their properties at formation, we exploit the evolutionary tracks derived from our suite of models. In particular, we focus on reconstructing the amount of primordial internal rotation of GCs and carry out a comparison with the structural parameters of MW GCs and high-z proto-GCs.

4.1 Reconstructing GC primordial rotation

The correlation found in Sect. 3.4 between rotation strength and current mass indicates that both the initial rotation and evolutionary processes cooperate in shaping this relation: GCs originally formed more massive and with stronger rotation would retain a stronger rotation after several Gyr of evolution.

To factor out the effects of initial rotation, we plot in Fig. 12 the evolution of the normalized rotation peak |V_peak/V_peak,t=0| (where V_peak,t=0 corresponds to the peak of the rotation profile at time t=0 Gyr) as a function of mass loss ∆M/M. All models follow very similar evolutionary tracks, indicating that mass loss is responsible for the decrease of rotational support. A similar result was highlighted in Tiongco et al. (2017) for three smaller rotating simulations.

Figure 12 indicates that by measuring the peak of rotation at the present time and the mass loss experienced by a cluster, it is possible to reconstruct the rotation value at earlier times. For this reason, we performed a fit of these evolutionary tracks using the following equation: $$\left|\frac{V_{\text{peak}}}{V_{\text{peak},t=0}}\right|=b\exp \left(-\frac{1}{2}\frac{{(\text{Δ}M/M)}^2}{a^2}\right),$$(5)

obtaining best-fit values of a and b of a = 0.26 and b = 0.87. Note that the formula is not physically motivated, but is found to provide a good description of all simulations throughout their entire long-term evolution. The plot also shows the regime where the mass loss rate of a cluster is dominated by stellar evolution processes or by stellar escapers: this transition happens when the models have lost on average a total of 42% of their total mass (see also Sect. 3.4 and Fig. 9).

Figure 12 demonstrates the strong impact of the long-term evolution in shaping the rotation strength. For a GC that underwent a mass loss of 50% after 12 Gyr of evolution (as is the typical case of our suite of simulations), the rotation amplitude has decreased by ≈84%. The rotation decreases by ≈75% if the GC has experienced 40% of mass loss. This implies that mass loss due to stellar evolution alone (∆M/M ≈ 40% throughout a cluster’s life; see Sect. 3.4 and Fig. 9) is responsible for a reduction in rotation strength of at least 75%. We can conclude that the rotation amplitudes observed in present-day GCs are a small fraction of their primordial rotation and Eq. (5) offers a means to quantify it.

4.2 Size and density of GCs at high redshift

Figure 13 shows the evolutionary tracks of our models in the surface density versus half-mass radius plane compared with current observations of real GCs. Specifically, we plot MW GCs (from the Baumgardt Galactic GCs Database¹⁴) and high-z observations of proto-GCs (Cosmic Gems, Adamo et al. 2024; Sunrise, Vanzella et al. 2023; Firefly Sparkle, Mowla et al. 2024). The evolutionary tracks of our simulations illustrate the effect of long-term evolution on GC structure: after 12 Gyr of evolution, the clusters have significantly expanded and reduced their surface density. Our models occupy the same parameter space as real GCs, suggesting that the observed high-z clusters could follow similar evolutionary tracks.

The right plot of Fig. 13 shows a normalized version of the left panel: despite the different initial conditions and tidal histories of each simulation, the evolutionary tracks are very similar. After 12 Gyr of evolution, the surface density of a GC has decreased by approximately 2 orders of magnitude (a factor of ≈1–6 × 10⁻²), while the half-mass radius has increased by a factor of ≈3–5. GCs with a higher initial filling factor have lost significantly more mass (see Fig. 9), and therefore they tend to have lower surface densities at a given half-mass radius. At 12 Gyr our models span a range in surface densities and sizes, which is the result of the specific initial condition and tidal history of each model. The only simulation experiencing a decrease in half-mass radius, at an approximately constant surface density, is the model 250k-A-R2-5, characterized by strong mass loss (as discussed in Sect. 3.2).

It is important to note that all our simulations are evolved in a simplified tidal field, specifically a circular orbit around a point-mass galaxy. More complex orbits within a realistic galactic potential would lead to a mass-loss rate that strongly depends on the orbit and local density. This would produce evolutionary tracks in the surface density versus half-mass radius plane that are nonlinear and could significantly enhance mass loss, particularly in the presence of strong tidal shocks—for example, those induced by disk crossings. Additional uncertainties may be due to the peculiar early evolution of high-z GCs due to the complexity of the time evolving tidal field (e.g., environment in the vicinity of their clumpy formation, buildup of host galaxy). Furthermore, all our simulations are initialized with half-mass radii of 1–4 pc, and are designed to reproduce low-density MW GCs. Higher initial densities could also affect the evolutionary tracks by allowing clusters to remain more compact throughout their evolution. Nevertheless, our comparison observations-simulations represents a first step toward understanding the evolution of GCs from high redshift to the present day, using realistic number of particles and incorporating all relevant dynamical and stellar evolutionary processes. Future comparisons with new data (e.g., Claeyssens et al. 2026) will provide a critical benchmark for validating and improving theoretical models of GC evolution.

Fig. 13 Left panel: evolution of the surface density calculated within the half-mass radius, Σ∗,<50%, as a function of the half-mass radius compared to observation of the present-day MW GCs (gray dots) and high-redshift observations of proto-GCs (black symbols). Orange crosses and blue dots indicate the initial and 12 Gyr snapshots of the simulations. Right panel: normalized version of the left panel where each simulation is color coded according to its initial tidal filling factor. This plot shows that in this parameter space, GCs follow a similar track: after 12 Gyr of evolution, the surface density decreases by a factor of ≈1–6 × 10−2, and the half-mass radius increases by a factor of ≈3–5.

5 Initial comparison with kinematic data

In this last section we carry out a first comparison with kinematic data of MW GCs. Figure 14 shows the V_peak/σ₀ parameter (see definition in Sect. 2.2) as a function of the half-mass relaxation time, the total mass and the velocity anisotropy. For each of our simulations we plot the snapshot at 12 Gyr with black circles. We compare our data with observations of the V/σ parameter of MW GCs from Bianchini et al. (2018a) and Leitinger et al. (2025). The V/σ parameter in Leitinger et al. (2025) is estimated using the global rotation amplitude in 3D and the associated velocity dispersion both in the plane-of-the sky and along the line-of-sight. Bianchini et al. (2018a) uses the same definition as in our current work, the peak of the rotation profiles divided by the central velocity dispersion, but their values refer only to proper motions in a reference system projected on the plane of the sky. Therefore, it should be emphasized that projection effects and differences in the kinematic definitions may introduce discrepancies between models and observations. It should also be remembered that our simulations are restricted to the parameter space of lower-density MW GCs, and therefore they do not encompass the full range of observed relaxation times and masses. For this reason, we regard the present comparison with the data as an initial step.

Nevertheless, the simulations and the data exhibit broadly consistent trends, with more massive clusters (i.e., clusters with longer relaxation time) showing higher V/σ. In the case of velocity anisotropy, intrinsic estimates of the β_<50% parameter are not currently available and estimates in the literature typically only refer to the ratio of the projected tangential and radial components of proper motion velocity dispersion. These anisotropy estimates depend in a nontrivial manner on projection effects and on the inclination angle of the rotation axis. We therefore refrain from comparing the anisotropy with observations, as a meaningful assessment would require projecting our simulations into the observational plane, which lies beyond the scope of this study. However, we note that our simulations predict a clear relation between V/σ and anisotropy. This correlation yields Pearson and Spearman coefficients of r_p = 0.68 and r_s = 0.51, respectively, corresponding to a strong linear association and a moderate monotonic trend. The only noticeable outlier of this correlation (empty circle) is the retrograde model 250k-A-R4-10-retr which is characterized by retrograde rotation coupled with a strong tidal field. An in depth analysis of this rotation versus anisotropy correlation offers a valuable testbed for future observations and motivates further efforts toward a comprehensive 3D kinematic analyses of GCs.

Fig. 14 Comparison of our simulations (at 12 Gyr, black circles) with recent observations (Bianchini et al. 2018a; Leitinger et al. 2025). From left to right: Vpeak/σ0 vs. relaxation time at the half-mass radius, total mass, and anisotropy within the half-mass radius. The dotted lines are the linear fits to the simulation data, for which the Pearson and Spearman correlation coefficient are indicated in the bottom right. The strongest correlation was obtained for Vpeak/σ0 vs. β<50%, for which data are not available. In this panel, we also note that the simulation 250k-A-R4-10-retr, characterized by retrograde rotation and strong tidal field, is a clear outlier (empty circle).

6 Conclusions

We have introduced a new suite of 25 direct N-body simulations of rotating GCs, the ROLLIN’ simulations, specifically designed to explore the impact of internal rotation in GCs for a wide range of initial conditions. Our models are characterized by a realistic number of particles (250k-1.5M stars) and are modeled on a star-by-star basis, including stellar evolution and interaction with the external environment, using the direct summation code NBODY6++GPU. With this method, all relevant evolutionary ingredients, inducing collisional dynamics, are taken into consideration. In this paper, we investigated the evolution of the kinematic properties of our models. Our analysis has demonstrated the following:

1. The early evolution of GCs is affected by the presence of internal rotation. Stronger rotating GCs undergo earlier and deeper core collapse in the first hundred million years of evolution, a finding consistent with the recent work of Kamlah et al. (2022b). This evolutionary phase is particularly efficient in segregating massive stars and stellar remnants in the center of clusters. Fast rotating GCs reach a higher level of mass segregation early in their evolution. This also suggests the importance of rotation in setting the process of segregation of BHs in a GC center. Further studies to disentangle the combined effects of rotation and GC structural properties (e.g., their concentration) are needed to fully understand the impact on the retention of compact objects. Such investigations could provide key constraints for interpreting observations of gravitational wave sources originating from dense stellar environments;

2. Mass loss, originating from both stellar evolution and tidal stripping, is the main driver of the evolution of kinematic properties. While the clusters evolve, their angular momentum dissipates, and the rotation strength, quantified as the peak of the rotation profile, decreases. Independently of the initial conditions and evolution, all clusters show a rotation peak around the half-mass radius, r_50%, as well as shapes and amplitudes of the rotation profiles that are fully consistent with observations of today’s MW GCs. After 12 Gy of evolution a correlation between cluster rotation and cluster mass builds up, with more massive clusters retaining higher rotation signatures than less massive clusters, consistent with rotation measurements of MW GCs (e.g., Bianchini et al. 2018a; Leitinger et al. 2025);

3. The tidal history of a GC is another key factor in the evolution of rotation profiles. Clusters with prograde rotation (with respect to their orbital angular momentum) reduce their rotation strength at a quicker pace than clusters with retrograde rotation. This is due to the preferential stripping of stars in prograde orbits, which are energetically less stable. Prograde GCs are characterized by a V/σ profile decreasing toward the outer regions with signatures of counter rotation. Retrograde clusters are instead characterized by a more extended and pronounced rotation profile and a flat (or increasing) V/σ profile. The differences between retrograde and prograde GCs are already visible well within the Jacobi radius. The shape of present-day rotation profiles may therefore contain important insights into the orbital history of a GC;

4. The evolution of velocity anisotropy is driven by mass loss. All clusters are initially close to isotropy and develop radial anisotropy, peaking around 40% mass loss (from both stellar escapers and stellar evolution), before progressing toward isotropy or tangential anisotropy for increasing mass loss. The GCs that enter the tangential anisotropy regime have undergone mass loss of ∼60–80%. These systems are characterized by large initial filling factors (r_50%/r_j ≈ 0.05, in the underfilling regime), making them more sensitive to the external tidal field and resulting in higher mass loss rates. The combined signature of tangential anisotropy and the lack of internal rotation can be indicative of a cluster approaching disruption;

5. The present-day rotation of GCs is at least a factor of five lower than the rotation at formation. From the evolutionary tracks of our simulations, we derived an empirical relation (Eq. 5) that allowed us to reconstruct a GC’s primordial rotation based on its present rotation strength and its mass loss. As an example, a mass loss of 40% (i.e., the typical value of mass loss due to stellar evolution alone throughout a cluster’s life) implies a decrease of the rotation peak of ≈75%. This indicates that internal rotation must have been a key dynamical ingredient already at the earliest phases of the life of a GC;

6. Our simulations provide a way to compare the structural properties of old MW GCs with those of proto-GCs observed at high-z. The evolutionary tracks of our models reveal that the combined effects of stellar and dynamical evolution cause clusters to expand by a factor of about three to five. This expansion is coupled with a decrease in surface density of up to two orders of magnitude (a factor of ≈1– 6 × 10⁻²), which is additionally driven by mass loss. These findings suggest that the observed high-redshift GCs could evolve into clusters with properties resembling those seen in the MW today.

In summary, we have shown that internal rotation is a key ingredient in the evolution of GCs, playing a fundamental role even at the earliest stages of their formation. Reconstructing the role of rotation in the life of a GC is therefore essential for gaining insights into GC origin, which remains largely unknown, even in the context of formation and kinematic evolution of multiple stellar populations.

The qualitative agreement between the ROLLIN’ simulations and the available kinematic observations of GCs indicates that these models provide a valuable benchmark for interpreting present-day as well as future observations. However, a thorough comparison between models and observations will require either investigating the projected properties of our simulations or performing a comprehensive 3D analysis of the kinematics of MW GCs—which is now possible thanks to the combination of large astrometric and spectroscopic surveys.

Finally, while our simulations mark a significant advance toward state-of-the-art direct models with a realistic number of stars and broad parameter coverage, a few key improvements remain to be made. In particular these involve (i) the inclusion of primordial binaries in the initial conditions, which is the current computational bottleneck; (ii) the adoption of more realistic time-dependent tidal fields to place GCs within a proper galactic and cosmological context; (iii) the use of less idealized initial conditions capable of capturing the complexity of cluster formation arising from the intricate interplay between gas and stars; and (iv) the inclusion of multiple stellar populations known to play a primary role already in the earliest stages of GC formation. Significant progress in these areas is already occurring (e.g., Arca Sedda et al. 2024; Cournoyer-Cloutier et al. 2024; Karam & Sills 2024; Lacchin et al. 2024; Webb et al. 2024; Lahén et al. 2025b,a; Aros et al. 2025; Giersz et al. 2025; Berczik et al. 2025), thus paving the way for a deeper understanding of the origin of GCs and their relation to high-redshift star formation and to galactic evolution.

Data availability

The simulation outputs from the ROLLIN’ suite are available from the corresponding author upon reasonable request.

Acknowledgements

The set of simulations was run for a total of ≈350 000 GPU hours. This corresponds to a carbon footprint of ≈9 tons of CO_2,eq, estimated using the methodology developed by the Labos 1point5 collective¹⁵. The authors thank the referee for their very useful insights. PB gratefully acknowledges Hugo Spitz and Gabriel Bounias for their valuable contributions to the visualization of simulations during their master’s internships. PB would also like to thank Rodrigo Ibata and Christian Boily for insightful discussions, and Jonathan Chardin for contributing to the early development of this work. This project was provided with computing resources by GENCI at IDRIS thanks to the following time allocations on the supercomputer Jean Zay (V100 partition): Grand Challenge-101470, A10-A0100412451, and A13-A0130412451. The authors would also like to acknowledge the High Performance Computing Center of the University of Strasbourg for supporting this work by providing access to computing resources. Part of the computing resources were funded by the Equipex Equip@Meso project (Programme Investissements d’Avenir) and the CPER Alsacalcul/Big Data. PB and AM acknowledge financial support by the IdEx framework of the University of Strasbourg. ALV acknowledges support from a UKRI Future Leaders Fellowship (MR/S018859/1; MR/X011097/1). AA acknowledges support for this paper from project No. 2021/43/P/ST9/03167 co-funded by the Polish National Science Center (NCN) and the European Union Framework Programme for Research and Innovation Horizon 2020 under the Marie Skłodowska-Curie grant agreement No. 945339. This research was also funded in part by NCN grant number 2024/55/D/ST9/02585. For the purpose of Open Access, the authors have applied for a CC-BY public copyright license to any author Accepted Manuscript (AAM) version arising from this submission.

References

Appendix A Evolution of Lagrangian radii and mean enclosed mass

Evolution of Lagrangian radii (Fig. A.1) and mean enclosed mass (Fig. A.2) for all the models in our suite of simulations.

Fig. A.1 Evolution of 50%, 10% and 1% Lagrangian radii (orange, gray, and black lines respectively) for all models, as in Fig. 4. The time of core collapse is indicated with a vertical line; when the core collapse is not clearly definable, no line is plotted.

Fig. A.2 Evolution of the mean mass within the 50%, 10% and 1% Lagrangian radii (orange, gray, and black lines respectively) for all models, as in Fig. 4. The time of core collapse is indicated with a vertical line; when the core collapse is not clearly definable, no line is plotted.

All Tables

All Figures

Fig. 1 Time evolution of the mass; half-light radius, r50%; and filling factor, r50%/rj, 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.

In the text

Fig. 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◦ and a field of view of 132 × 132 arcsec2. A video illustrating the time evolution of this simulation is available at this link.

In the text

Fig. 3 Visualization of the 1.5M-A-R4-10 model at 12 Gyr as if it were observed by the JWST NIRcam (with the F070W, F115W, and F356W filters) with a field of view of 132 × 132 arcsec2 at two different distances in an extragalactic and galactic context (d=800 kpc in the left image, and d=20 kpc in the right image). The GC is observed face on, along the direction of the rotation axis, which is perpendicular to the orbital plane.

In the text

Fig. 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.

In the text

Fig. 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ϕ(r), the Vϕ/σ(r) profile indicating the ratio between ordered and random motion, the z-component of the angular momentum per unit mass lz(r), and the anisotropy profile β(r). The time resolution between the plotted lines is ≈64 Myr.

In the text

Fig. 6 Kinematic profiles, as in Fig. 5, for the entire set of simulations at 12 Gyr. The profiles have been normalized by the respective half-mass radii r50% (vertical dashed lines) and are color coded by the value of their initial tidal filling factor (r50%/rj) tracing the tidal influence on a GC. In the rotation profile, Vϕ(r), we indicate the rotation peaks with red dots and their average radial position with a vertical red line.

In the text

Fig. 7 Kinematic profiles at 12 Gyr of three simulations with prograde rotation (with respect to the orbital angular momentum) compared to their retrograde counterparts. For the retrograde profiles, we plot the absolute values. Each profile has been normalized by its corresponding Jacobi radius, rj. Retrograde GCs are indicated as solid lines, while prograde ones are shown as dotted lines. Although prograde and retrograde simulations pairs have identical initial structures and rotation strengths, they significantly differ after 12 Gyr of evolution due to the complex coupling of internal rotation and the tidal field. Retrograde GCs retain a higher amount of angular momentum and show significant differences in rotation both in their outskirts and in the intermediate regions (∼0.5 rj).

In the text

Fig. 8 Time evolution of the global properties of all our simulations: absolute value of the peak of the rotation profile Vpeak, |Vpeak|/σ0 parameter with σ0 (the central velocity dispersion), absolute value of the z-component of the angular momentum per unit mass calculated within the half-mass radius, and anisotropy parameter β within the half-mass radius. The initial values are indicated with orange crosses, and the curves are color coded according to their initial tidal filling factor, as in Fig. 6.

In the text

Fig. 9 Left panel: globular cluster mass loss ∆M/M vs. age. Right panel: globular cluster mass loss ∆M/M vs. number of relaxation times experienced by the GC, nrel = age/trh. Each GC curve is color coded according to the initial tidal filling factor. The red lines indicate the mass loss due to stellar evolution alone, common for all simulations and dominating in the first gigayear of evolution. The horizontal dashed line indicates the point where the mass loss starts to be dominated by tidal effects (stellar escapers) rather than stellar evolution processes. The plot shows the importance of the initial filling condition in defining the amount of mass loss cumulatively experienced by a GC after 12 Gyr of evolution. The nrel plot clearly shows that GCs formed with initial filling factors higher than 0.035 lose significantly more mass given the same relaxation condition.

In the text

Fig. 10 Correlation between the peak of the velocity profiles |Vpeak| of all simulations versus various quantities: the relaxation time trh at a given snapshot; the number of relaxation times nrel; the total mass M; the mass loss ∆M/M; the tidal filling factor r50%/rt; and the cumulative tidal strength experienced by a cluster (defined in Sect. 3.4). Orange crosses indicate the initial conditions, blue circles the snapshots at 12 Gyr, and the gray lines the evolutionary tracks. We highlight in black the track corresponding to the simulation 1.5M-A-R4-10. For each panel we performed a Pearson and a Spearman correlation test for the data points at 12 Gyr, and we indicate the corresponding coefficients (rp and rs), and plot the corresponding linear fit (dotted blue line) for the panel showing the strongest correlation, |Vpeak| vs. M.

In the text

	Fig. 11 Correlation between the anisotropy β<50% calculated within the half-mass radius and the same quantities as in Fig. 10. The strongest correlation at 12 Gyr was found for β<50% vs. ∆M/M.
In the text

Fig. 12 Evolution of the absolute value of the rotation peak Vpeak normalized by the initial value at time t=0, Vpeak,t=0 as a function of mass loss ∆M/M. All simulations follow a similar track with a small variation, indicating that mass loss is the key parameter driving the decrease of rotation strength. The orange line is a fit of a simple function able to reproduce the shape of this evolutionary track. By measuring Vpeak and mass loss at the present time, one can estimate the initial rotation peak (and the rotation peak at any given time). The gray shaded area indicates the phases where mass loss is dominated by stellar evolution or by tidal stripping (see also Fig. 9).

In the text

Fig. 13 Left panel: evolution of the surface density calculated within the half-mass radius, Σ∗,<50%, as a function of the half-mass radius compared to observation of the present-day MW GCs (gray dots) and high-redshift observations of proto-GCs (black symbols). Orange crosses and blue dots indicate the initial and 12 Gyr snapshots of the simulations. Right panel: normalized version of the left panel where each simulation is color coded according to its initial tidal filling factor. This plot shows that in this parameter space, GCs follow a similar track: after 12 Gyr of evolution, the surface density decreases by a factor of ≈1–6 × 10−2, and the half-mass radius increases by a factor of ≈3–5.

In the text

Fig. 14 Comparison of our simulations (at 12 Gyr, black circles) with recent observations (Bianchini et al. 2018a; Leitinger et al. 2025). From left to right: Vpeak/σ0 vs. relaxation time at the half-mass radius, total mass, and anisotropy within the half-mass radius. The dotted lines are the linear fits to the simulation data, for which the Pearson and Spearman correlation coefficient are indicated in the bottom right. The strongest correlation was obtained for Vpeak/σ0 vs. β<50%, for which data are not available. In this panel, we also note that the simulation 250k-A-R4-10-retr, characterized by retrograde rotation and strong tidal field, is a clear outlier (empty circle).

In the text

	Fig. A.1 Evolution of 50%, 10% and 1% Lagrangian radii (orange, gray, and black lines respectively) for all models, as in Fig. 4. The time of core collapse is indicated with a vertical line; when the core collapse is not clearly definable, no line is plotted.
In the text

	Fig. A.2 Evolution of the mean mass within the 50%, 10% and 1% Lagrangian radii (orange, gray, and black lines respectively) for all models, as in Fig. 4. The time of core collapse is indicated with a vertical line; when the core collapse is not clearly definable, no line is plotted.
In the text