© The Authors 2025

1 Introduction

The dynamics of stellar systems is a major research field of astrophysics. It focuses on the gravitational interactions of stars in star clusters and can provide crucial insights into stellar evolution and Galaxy formation (see, e.g., Chandrasekhar 2005; Binney & Tremaine 2011; Krumholz et al. 2019). It is also essential for the study of the constituents of stellar systems, such as binaries, black holes, and neutron stars, and to the characterization of their dynamical and observational effects on dense stellar systems (see, e.g., Sridhar & Touma 1999; Ivanova et al. 2006; Alexander 2017; Ishchenko et al. 2024; Berczik et al. 2025). The study of stellar dynamics can provide constraints on the conditions of star formation and can be used to examine the effects of properties on the system. One of such properties is bulk rotation, which is the rotational motion of a star cluster around its center of mass. Star clusters inherit their initial rotation from molecular clouds through the star formation process (Mapelli 2017; Lahén et al. 2020b). At later times, bulk rotation can evolve through two-body relaxation, tidal interactions with other star clusters or molecular clouds, and tidal interactions with the Galactic tides. Over time, bulk rotation can significantly influence the morphology, core collapse process, and long-term evolution of star clusters.

Numerous observations have confirmed the presence of bulk rotation in star clusters. Bianchini et al. (2018) observed 81 globular clusters and found clear evidence of rotation in 21 of them. Furthermore, Pang et al. (2021) found velocity anisotropy among the elongated clusters. They suggest that velocity anisotropy might be induced by global rotation. There is also evidence that globular clusters in Gaia Data Release 3, NGC 104 and NGC 2904, show clear evidence of a rotation-mass relation (Scalco et al. 2023). Jadhav et al. (2024) observed the proper motions and radial velocities of 1379 open clusters and identified rotation in ten clusters, and suggested 16 cluster candidates with rotation. Hao et al. (2024) discovered a 3D rotation signature in three open clusters, Pleiades, Alpha Persei, and Hyades, using Gaia DR3 proper motions and radial velocities.

Initial bulk rotation has been numerically studied using hydrodynamical simulations of star formation in molecular clouds. Mapelli (2017) found that initial bulk rotation originates from molecular clouds as a result of interactions with gas filaments and clumps. They indicated that rotation should be commonly noticeable in massive embedded star clusters with masses of M ≈ 4 × 10³ M_⊙. Lahén et al. (2020a) investigated star clusters formed in metal-poor starburst dwarf galaxies and spotted significant angular momentum upon formation in their simulations. They found that more massive clusters tend to have a greater specific angular momentum, and that the angular momentum is not always correlated with flattening. Ballone et al. (2020) further confirmed that simulated clusters initialized with significant rotation develop substructuring and show a high degree of fractality.

The dynamical evolution of rotating stellar systems can also be studied using 2D orbit-averaged Fokker-Planck (FP) models (Einsel & Spurzem 1999). FP models are extensions of the King model (King 1966) and based on the numerical solution of the orbit-averaged FP equation (Goodman 1983; Lupton & Gunn 1987; Longaretti & Lagoute 1996). These models explained the influence of rotation and its effects on the dynamical evolution of star clusters. Further works based on FP models include studies on binary heating and a stellar mass spectrum (Kim et al. 2002, 2004, 2008; Fiestas et al. 2006). Fiestas et al. (2006) stated that globular clusters in the Milky Way were initially flattened due to rotation and that their current flattening can also be explained by bulk rotation. This also suggests that open clusters may form with a certain amount of bulk rotation.

In addition, N-body simulations have also been used to investigate rotating star clusters. Ernst et al. (2007), Kim et al. (2008), and Hong et al. (2013) used direct N-body simulations to evaluate the accuracy of the numerical solutions of the FP model. They generally found a fair agreement between both models in all three phases of cluster evolution (core collapse, post-collapse, final dissolution). Rotation in nuclear star clusters has been studied with N-body models and FP models in Fiestas & Spurzem (2010) and Fiestas et al. (2012), although their main focus was the evolution of black holes. The N-body models of Tiongco et al. (2022) and Livernois et al. (2021, 2022) were all carried out with a larger particle number (N = 50 000) than earlier studies. These studies, however, initialized their models with a limited stellar initial mass spectrum (only stars with masses less than one solar mass). Kamlah et al. (2022b) performed N-body simulations of rotating star clusters with and without stellar evolution. In comparison to previous works, Kamlah et al. (2022b) used a more realistic initial mass function with stellar masses up to 150 M_⊙ and 1.1 × 10⁵ particles. They examined the effects of different initial bulk rotation speeds on systems with an initial binary fraction of 10%. In their analysis, they focused on the impact of the coupled phenomenon of the gravothermal-gravogyro catastrophe and its relation to rotating systems.

Despite the significant efforts to study the impact of bulk rotation on the evolution of star clusters and their constituents, there are still areas where understanding can be improved by further theoretical analysis. Hydrodynamical simulations are the most detailed and realistic, but they are computationally expensive and only relevant during the early stages of star cluster evolution. FP simulations and N-body simulations, on the other hand, produce similar results for the models with moderate bulk rotation, but have slight differences in models with high rotation (Hong et al. 2013). Several of the N-body studies have used an equal-mass spectrum or an initial mass function with stellar masses below 1 M_⊙ mass, which affects the overall dynamics (Tiongco et al. 2018, 2022; Livernois et al. 2021, 2022; Scalco et al. 2023). With the exception of Kamlah et al. (2022b), previous studies do not account for the presence of primordial binaries. The impact of binary systems is particularly important, considering the gravothermal-gravogyro catastrophe that leads to contraction and core collapse. Binaries can delay the contraction and affect the dynamical evolution of rotating systems. Evaluating the impact of binarity is therefore important for a better understanding of the evolution of rotating star clusters. Finally, earlier numerical studies tend not to compare the outcomes of the simulations with observations (e.g., absolute magnitudes and spectral energy distributions) which is the final objective of theoretical studies such as these.

In this paper, we investigate the dynamical evolution of star clusters with different primordial binary fractions, metallicities, and initial bulk rotations. We identify the rotation signature and discuss the effect on the dynamical evolution of the star clusters. This work highlights the impact of primordial binaries and metallicities on star cluster dynamics and the evolution of the 3D morphology of star clusters. Additionally, we present our data in the form of synthetic observations for future comparison with observations.

The paper is organized as follows. In Section 2, we describe the initial conditions and the N-body simulations. In Section 3, we present and discuss the results in terms of rotation signatures (Section 3.1), dynamical evolution (Section 3.2), morphology (Section 3.3), and synthetic observations (Section 3.4). In Section 4, we analyze the star cluster properties simultaneously and discuss their correlations. Finally, in Section 5, we summarize and conclude the findings of our work.

2 Methodology

We simulated 20 star cluster models using NBODY6++GPU, a high-performance GPU-accelerated code for gravitational N-body simulations (Kamlah et al. 2022a). NBODY6++GPU is a successor of the direct force integration N-body codes originally written by Sverre Aarseth (Aarseth 1985, 1999, 2003; Spurzem 1999; Aarseth et al. 2008). Using these direct N-body simulations, we studied the impact of initial rotation, stellar evolution, binary fraction, and the combination of catastrophes. Our models use level C stellar evolution in the form presented by Kamlah et al. (2022a).

The models were initialized using McLuster (Küpper et al. 2011; Kamlah et al. 2022a; Leveque et al. 2022). McLuster can set up initial conditions of N-body simulations or generate star clusters from observations (Küpper et al. 2011).

We used FOPAX (Einsel & Spurzem 1999; Kim et al. 2002, 2004, 2008) to generate the 2D FP initial models that add rotation to initialized N-body simulations, following a similar approach as in Einsel & Spurzem (1999) and Hong et al. (2013). The code produces a 2D mesh from the output of density ρ and velocity dispersions σ as a function of r and z of a rotating King model f (E, J_z). The density model is based on the King model (King 1966) that can be expressed as $$f_{rk}\propto (e^{\beta E}-1)\times e^{-\beta {\mathrm{\Omega }}_0{J}_z}.$$(1)

Here, β = (mσ²)⁻¹ is the anisotropy parameter (where σ is the velocity dispersion and m is the stellar mass), ω₀ = Ω₀ ${\omega }_0={\mathrm{\Omega }}_0\sqrt{9\pi G{n}_c/2}$ is the angular velocity, and W₀ = βm(ψ - ψ_t) is the King parameter, in which ψ and ψ_t are King potentials at the center and at the truncation radius, respectively (Lupton & Gunn 1987; Einsel & Spurzem 1999). King models are usually described by the pair of parameters (W₀, ω₀) (Einsel & Spurzem 1999). For our models, we adopted W₀ = 6.0 and ω₀ ∈ [0.0,0.6,1.2,1.8].

We used a Monte Carlo rejection technique to generate the discrete system of N particles that follow the distributions of r and z. All the output is in N-body format (one line per particle, containing the mass, 3D position, and 3D velocity data). We combined the generated N-body velocity distribution with the position distribution generated by MCLUSTER (see above). All the data were then scaled to standard Hénon units. Binary orbital parameter distributions were retained from MCLUSTER. Finally, we obtained the rotating King model with binary systems and the Kroupa (2001) initial mass function.

Initial conditions of the simulations are summarized in Table 1. The initial total number of single stars and binary systems in each system is n = N_single + N_bin = 10 000. We adopted a binary fraction f_bin, which is the ratio of the cluster members that are binary systems: $$f_{\mathrm{b}\mathrm{i}\mathrm{n}}=\frac{N_{\mathrm{b}\mathrm{i}\mathrm{n}}}{N_{\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{e}}+N_{\mathrm{b}\mathrm{i}\mathrm{n}}}\times 100\mathrm{\%}.$$(2)

Here, N_bin and N_single are the number of binary systems and single stars, respectively. We modeled three binary fractions f_bin = 0%, 10%, and 50%, such that the total number of stars in these models is N = N_single + 2N_bin = 10 000, 11 000, and 15 000, respectively. The corresponding initial total cluster masses are 5.8 × 10³ M_⊙, 6.4 × 10³ M_⊙, and 8.79 × 10³ M_⊙. Binary systems are generated using a uniform mass ratio distribution (0.1 ≤ q = m₂/m₁ ≤ 1.0) for m₁ > 5 M_⊙, and random pairing is adopted for the remaining binary systems. The semimajor axes are distributed uniformly in log a between the sum of the two stellar radii and 100 au. We adopted a thermal eccentricity distribution, f (e) = 2e. We adopted a Kroupa (2001) initial mass function with stellar masses between 0.08 M_⊙ and 150.0 M_⊙. For each of these models, we studied two different metallicities: metal-poor clusters (Z = 0.01) and clusters with solar metallicity (Z = 0.02). Finally, we modeled four degrees of bulk rotation ω₀ = 0.0, 0.6, 1.2, and 1.8.

We initialized the models with N = 10 000, 11 000, and 15 000 stars with half-mass radii of r_h = 2.25, 2.03, and 2.09 pc, respectively. The initial distributions of King models become more concentrated for larger values of ω₀ (see Fig. 1 in Einsel & Spurzem 1999). After initializing the rotating clusters with the input parameters of the King model with MCLUSTER, structural parameters from McUSTER such as r_h are slightly changed. Thus, the initial half-mass radii of the model realizations are r_h ≈ 2-2.5 depending on the choice of ω₀ (the second column from the left in Table 1). The reader could notice that r_h of models with f_bin = 0% are larger compared to those with binary fractions. This is because clusters with binaries tend to be smaller compared to those without. However, when the clusters rotate very fast, the binding effect of binaries wears off. Thus, r_h is larger in the modes with high ω₀. This affects clusters with f_bin = 0% too, but it is mostly noticeable in the model with ω₀ = 1.8. The stochastic nature of the initial mass function sampling and random selection of particle coordinates also contribute to an average 0.1 pc deviation in r_h. The tidal field is modeled as a point mass galaxy with a mass ≈1.4 × 10¹¹ M_⊙, and the clusters are placed at a circular orbit with a distance ≈8.178 kpc from the Galactic center. The models with 10k, 11k, and 15k stars have initial tidal radii of 21.65 pc, 22.37 pc, and 24.81 pc, respectively. We evolved all models until dissolution, which typically occurs after 1−2.5 Gyr.

Particles are removed from the simulation and labeled as escapers once they reach a distance equal to twice the tidal radius at that time. These escapers are excluded from the calculations of any of the star cluster global parameters, such as r_t or m_c.

Fig. 1 Angular momentum of all stars (panel a), binary stars (panel b), and single stars (panel c) of the simulations with 15k stars, fbin = 50%, and Z = 0.01. Colors represent simulations with different ω0. The line styles represent the angular momentum of stars within the inmost, intermediate, and outmost regions.

3 Results

3.1 Rotation signatures

3.1.1 Angular momentum

We quantify the rotation signature using the total angular momentum vector of the star cluster, ${\mathbf{L}}_{\text{total}}=\sum _i{M}_i({\mathbf{r}}_i\times {\mathbf{V}}_i)$. For single stars, M_i is the mass, r_i is the position relative to the stellar density center of the cluster, and V_i is the velocity, relative to the stellar density center of the cluster. For binary systems, M_i is the sum of the two stellar masses, and r_i and V_i correspond to position and velocity relative to the center of mass of each binary system. The center of mass is adjusted after every energy check in the simulation. We also calculate the angular momentum of the single stars and binary systems within different Lagrangian radii, φ, and examine the angular momentum of the inmost (φ = 0.0-0.3), intermediate (φ = 0.3-0.6), and outmost (φ = 0.6-1.0) regions of the system. We note that we used the total cluster mass at the time in the calculations of Lagrangian radii.

Fig. 1 shows the angular momentum of all stars (panel a), of the binary stars (panel b), and of the single stars (panel c), for the simulation with 15k stars, f_bin = 50%, and Z = 0.01 as an example. The angular momentum is computed considering the z-axis as the rotation axis. We note that only stars M_ZAMS < 6M_⊙ were used for the calculation of the angular momentum of the total system, the population of binaries, and the population of single stars. This is the optimal subset of stars for calculations since M_ZAMS ≥ 6 M_⊙ stars undergo supernova explosions due to the stellar evolution (further explanations are in Sect. 3.1.2).

Panel a of Fig. 1 shows the total angular momentum of the models as a function of time. All the parts of the rotating systems experience moderate changes up until 10 Myr. The model without bulk rotation has angular momentum close to zero in all its parts, and changes in the outer part to a negative direction after 7 Myr. The angular momentum of the inmost parts of all models is seen to be decreasing, and the rate of decrease is higher for models with more bulk rotation. The angular momentum of the intermediate parts of the systems is also seen to be decreasing, but the model with ω₀ = 0.6 shows some spike increase at ≈1.5 Myr. There is no significant change in the angular momentum for the outmost parts of all models, except for the model with ω₀ = 1.8. Such changes of the angular momentum for different parts of the cluster indicate that angular momentum is transported from the inner region to the outer parts of the cluster. Thus, we observed a decrease in the angular momentum of the inmost and intermediate parts. The outmost part receives the momentum from the inner and intermediate parts, but it can be clearly shown only in the model with ω₀ = 1.8. The reason for this is that models with high rotation tend to have larger core mass and core radius (see Fig. 6). Thus, the retention of stars in outmost parts of models with ω₀ = [0.6, 1.2] is lower. Smaller retention of stars does not allow clusters with ω₀ = [0.6, 1.2] to hold the angular momentum in the outmost part for long, and we did not see an increase as in the model with ω₀ = 1.8. The process of angular momentum transport results in a core collapse that can be seen in the decrease of core radius in Fig. 6. This is a gravogyro catastrophe described in the work of Hachisu (1979). After 10 Myr, the angular momentum of the inmost and intermediate parts reaches the constant values that they will hold on until ≈400 Myr. The angular momentum of outmost parts in the models is decreasing. This is further followed by the slight increase in momentum at ≈45 Myr in the models with ω₀ = [1.2, 1.8]. Further, the outmost part shows a constant decrease and reaches even negative values, and the inmost and central parts reach 0 pc² Myr⁻¹ by the time of dissolution. A similar decrease in angular momentum was seen in the work of Livernois et al. (2022). However, negative angular momentum is not seen in the study, since the runtime of their simulations is until the core collapse.

The panels in rows b and c of Fig. 1 show the angular momentum of the single stars and binary systems of the models. Trends seen in total angular momentum are also seen in the angular momentum of binaries and singles. Angular momentum of binaries is a larger contributor to the total angular momentum compared to single stars. This is mostly because binaries contain twice the mass of single stars in this model with f_bin = 50%.

We transformed the coordinates into the angular momentum vector’s direction at each time and inspect stars by their velocities such that the net angular momentum of the star cluster points in the direction of the z⁰-axis. Fig. 2 presents the transformed spatial coordinates for models with ω₀ = 1.8 at 0Myr (panels a), 10 Myr (panels b), and 50 Myr (panels c). The color map shows the velocity of the stars towards the transformed direction. We note that we limit displayed velocities from ≈10 km s⁻¹ to 10 km s⁻¹ to lessen the impact of escapers. From this, we can see that the model with ω₀ = 1.8 initially has a noticeable rotation signature at 0Myr and 10, Myr. v′_x and v′_y are seen in two distinctive regions with positive and negative velocities. As for the v′_z, it does not display similar positive and negative regions. This is because the z′-axis is still the rotation axis. These results are consistent with the hydrodynamical simulations of Mapelli (2017) and Ballone et al. (2020), who show similar results but for clusters from a simulated molecular cloud at 3 Myr and 5 Myr. However, we can also show the velocities at the later stages of evolution. At 50 Myr, it is seen that the regions are not as distinguishable, and the velocity differs in the range of ≈2 km s⁻¹ to 2 km s⁻¹. This is due to the reduction of the rotation signature over time, and it continues to decrease further. The two-body relaxation timescale in the models is ≈96-119 Myr, with clusters with a higher ω₀ having a longer relaxation time. The significant reduction in rotation signature by 50 Myr and further reduction at later times suggests that the effects of rotation are reduced when the clusters reach one relaxation time. The changes of the rotation signature from a few to 0.1 km/s are consistent with the current observation of open clusters by Hao et al. (2024), who measured the internal rotation of open clusters at the level of 0.1 km/s.

Fig. 2 Spatial coordinates for the model with ω0 = 1.8 at 0Myr (panel a), 10 Myr (panel b), and 50 Myr (panel c). The color map shows the component of the stellar velocity along the x′-, y′-, and z′-axes (where the z′ axis is the direction of the net angular momentum).

3.1.2 Escapers

In order to interpret the rotation signature, it is crucial to investigate the impact of escapers on the stellar system. As a star cluster evolves, angular momentum is transported from the inner regions to the outer regions of the cluster and is subsequently removed from the system through escaping stars and binary systems.

Aside from the escapers that stop being gravitationally bound to a system, the stellar evolution of the stars also affects the evolution of the system. In our analysis, we divide the stars into four groups by the zero-age main-sequence mass of stars (ZAMS): M_vlm. M_lm, M_mm, and M_hm, which correspond to the following mass ranges: $$\begin{array}{l}{M}_{\text{vlm}}:0.08M_{\odot }\le m_{\text{ZAMS}}<0.9M_{\odot }\text{}\\ M_{\text{lm}}:0.9M_{\odot }\le m_{\text{ZAMS}}<6M_{\odot }\text{}\\ M_{\text{mm}}:6M_{\odot }\le m_{\text{ZAMS}}<15M_{\odot }\text{}\\ M_{\text{hm}}:15M_{\odot }\le m_{\text{ZAMS}}<150M_{\odot }.\end{array}$$

Here, m_ZAMS is the ZAMS mass of a single star. The mass ranges are chosen to separate stars with different evolutionary tracks. Stars with masses M_vlm will remain main-sequence stars throughout the entire simulation. Most of those with M_lm will become white dwarfs (WDs), those with M_mm become neutron stars, and those with M_hm become black holes. However, we note that these definitions are not strict, and exceptions do exist.

Fig. 3 shows the sum of the angular momentum (top panel) and number of stars (bottom panel) of the four groups of stars as a function of time of the simulations with 15k stars, f_bin = 50%, and Z = 0.01. The M_vlm stars follow the global trend of the angular momentum of Fig. 1. This is the largest group of stars in the cluster, as seen in the bottom panel of Fig. 3. The 1500 stars in the M_lm group show some fluctuations in the angular momentum, but the trend of angular momentum is unaffected. For the stars in the M_mm and M_hm groups, there are strong fluctuations around 20-40 Myr. Considering the comparatively low number of stars (less than 100) in the groups M_mm and M_hm, the statistical fluctuations are considerably high. The number of stars in these two groups suffers from a sharp decrease during this period, which is due to stellar evolution. These stars undergo supernova explosions, that lead to observed extreme values in angular momentum (spikes in the upper panel). These fluctuations might bias our analysis of the angular momentum of the cluster. Therefore, we excluded stars with M_ZAMS ≥ 6 M_⊙ from the angular momentum calculations from Sect. 3.1.

Fig. 4 shows the properties of escaping stars such as the number of stars N_esc (panel a), the total mass of escapers M_esc (panel b), mean kick velocity v_k (panel c), and angular momentum L_esc (panel d) as a function of time. Unlike in previous calculations, the total angular momentum of escapers is calculated by $L_{\mathrm{e}\mathrm{s}\mathrm{c}}=\sum _i{M}_i{r}_i{v}_{t,i}$, where r_i and v_t,i is the distance to the center of mass and tangential velocity (v_t,i = v_k,i sin α, where v_k,i is the kick velocity, and α is the angle between the velocity vector and radius vector). The vectors are in the (x, y, z) coordinate system, of which the origin is the cluster center and the initial rotation of the cluster is oriented along the z direction. For the analysis, we bin the escapers into consecutive time intervals of of 2.5 Myr (until 10 Myr), 10 Myr (until 100 Myr), and 30 Myr (until dissolution; e.g., 0-2.5 Myr, 2.5-5 Myr). Each escaping star was assigned to its respective interval based on its escape time, and all necessary mean and total values were calculated accordingly. After 10 Myr, there are significant changes in the properties of the escapers. Firstly, the escape rate increases, which is accompanied by increases in Mesc, v̄_k, and L_esc) starting from 20 Myr. This occurs in all models, with the model with ω₀ = 1.8 having the largest number of escapers. The escapers at this time show high Mesc and v̄_k regardless of initial rotation. Extreme values are seen in their angular momentum that reaches over L_esc = 10³M_⊙pc²Myr⁻¹ or L_esc = 10⁻³M_⊙pc²Myr⁻¹ with 50-300 escaping stars. By comparing it to the bottom panel of Fig. 3 and considering the L_esc values, it is seen that those are influenced by the M_mm and M_hm stars. On further evolution, the parameter values of escapers show a stable decreasing rate for models with ω₀ = [1.2; 1.8] and an increasing rate for models with ω₀ = [0.0; 0.6] until dissolution.

Fig. 3 Evolution of the angular momentum (top panel) and number of stars (bottom panel) of four groups, for the simulations with 15k stars, fbin = 50%, and Z = 0.01. Colors scheme as in Fig. 1. The line styles represent groups of stars by MZAMS.

3.1.3 Change of rotation direction

We also examine the change in the angle, θ, between the cluster’s net angular momentum vector and the z-axis of the Galactic plane. Fig. 5 shows the angle, θ, as a function of time for the models with 15k stars, f_bin = 50%, and Z = 0.01.

The rotating models (ω0 > 0.0) start with an angle very close to 0°, since the rotation axis is initialized along the z-axis. The evolution continues at a similar rate in all the rotating models, but the model with ω₀ = 1.8 experiences moderate changes to θ ≈ 25° at 400 Myr. This is followed by the increase of angle θ until it reaches ≈175° by the age of 0.9-1 Gyr. The evolution in models with ω₀ = [0.6, 1.2] proceeds similarly, but the changes occur later compared to a model with ω₀ = 1.8. Such a change is the indication that the models keep rotating along the z-axis, but in the opposite direction. This is the effect of the gravitational tidal field and the low number of stars in the system. These results explain the negative angular momentum seen in Fig. 1.

The value of θ is not relevant for the nonrotating model. Even though the θ for the nonrotating model (ω₀ = 0.0) starts at θ = 38°, its total angular momentum is induced by the random motions of stars, which gives an arbitrary angle higher than 0°. The value of θ significantly increases after 10 Myr, ultimately reaching θ ≈ 175° by 100Myr.

General minor changes in θ in rotating models are somewhat consistent with the previous work of Tiongco et al. (2018), but the change to the opposite side of the rotation axis is not seen there. The possible reason is that their simulations end earlier (upon reaching the 75% mass loss) compared to the duration of our simulations (see Fig. 6). Also, their simulations have 32 768 particles with equal mass.

Fig. 4 Properties of escaping stars over time: Evolution of the number of stars Nesc (panel a), total mass Mesc (panel b), mean kick velocity v̄k (panel c), and angular momentum Lesc (panel d). The color scheme is the same as that of Fig. 1.

Fig. 5 Angle, θ, between the angular momentum vector and the z-axis as a function of the time of the models with 15k stars, fbin = 50%, and Z = 0.01. The color scheme is the same as in Fig. 1.

3.2 Dynamical evolution of binary systems and single stars

In this section, we investigate the impact of the binary fraction, metallicity, and initial bulk rotation on the dynamical evolution of star clusters. We examine the change in global parameters and the rate of change of the number of objects present in the system.

Figure 6 shows the evolution of cluster parameters such as the tidal radius r_t (panel (a)), the core radius r_c (panel (b)), the mass of the core m_c (panel (c)), and the half-mass radius r_h (panel (d)) in four panels from top to bottom. The core radius is the fundamental global parameter used for the characterization of the central structure of star clusters. The core mass is the mass of the stars within the core radius. The example models are the 4 models with 15k stars, f_bin = 50%, Z = 0.01 and different ω₀. Colors represent the bulk rotation of the models.

First, there is the time evolution of r_t in panel a of Fig. 6. r_t decreases over time, and it is seen that the models with ω₀ = 1.2, 1.8 show a faster decrease over time compared to models with ω₀ = [0.0, 0.6]. However, during further evolution, models with ω₀ = [0.0, 0.6] reach and even surpass the rate of change of r_t of models with high rotation. This leads to earlier cluster dissolution.

Panels b-d of Fig. 6 show the time evolution of the evolution of r_c, m_c, and r_h. Initially, r_c and m_c are higher among the models with greater rotation. Further, at the age of ≈2-5 Myr, there is a huge decrease in r_c and m_c among all the models. This is the core collapse. Noticeably, the models with higher rotation experience much stronger core collapse. At the age of 5-10 Myr, there is another increase of r_c and m_c, especially seen in the models with high rotation. The second core collapse further follows this at ≈10-15 Myr in the models with high rotation. The early evolution of r_h up to 400 Myr is similar in all the models. The higher r_h is seen in models with high rotation. One particular increase is seen at ≈7-10 Myr in the model with ω₀ = 1.8. This time corresponds to the increase in r_c and m_c . Further in the evolution, the clusters with higher bulk rotation have decreased r_h . This is due to the high mass loss of these models (see Fig. 7).

It is known from previous studies that binary stars act as gravitational energy sources (Aarseth 1985; Heggie 1975, 1984; Hénon 1975; Trenti et al. 2007). The dynamical formation of binaries primarily occurs near the center of the cluster (Bettwieser & Sugimoto 1984). Additionally, binary stars should supposedly halt the contraction, affecting the gravogyro catastrophe (Aarseth 1985; Heggie 1975, 1984). The high binary fraction can also delay or prevent the core collapse (Kamlah et al. 2022b). Additionally, clusters with high primordial binary fractions have high gravitational binding energy, which means a higher chance of retaining stars and a longer cluster lifetime. For all these reasons, we observe the impact of primordial binaries on rotating clusters.

We investigate the rate of change of the number of objects. Fig. 7 shows the number of stars and binary systems (n = N_single + N_bin; panel a), the fraction of remaining binary systems (N_bin/n₀; panel b), single stars (N_single/n₀; panel c), and hot stars with T_eff > 20 000 K (N_hot/n₀; panel d) as a function of time for the models with ω₀ = 0.6 and ω₀ = 1.8. n₀ is the sum of the initial number of single stars and the initial number of binary systems. Finally, we define stars with T_eff > 20 000 K as hot stars.

From the rate of change in the total number of particles (panel a), binary systems (panel b), and single stars (panel c) of Fig. 7, It is seen that the models with high f_bin and solar metallic-ity survive the longest. The models with similar Z show almost the same rate of change of binaries and single stars. Nonetheless, models with ω₀ = 1.8 show a decrease in the number of particles starting at 40 Myr, while models with ω₀ = 0.6 start to decrease after 100 Myr. A decrease is seen in both single stars and binary systems. In the models with faster rotation ω₀ = 1.8 (solid curves), the fraction of both binaries and singles decreases quicker than the models with slower rotation ω₀ = 0.6 (dashed curves).

The hot stars (panel d of Fig. 7), are stars with T_eff > 20 000 K. Models with an identical particle number but different abundance Z have a similar number of massive stars. The model with 15k stars, f_bin = 50%, and Z = 0.01 has the highest number of hot stars, and the model with 11k stars, f_bin = 10%, and Z = 0.02 has the lowest amount of hot stars. Models with higher Z tend to have fewer hot stars compared to the clusters with low Z . The number of hot stars is higher for models with a higher initial binary fraction. For instance, the initial ratio of hot binary stars is more than 65% for models with N = 15k, f_bin = 50%, and both Z options. It is seen that the number of hot stars reaches its lowest value at around 30-40 Myr for all models. This is followed by an increase of up to 200 Myr and a constant decrease until the end of the dynamic evolution. The hot stars are mostly the most massive stars, and depletion by 30-40 Myr can be explained by the short lifetime of massive stars and binary evolution. After this, they evolve into compact objects, such as neutron stars. The second peak of hot stars occurs due to the formation of carbon/oxygen WDs (COWD) and oxygen/neon (ONWD) WDs.

Fig. 6 Evolution of several star cluster parameters, tidal radius rt (panel (a)), the core radius rc (panel (b)), the core mass mc (panel (c)) and the half-mass radius rh (panel (d)), for simulations with 15k stars, fbin = 50%, Z = 0.01. Different colors represent simulations with different ω0.

Fig. 7 Number of single stars and binary systems (n = Nsingle + Nbin; panel a), the fraction of remaining binary systems (Nbin/n0; panel b), single stars (Nsingle/n0; panel c), and hot stars with Teff > 20 000 K (Nhot/n0; panel d) as a function of time for the models with ω0 = 0.6 and ω0 = 1.8. The term n0 is the initial number of single stars and binary systems. Nbin, Nsingle, Nhot are the number of remaining binary systems, single stars, and hot stars. The colors represent models with different N (total initial number of stars), fbin (initial binary fraction) and Z (metallicity). Line styles represent simulations with different ω0.

3.3 3D morphology evolution

In this section, we investigate the shape of the star clusters in the simulations and the observations and compare our results with the observational data of 90 clusters obtained from Gaia DR3 (Gaia Collaboration 2023) with the membership from Pang et al. (2021, 2022b, 2024). We investigate the morphology by fitting an ellipsoid to the cluster’s shape and use its principal axes: the major axis a, the intermediate axis b, and the minor axis c, to quantify the cluster’s structural and morphological properties at the different times. The principal axis ratios of the intermediate-to-major axial ratio b/a, the minor-to-major axial ratio c/a, and the minor-to-intermediate axial ratio c/b characterize the degree of elongation and ellipticity of the cluster towards the y and z directions. We do note that we align a, b, c axes to x, y, z axes. This would show the development of morphological structures in the x-y plane while z is the rotational axis.

The simulated clusters of Kamlah et al. (2022b) develop a bar-like structure during the gravogyro catastrophe. To compare our results to those of Kamlah et al. (2022b), we fit the simulated cluster with a 3D ellipsoid model. The code of the ellipsoid fitting is adapted from the Ellipsoid-Fit code published on Github by R. M. Sandu¹. The fitted axis values depend on the number of stars that are included for the fitting. Stars in the inner part of the cluster do not contribute much to the ellipsoid fitting. Therefore, we only use the outer 10% of stars within one tidal radius for the ellipsoid fitting. The ellipsoid fitting does measure the shape in the outer region of the cluster, which is influenced by the tidal field. However, the effects of internal rotation induce global changes in morphology, which also affect the outer part of the cluster, as can be seen from the decline of the axis ratio. For illustration, we display the ellipsoid fitting of the Pleiades cluster, whose members are taken from Pang et al. (2022b), as an example in Fig. 8. Blue and black scatter points are the member stars and the stars used for the ellipsoid fitting (outer 10% of stars within the tidal radius), respectively.

Fig. 9 presents the evolution of the principal axis ratios b/a (upper row), c/a ratios (middle row), and c/b ratios (bottom row) of the simulated clusters and observations. The axis ratios of simulated clusters are the results of the bootstrapping procedure. We generate 100 bootstrap datasets² from 1 snapshot and run ellipsoid fitting on each of the created datasets. We use mean values and standard deviation to construct the uncertainties in the principal axis ratios for the selected snapshot. We choose 3 snapshots between 0 to 10 N-body time units (0, 3, and 7 N-body time), 9 between 10 to 100 N-body time units (at a 10 N-body time interval), and 9 between 100 to 1000 N-body time units (at a 100 N-body time interval). Models have different scaling factors of N-body units depending on N . For the models with N = 15 000, one N-body time unit corresponds to a physical timescale of ≈0.63 Myr. For N = 10000 and N = 11000, the scaling factors are ≈0.76 and ≈0.73, respectively.

Data points with error bars in Fig. 9 are the observed clusters from Pang et al. (2021, 2022b, 2024). Error bars are obtained by subtracting the difference between calculated principal axis ratios of clusters with and without corrected distances of Pang et al. (2021, 2022b, 2024). The color pattern of observations is governed by the sample size used for the ellipsoid fitting. We note that the observation scatter points are identical in all panels.

Fig. 9 shows that models have different axis ratios depending on the initial conditions (especially ω₀) and on time. During the early evolution of the clusters, the largest differences in bootstrap ranges for models with different ω₀ are present in the models with N = 10k, f_bin = 0% and Z = 0.01. These models have similar axis ratios from the start, ranging between ≈85-94% for b/a, ≈73-84% for c/a, and ≈78-94% for c/b. As by definition (c ≤ b ≤ a), all models have c/a ≤ b/a. The model without rotation (ω₀ = 0.0) starts with a more or less spherical morphology, with c/b - c/a - b/a - 0.9 ± 0.042. The initial morphology in the models with stellar evolution from the work of Kamlah et al. (2022b) did show similar triaxiality (b ≠ a ≠ c) in the beginning.

Upon further evolution up until 10 Myr, all the models experience a decrease in axis ratios, which is the effect of the core collapse. The models without bulk rotation also have a moderate decrease in the axis ratios, as they develop tidal tails under the influence of the Galactic tidal field. For the models with bulk rotation, the decrease is earlier and sharper, especially for the models with high bulk rotation. This is followed by differing increases and further decreases of axis ratios in all the models up until 100 Myr. The models without bulk rotation and ω₀ = 0.6 have a increase in b/a, c/a and c/b by the age of 20 Myr. The models with ω₀ = 1.2 and ω₀ = 1.8 show the same increase and decrease but with earlier peak and sharp decrease by the 2025 Myr. This is especially apparent for the model with ω₀ = 1.8 that reaches its all-time low in c/a. This kind of dynamic persists to 200 Myr. Clusters achieve high axis ratios until they are close to dissolution. The major change in c/a indicates the existence of elongation and bar-like structure at the given moment. This kind of morphology was also apparent in the work of Kamlah et al. (2022b). Another similarity from Kamlah et al. (2022b) is the change from initial triaxiality to axisymmetry (b = a, but c ≠ a and c ≠ b) in further evolution. The most drastic elongation reaching over c/a - 0.4 can be seen for the model with ω₀ = 1.8. This model has the highest angular momentum of the system and shows high ellipticity up until 50 Myr. This is consistent with the findings of Lahén et al. (2020a), who investigated the structure and rotation of massive star clusters at 20 Myr. They found that for clusters with masses above 30 000 M, the higher the specific angular momentum, the higher the ellipticity (lower c/a).

All the other models with more particles and binaries perform similarly but have fluctuations at different timescales and with differing degrees. The second decrease in c/a, supposedly the second core collapse, is sharper and more apparent for the models without binaries. However, a similar decrease with a lesser scale is seen in the models with Z = 0.01 as shown in the 2nd column and 4th columns from left in Fig. 9. The decrease is still present in models with Z = 0.02, but it is even less apparent compared to models with Z = 0.01. The decrease is sharpest for the models with ω₀ = 1.8 with the only exception for the models with N = 11k and f_bin = 10%. The model with N = 11k, f_bin = 10%, Z = 0.01, and ω₀ = 1.2 shows the sharpest decrease in both b/a and c/a at 4Myr among all the models.

In terms of the observations, most of the selected observable clusters are inside the bootstrap ranges of the simulations. The ones that have higher axis ratios are the clusters with a low number of input stars for ellipsoid fitting. The exceptions are the NGC 1980, NGC 2516, and the Pleiades, which have more than 100 input samples and show both axis ratios above 0.9. The discrepancy of the observed clusters is ≈40%. Nonconformant clusters majorly tend to have input samples of less than 100 stars. The only exception is the NGC 2516, which has b/a, c/a and c/b of < 0.97, despite having input stars of ≈200. Those that have lower c/a compared to the simulation bootstrap ranges are mostly filamentary clusters such as Collinder 132 gp3, Collinder 132 gp5, and NGC 6405. They are young and elongated, thus showing c/a <0.6. They already formed the tidal tails and were considered dynamically older compared to the simulations at a similar age.

Fig. 8 Ellipsoid drawn for the Pleiades cluster. The ellipsoid is represented by the green wireframe. Blue and black scatter points are the member stars and input stars of Pleiades used for the ellipsoid fitting. Orange, pink, and red lines indicate the principal axes a, b, and c. The stellar members of the Pleiades are obtained from Pang et al. (2022b).

3.4 Synthetic observations for simulated clusters

Here, we present our simulations in three observational magnitude filters and spectral energy distribution (SED). We used GalevNB (Pang et al. 2016) to convert the physical properties of N-body simulation data into filter bands of Gaia, China Space Station Telescope (CSST), and Hubble Space Telescope/ Space Telescope Imaging Spectrograph (HST/STIS). This allows for direct comparison between numerical studies and observations. The SED is needed to investigate excess ultraviolet (UV) emission. We converted the data into spectra spanning far-UV to near-IR wavelengths. The SED was created by summing the fluxes of the stellar population. We further separated the summation of fluxes into two categories: single stars and binary systems.

Fig. 10 shows the CMDs of the simulated clusters in three observational filters for ages between 0 Myr to 1.3 Gyr for the model with 15k stars, f_bin = 50%, Z = 0.01, and ω₀ = 1.2. As can be seen from the figure, the magnitude ranges of CSST and HST are very similar to that of Gaia. We note that for a binary to be classified as a WD binary, at least one of the two components should be a WD. Initially, the majority of the single stars and binaries are on the main sequence. The number of binary systems is equal to the number of single stars. After 100 Myr, there are ONWDs that are mostly single stars. However, in further evolution, there are mostly COWDs. Single-star WDs formed the WD sequence on the bottom left side of the CMD. Binaries with WDs are either near the main sequence or between the main sequence and the WD sequence. The ones between them are the cataclysmic variable star binaries, which are interacting WD binaries that have an MS star companion with a mass comparable to the WD (Warner 1995; Pang et al. 2022a).

Fig. 11 shows SEDs for ages from 0Myr to 1.3 Gyr for the model with 15k stars, f_bin = 50%, Z = 0.01, and ω₀ = 1.2. The curves represent the total flux (blue), single star flux (green), and binary system flux (orange). In the models with f_bin = 50%, binaries have twice the flux the of the single stars.

Fig. 9 Evolution of principal axis ratios of the simulated clusters with age. The areas shaded in light colors enclose the minimum and maximum values from bootstrapping as a function of time, for intermediate-to-major axis ratio b/a (upper row), the minor-to-major axis ratio c/a (middle row), and the minor-to-intermediate axis ratio c/b (bottom row). Columns show the results of the different models (indicated at the top). The colors of the shaded areas indicate the bootstrap range for models with different ω0. Symbols with error bars indicate the principal axis ratios for observations calculated from the members of the star clusters studied in Pang et al. (2021, 2022b, 2024).

4 Discussion

In this section, we discuss and interpret our findings by examining several parameters simultaneously. Furthermore, we compare our results with findings from the literature.

Fig. 12 shows a comparison of the Lagrangian radii r_lagr (panels (a)), core radii r_core (panels (b)), total mass M_total and core mass M_core (panels (c)), principal axes ratios b/a and c/a (panels (d)), and angular momentum on the z-axis (panels (e)) of the model with 15k stars, f_bin = 50%, Z = 0.01, for different values of ω₀.

The Lagrangian radii (Fig. 12; top row) show that 90% of the cluster mass remains within a similar distance with a slight expansion over the course of the simulation. The radii that change a lot are those with φ = 0.9 and φ = 1.0. Notably, the φ = 1.0 radius increases drastically, and then remains more or less constant until dissolution. The expansion of this part starts as early as 3 Myr, but it is seen that it occurs sooner for the models with high ω₀. For example, for the model with ω₀ = 1.8 it starts expanding at 1 Myr and reaches its peak of ≈40 pc around ≈2.5–3 Myr.

Fig. 12 shows the evolution of r_core (row b), and M_total and M_core (row c). These panels allowed us to identify core collapses, as well as the changes in total cluster mass. A significant downward change in r_core and M_core occurs in all models. The change is greater in the models with ω₀ = [1.2, 1.8]. Simultaneously, M_total decreases significantly. It is seen that rapidly rotating models experience additional core collapses after 50 Myr. Another notable feature is that models with higher rotation have a higher core mass and core radius at the start of the simulations. This is consistent with the hydrodynamic simulations, where it was seen that the rotation signature is more apparent for the high-mass clusters (Mapelli 2017).

The effects of core collapse are also seen in the principal axes ratio evolution shown in row d of Fig. 12. At the same timescales of the core collapse shown in row b, there are major changes in c/a, showing the elongation among the models. The higher the bulk rotation, the more apparent the elongation at the exact time periods of angular momentum transport. Additionally, rapidly rotating models have another significant elongation at the timescale of the second core collapse. Elongation at this time is more prominent in the models without primordial binaries (see Section 3.3).

The effects of core collapse are seen in the angular momentum shown in the bottom row e of Fig. 12. Colors in row 5 represent angular momenta at different Lagrangian radii, where blue, green, yellow, and red correspond to angular momentum in the z-axis for all of the stars, the inmost part (φ = 0.0-0.3), central part (φ = 0.3-0.6), outmost part (φ = 0.6-1.0), and the entire system (φ = 0.0-1.0). Prior to the first core collapse, there is angular momentum transport from the inner part and the central parts to the outer regions of the clusters. This is more apparent in the model with ω₀ = 1.8. Further, the momentum of the outer part decreases until 50 Myr, where there is another increase that leads to the second core collapse. The consequences of this process can also be observed in other parameters, such as r_core and c/a.

All this shows the effects of the gravothermal-gravogyro catastrophe on the dynamic evolution of the star clusters. The gravogyro catastrophe can be noticed by the angular momentum transport. The culmination of the gravothermal-gravogyro catastrophe is the 1st core collapse that happens in every observed model. If the cluster happens to have high bulk rotation, then there might be a second core collapse that can be seen through the instability in angular momentum and 3D morphology. If the observed clusters also happen to have primordial binaries, it can somewhat hold out the second core collapse.

Fig. 10 Color-magnitude diagrams for three filter bands at ages between 0 Myr to 1.3 Gyr, for the cluster model with 15k stars, fbin = 50%, Z = 0.01, and ω0 = 1.2. Rows represent CMDs for different filters. Columns represent different ages. Blue and red circles represent the magnitude of single stars and binary systems. Green and orange crosses inside the circles are the COWDs and ONWDs.

Fig. 11 Spectral energy distributions for ages between 0Myr to 1.3 Gyr, for the model with 15k stars, fbin = 50%, Z = 0.01, and ω0 = 1.2. The colors of the curves represent the total flux, single star flux, and binary system flux as dashed blue, green, and orange lines respectively.

Fig. 12 Lagrangian radii rlagr (panels (a)), core radii rcore (panels (b)), total mass Mtotal and core mass Mcore (panels (c)), principal axes ratios b/a and c/a (panels (d)), and angular momentum on the z-axis (panels (e)) of the cluster model with 15k stars, fbin = 50%, Z = 0.01, and different ω0. Rows represent respective parameters, and columns represent the initial degree of rotation. The upper row shows Lagrangian radii, where color map colors correspond to the radius of different mass percentages. The colors in the 4th row show the principal axis for the stars within one tidal radius. Colors in the bottom row show the angular momentum in z-axis for all of the stars within the inmost part (φ = 0.0-0.3), central part (φ = 0.3-0.6), outmost part (φ = 0.6-1.0), and the entire system (φ = 0.0-0.1).

5 Conclusion

We examined the dynamical evolution of the rotating star clusters with differing primordial binary fractions, as well as different metallicities. We examined the evolution of the rotation signature, the impact of binaries and single stars, and the 3D morphology of the rotating star clusters. We also converted our simulations to the current observational filters for future comparison with observations. Our main findings can be summarized as follows:

• The models with a high bulk rotation experience two core collapses before the age of 20 Myr. This is due to the gravogyro-gravothermal catastrophes seen by angular momentum transport from the inner parts to the outer. There is an increase in the number of escapers at this period, and they are induced from M_mm and M_hm stars under stellar evolution;

• The effect of bulk rotation is noticeable during the early phases of evolution, but disappears by the timescale of two-body relaxation;

• There is a change in the direction of the angular momentum vector, which occurs early in the evolution for models with no rotation, and near the time of dissolution for models with rotation. This change is caused by the tidal field. The effects are stronger in nonrotating clusters, and rotating clusters are affected after the disappearance of the rotation signature and eventual mass loss;

• Bulk rotation affects the long-term dynamic evolution and strongly affects the morphology at the early stages of evolution. The effect is seen in high elongation at the timescales of angular momentum transport;

• We presented synthetic observations of the simulated clusters for future comparison with observations. These results can be used for the interpretation of future observations of Gaia, CSST, and HST.

Our study has several limitations that will be addressed in future studies. We have analyzed simulations with ~10₄ stars. We plan to extend our study with simulations of 50 000 particles to assess the effects of rotation and primordial binaries on more massive stellar systems. We will examine their dynamic evolution and attempt to confirm the relationship between the size of the cluster and the ellipticity and angular momentum transport described in Lahén et al. (2020a). We will similarly examine the dynamical evolution and convert the simulation results to different observational filter bands.

Data availability

The snapshot data of the models used in this study are available at Zenodo, https://zenodo.org/records/15637954. Uploaded snapshot files have basic data of single stars, binaries, and mergers at different ages, with an interval of 10 Myr before the age of 100 Myr. After 100 Myr until the end of simulations, the interval of snapshots is increased to 100 Myr. Refer to README file on Zenodo for instructions.

If more snapshot data or other output data is required, it can be requested from the corresponding author.

Acknowledgements

We thank the referee for the constructive comments that helped us to improve this paper in many ways. We thank Dr. Michela Mapelli and Dr. Jiadong Li for helpful discussions. Xiaoying Pang acknowledges the financial support of the National Natural Science Foundation of China through grants 12173029 and 12233013. We acknowledge the science research grants from the China Manned Space Project with No. CMS-CSST-2021-A08. This research is supported by the Science Committee of the Ministry of Science and Higher Education of the Republic of Kazakhstan (Grant No. AP19677351 and AP13067834). R.S. acknowledges NAOC International Cooperation Office for its support in 2023, 2024, and 2025; and the National Science Foundation of China (NSFC) under grant No. 12473017. This research was supported in part by grant NSF PHY-2309135 to the Kavli Institute for Theoretical Physics (KITP). R.S. acknowledges support by the German Science Foundation (DFG, project Sp 345/24-1). M.G. was supported by the Polish National Science Center (NCN) through the grant 2021/41/B/ST9/01191. A.A. 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. P.B. thanks the support from the special program of the Polish Academy of Sciences and the U.S. National Academy of Sciences under the Long-term program to support Ukrainian research teams grant No. PAN.BFB.S.BWZ.329.022.2023. P.B. also thanks the support under the project No. BR24992759 “Development of the concept for the first Kazakhstani orbital cislunar telescope - Phase I”, financed by the Ministry of Science and Higher Education of the Republic of Kazakhstan. Software: Astropy (Astropy Collaboration 2013, 2018, 2022), SciPy (Millman & Aivazis 2011), NBODY6++GPU (Kamlah et al. 2022a), McLuster (Küpper et al. 2011; Kamlah et al. 2022a), and FOPAX (Einsel & Spurzem 1999; Kim et al. 2002, 2004, 2008).

References

All Tables

All Figures

	Fig. 1 Angular momentum of all stars (panel a), binary stars (panel b), and single stars (panel c) of the simulations with 15k stars, fbin = 50%, and Z = 0.01. Colors represent simulations with different ω0. The line styles represent the angular momentum of stars within the inmost, intermediate, and outmost regions.
In the text

	Fig. 2 Spatial coordinates for the model with ω0 = 1.8 at 0Myr (panel a), 10 Myr (panel b), and 50 Myr (panel c). The color map shows the component of the stellar velocity along the x′-, y′-, and z′-axes (where the z′ axis is the direction of the net angular momentum).
In the text

	Fig. 3 Evolution of the angular momentum (top panel) and number of stars (bottom panel) of four groups, for the simulations with 15k stars, fbin = 50%, and Z = 0.01. Colors scheme as in Fig. 1. The line styles represent groups of stars by MZAMS.
In the text

	Fig. 4 Properties of escaping stars over time: Evolution of the number of stars Nesc (panel a), total mass Mesc (panel b), mean kick velocity v̄k (panel c), and angular momentum Lesc (panel d). The color scheme is the same as that of Fig. 1.
In the text

	Fig. 5 Angle, θ, between the angular momentum vector and the z-axis as a function of the time of the models with 15k stars, fbin = 50%, and Z = 0.01. The color scheme is the same as in Fig. 1.
In the text

	Fig. 6 Evolution of several star cluster parameters, tidal radius rt (panel (a)), the core radius rc (panel (b)), the core mass mc (panel (c)) and the half-mass radius rh (panel (d)), for simulations with 15k stars, fbin = 50%, Z = 0.01. Different colors represent simulations with different ω0.
In the text

Fig. 7 Number of single stars and binary systems (n = Nsingle + Nbin; panel a), the fraction of remaining binary systems (Nbin/n0; panel b), single stars (Nsingle/n0; panel c), and hot stars with Teff > 20 000 K (Nhot/n0; panel d) as a function of time for the models with ω0 = 0.6 and ω0 = 1.8. The term n0 is the initial number of single stars and binary systems. Nbin, Nsingle, Nhot are the number of remaining binary systems, single stars, and hot stars. The colors represent models with different N (total initial number of stars), fbin (initial binary fraction) and Z (metallicity). Line styles represent simulations with different ω0.

In the text

	Fig. 8 Ellipsoid drawn for the Pleiades cluster. The ellipsoid is represented by the green wireframe. Blue and black scatter points are the member stars and input stars of Pleiades used for the ellipsoid fitting. Orange, pink, and red lines indicate the principal axes a, b, and c. The stellar members of the Pleiades are obtained from Pang et al. (2022b).
In the text

Fig. 9 Evolution of principal axis ratios of the simulated clusters with age. The areas shaded in light colors enclose the minimum and maximum values from bootstrapping as a function of time, for intermediate-to-major axis ratio b/a (upper row), the minor-to-major axis ratio c/a (middle row), and the minor-to-intermediate axis ratio c/b (bottom row). Columns show the results of the different models (indicated at the top). The colors of the shaded areas indicate the bootstrap range for models with different ω0. Symbols with error bars indicate the principal axis ratios for observations calculated from the members of the star clusters studied in Pang et al. (2021, 2022b, 2024).

In the text

	Fig. 10 Color-magnitude diagrams for three filter bands at ages between 0 Myr to 1.3 Gyr, for the cluster model with 15k stars, fbin = 50%, Z = 0.01, and ω0 = 1.2. Rows represent CMDs for different filters. Columns represent different ages. Blue and red circles represent the magnitude of single stars and binary systems. Green and orange crosses inside the circles are the COWDs and ONWDs.
In the text

	Fig. 11 Spectral energy distributions for ages between 0Myr to 1.3 Gyr, for the model with 15k stars, fbin = 50%, Z = 0.01, and ω0 = 1.2. The colors of the curves represent the total flux, single star flux, and binary system flux as dashed blue, green, and orange lines respectively.
In the text

Fig. 12 Lagrangian radii rlagr (panels (a)), core radii rcore (panels (b)), total mass Mtotal and core mass Mcore (panels (c)), principal axes ratios b/a and c/a (panels (d)), and angular momentum on the z-axis (panels (e)) of the cluster model with 15k stars, fbin = 50%, Z = 0.01, and different ω0. Rows represent respective parameters, and columns represent the initial degree of rotation. The upper row shows Lagrangian radii, where color map colors correspond to the radius of different mass percentages. The colors in the 4th row show the principal axis for the stars within one tidal radius. Colors in the bottom row show the angular momentum in z-axis for all of the stars within the inmost part (φ = 0.0-0.3), central part (φ = 0.3-0.6), outmost part (φ = 0.6-1.0), and the entire system (φ = 0.0-0.1).

In the text