© The Authors 2026

1 Introduction

Numerous observational studies in the last few decades have provided new insights into the stellar content of globular clusters and revealed that these systems host multiple stellar populations differing in the abundance of several light elements (such as Na, O, Al, Mg, C, N, and He; see e.g. Carretta et al. 2009a,b; Marino et al. 2019; Carretta 2019; see also Gratton et al. 2012 and references therein). About 20% of Galactic clusters have been found to be even more complex and also display a spread in Fe (see e.g. Milone et al. 2017). Despite the extensive investigations and progress made in this field, many questions concerning the origin of multiple populations and the formation history of globular clusters are still unanswered.

The various formation pathways proposed in the literature consider different sources of processed gas for the formation of the chemically anomalous population including asymptotic giant branch stars, massive binary stars, rapidly rotating stars, and supermassive stars (see e.g. de Mink et al. 2009; D’Ercole et al. 2010; Krause et al. 2013; Denissenkov & Hartwick 2013; Wang et al. 2020; Nguyen & Sills 2024). Models also differ in the star formation history of multiple populations with some scenarios proposing separate star formation events (see e.g. D’Ercole et al. 2008; de Mink et al. 2009; Denissenkov & Hartwick 2013; D’Antona et al. 2016) and others in which the different populations are instead formed at the same time (see e.g. Bastian et al. 2013; Gieles et al. 2018). Hereafter, we refer to chemically anomalous stars formed from processed gas as second-generation (SG) stars and to stars with a chemical composition similar to that of halo field stars as first-generation (FG) stars.

Despite these differences, most models predict that SG stars would form centrally concentrated in the inner regions of a more spatially extended FG system (see e.g. D’Ercole et al. 2008; Calura et al. 2019). A few studies have investigated the SG formation in a FG cluster characterised by internal rotation and have shown that in this case the two populations also differ in their kinematic properties, with the SG system being characterised by a more rapid rotation than the FG one (see e.g. Bekki 2010, 2011 ; Lacchin et al. 2022).

In order to link the present-day dynamical properties of FG and SG stars with those imprinted at the time of their formation, it is necessary to study how various early and long-term evolutionary processes affect the structure and kinematics of globular clusters and their multiple populations. A number of studies have been focused on the theoretical investigation of the evolution of the initial dynamical differences of FG and SG subsystems, their dependence on the initial conditions, and the role played by various dynamical processes (see e.g. Vesperini et al. 2013; Mastrobuono-Battisti & Perets 2013, 2016; Hénault-Brunet et al. 2015; Vesperini et al. 2018; Tiongco et al. 2019; Livernois et al. 2024; Hypki et al. 2025; Berczik et al. 2025; Aros et al. 2025). On the observational side, many studies have detected spatial and kinematic differences in the present-day properties of several Galactic globular clusters (see e.g. Bellini et al. 2015; Milone et al. 2018; Dalessandro et al. 2019, 2021, 2024; Libralato et al. 2023; Leitinger et al. 2023; Martens et al. 2023; Cordoni et al. 2025; Ziliotto et al. 2025; Griggio et al. 2025) and revealed trends between the extent of these differences and the clusters’ dynamical ages (see e.g. Dalessandro et al. 2018, 2019, 2024; Libralato et al. 2023; Cordoni et al. 2025).

The goal of this paper is to investigate the long-term dynamical evolution of multiple generations in globular clusters with internal rotation through the use of N-body simulations. Specifically, we seek to explore the spatial and kinematic properties of these multiple generations, how they depend on stellar mass, the role of the initial orientation of the internal rotation axis relative to the cluster’s orbital angular momentum, and what this can reveal about the cluster’s initial conditions.

The outline of this paper is as follows. In Section 2, we describe the methods and initial conditions of our simulations. In Section 3, we describe the results of our analysis. Finally, in Section 4, we summarise the conclusions drawn from our analysis.

2 Methods and initial conditions

We performed a set of N-body simulations using NBODY6++GPU (Wang et al. 2015) on the Big Red II and Big Red 200 supercomputers at Indiana University. The models in our simulations move on circular orbits around the centre of their host galaxies, and the potential of the host galaxy is modelled as that of a point mass.

Each of our models begins with 100 000 stars with masses between 0.1 and 1 M_⊙ distributed using a Kroupa (2001) initial mass function. Our investigation is focused only on the longterm evolution driven by two-body relaxation and the effects of the host galaxy tidal field; effects due to mass loss associated with stellar evolution are not included, and, given the mass range included in our models, they would be expected to be small. Our N-body models begin with an equal number of FG and SG stars; in a number of scenarios for the formation of multiple stellar populations, the initial number of SG stars is smaller than the number of FGs, but the SG-to-FG ratio rapidly increases to values close to the one adopted here or larger during the cluster’s early evolution (see e.g. D’Ercole et al. 2008; Vesperini et al. 2021; Sollima 2021; Hypki et al. 2025).

The initial structural properties of our models are characterised by a rapidly rotating, flattened, and centrally concentrated SG subsystem embedded in a more diffuse, slowly rotating FG system. These structural properties are based on the results of a number of studies of hydro and N-body simulations following the formation of multiple stellar populations in rotating clusters (Bekki 2010, 2011; Lacchin et al. 2022).

Our initial conditions have been set using the code MaGalie (Boily et al. 2001): the SG subsystem is modelled as a rotating exponential disc, while the FG system follows a spherical distribution with a Hernquist density profile (Hernquist 1993). Rotation is included in the FG system by following the Lynden-Bell demon procedure (Lynden-Bell 1960): we reverse the rotational velocity around the rotation axis of the cluster of 50% of retrograde stars in the FG such that 75% of stars in the FG rotate on prograde orbits.

The FG system extends to the cluster’s tidal radius and has a ratio of the 3D half-mass radius to the tidal radius equal to about 0.085. The SG subsystem is initially centrally concentrated, and the ratio of the FG to the SG maximum (half-mass) radii is initially equal to 5 (5.93). Stars are removed from a simulation once they pass beyond a radius equal to twice the cluster’s Jacobi radius.

In order to explore the effects of the external tidal field and its torque on the internal rotational kinematics of the multiple stellar populations, we ran a set of simulations with different orientations of the cluster’s initial total internal angular momentum vector relative to the orbital angular momentum vector associated with the cluster’s orbit around the host galaxy (see Tiongco et al. 2019, 2022 for the first studies showing the complex kinematic features that may emerge from such a misalignment). Throughout this work, we utilised two parameters when describing the internal angular momentum. These are the internal angular momentum, which describes the angular momentum of individual stars in the cluster, and the total internal angular momentum, which describes the combined internal angular momentum of a group of stars in a specified region.

To explore the full range of dynamic possibilities resulting from this misalignment, we ran four different simulations corresponding to four different initial angles between the internal and orbital angular momentum vectors equal to 0, 45, 90, and 180 degrees (hereafter, we refer to those models as Delta0, Delta45, Delta90, and Delta180; see Table 1).

All the kinematic quantities presented in this study have been calculated in an inertial non-rotating reference frame. At any given time, all dynamical quantities have been calculated by combining five sequential snapshots of a system at times centred around the chosen time to improve statistics of the calculated quantities. The only exception to this is when looking at quantities calculated at a time of t = 0, where the initial conditions are used without combining multiple snapshots.

Our analysis of a cluster’s spatial and kinematic properties is performed in a cylindrical coordinate system defined as follows: the z-axis of this system is aligned with the direction of the total internal angular momentum vector of stars within the 2D half-mass radius of the aggregate snapshot; the other two coordinates are the standard radial and angular coordinates on the plane perpendicular to the z-axis, with the radius on this plane corresponding to the standard projected distance from the cluster’s centre; radial and tangential velocities used in the rest of the paper refer to the velocities in these two components. In particular, for the analysis of the internal angular momentum distribution and evolution, we calculated the rotational velocity profiles using the tangential velocities on the plane defined above. This calculation corresponds to the observational determination of rotational velocity using proper motion measurements adopting an ideal line of sight parallel to the cluster’s total internal angular momentum. The calculation of the total internal angular momentum for this purpose only includes stars within the cluster’s 2D half-mass radius; this was done to avoid the possible effect of stars in the outer regions of the system where the torque of the external tidal field may cause the internal angular momentum to be orientated differently from that of the stars in the inner regions (see Tiongco et al. 2019, 2022 for a discussion of this effect).

For the analysis of the time evolution of the kinematic properties, time is normalised to the half-mass relaxation timescale (see e.g. Heggie & Hut 2003), which is defined as $$t_{\mathrm{r}\mathrm{h}}=\frac{0.138M^{1/2}{r}_{\mathrm{h}}^{3/2}}{ln(0.02N)G^{1/2}<m>},$$(1)

where N is the number of stars in the system, M is the mass of the system, < m > is the mean stellar mass, r_h is the 3D half-mass radius, and G is the gravitational constant.

In some of our analysis, we fit the rotation curve radial profiles (see Sections 3.2 and 3.3) using the following analytical expression (Lynden-Bell 1967): $$V_{\mathrm{r}\mathrm{o}\mathrm{t}}=\frac{2R{a}_{\mathrm{r}\mathrm{o}\mathrm{t}}}{R_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}[1+(\frac{R}{R_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}}{)}^2]}$$(2)

Here, V_rot is the rotational velocity, R is the projected distance from the centre of the cluster, R_peak is the location of the peak velocity, and a_rot is the peak rotational velocity.

3 Results

3.1 Evolution of internal kinematics

We begin our analysis by exploring the simplest case, in which the total internal angular momentum vector of stars within the cluster’s Jacobi radius and the orbital angular momentum vector are initially aligned (referred to as the Delta0 model). It is important to note, however, that the results presented in this section are generally similar for all the models we studied.

The two populations have initially different spatial and kinematic properties, as dictated by the initial conditions. Fig. 1 shows how the mass density of the different generations evolves through the effects of two-body relaxation, and the spatial characteristics become more similar over time.

Figure 2 shows the radial profile of the rotational velocity. This figure shows that the rotational velocity decreases as a result of two-body relaxation as the system evolves and stars are ejected from the cluster, carrying away angular momentum. This is in agreement with previous results (Einsel & Spurzem 1999; Ernst et al. 2007; Kim et al. 2007; Tiongco et al. 2017). As the system evolves and the overall rotational velocity decreases, the rotational velocities of the two populations also become more similar. We point out that while the rotation curves of the two populations still differ after several relaxation times, the strength of this difference rapidly decreases, suggesting that even relatively small present-day differences would be indicative of much stronger primordial differences imprinted at the time of the cluster’s formation. We further discuss this point in Section 3.2 (see also the discussion of Dalessandro et al. 2024). It is interesting to point out that Berczik et al. (2025) also find that, in N-body models evolving on circular orbits similar to those used in this work, memory of the initial rotation may still be found in present-day clusters, while initial differences between the rotational velocities of multiple populations are erased more rapidly for clusters on eccentric orbits.

As the system evolves, the initial differences between the spatial and rotational properties of the two populations are gradually erased. However, the extent of the differences between the spatial and rotational profiles in the later stages of evolution differ. If, for example, we only consider stars within the inner two half-mass radii (similar to what is typically available in observational data), the spatial properties are seen to mix much earlier than the rotational properties. For example, by approximately 10 t_rh,i, the spatial properties in Fig. 1 are well mixed within the inner two half-mass radii, while the rotation profiles of the two generations in Fig. 2 are still distinct (see Dalessandro et al. 2019 for an observational analysis showing the loss of spatial differences between multiple populations over time). It is important to point out that the actual observational detection of differences between the rotational properties of FG and SG stars may be complicated by a variety of factors including the specific kinematic dataset available (line-of-sight and/or propermotion velocities), the radial range covered by the data, and the angle between the line of sight and the cluster’s rotation axis (see Tiongco et al. 2019 for a discussion of these points).

In addition to the rotational velocity, we also assessed the degree of anisotropy in the velocity space of the models, and its time evolution. Fig. 3 shows the anisotropy profiles at the same times as shown in Figures 1 and 2. To quantify the velocity anisotropy we used the ratio of tangential-velocity dispersion to the radial-velocity dispersion, σ_T/σ_R, where the tangentialand radial-velocity dispersions are calculated using cylindrical bins from velocities on the plane perpendicular to the cluster’s rotation axis.

As the system evolves, the outward migration of stars on more radial orbits causes the development of a significant radial anisotropy in the cluster’s outer regions; this is due to the progressive formation of a core-halo structure in collisional stellar systems due to two-body relaxation (see e.g. Heggie & Hut 2003). While for the initial conditions we adopted both populations develop radial anisotropy in the outer regions, the SG, which is initially more centrally concentrated, is characterised by a stronger radial anisotropy (see e.g. Bellini et al. 2015; Tiongco et al. 2019; Vesperini et al. 2021; Sollima 2021; Aros et al. 2025). The preferential loss of stars on radial orbits causes the clusters to eventually become tangentially anisotropic towards their outermost regions (see e.g. Giersz & Heggie 1997; Baumgardt & Makino 2003; Tiongco et al. 2016 for a similar effect in the context of single-population clusters). Differences between the velocity anisotropy of FG and SG stars, in agreement with those found in our simulations and in previous numerical investigations of multiple-population clusters (see e.g. Bellini et al. 2015; Vesperini et al. 2021; Sollima 2021; Aros et al. 2025), have also been found in several observational studies (see e.g. Bellini et al. 2015; Libralato et al. 2023; Cadelano et al. 2024; Dalessandro et al. 2024; Cordoni et al. 2025).

Fig. 1 Surface-mass-density profiles as a function of projected distance from cluster’s centre normalised by the projected half-mass radius (Rh) of the total cluster (see Sect. 2 for a complete definition). The profiles are shown at a few representative times in the Delta0 model. All densities are normalised by the central bin of the total surface-mass density (FG + SG). Time is normalised by the initial half-mass relaxation time of the cluster, trh,i. The difference between the profiles of the two generations is gradually erased as the system evolves.

Fig. 2 Radial profiles of rotational velocity are shown at representative times for the Delta0 model. The rotational velocity is normalised by the central velocity dispersion (inner 1% of stars) of all stars (FG + SG) (see Sect. 2 for details about the rotational velocity calculation). The radius is normalised by the projected half-mass radius (Rh) of all stars (FG +SG). Time is normalised by the initial half-mass relaxation time of the cluster, trh,i.

3.2 Dependence of kinematic properties on the stellar mass

We continue to focus the presentation of our results on the Delta0 model since, for this aspect of our study, we also find the results to be independent of the angle between the cluster’s initial rotation axis and the cluster’s orbital plane. Fig. 4 shows the rotational velocity profiles for different stellar masses for both the FG and SG stars. Within the individual generations, mass groups initially start out with the same rotation profile; however, as the system evolves, both populations develop a relationship between rotational velocity and stellar mass, where high-mass stars tend to rotate more rapidly than low-mass stars.

While two-body relaxation causes the overall rotation to decrease for each mass bin, the lowest mass bin loses rotational velocity more efficiently. This result agrees with the recent results of Livernois et al. (2022), which finds the same behaviour in the context of single-population clusters (see also Kim et al. 2004; Hong et al. 2013; Livernois et al. 2021; see Scalco et al. 2023 for the first observational evidence of this trend).

Figure 4 also shows that the difference between generations is erased at different rates for different mass groups. While higher mass stars maintain a higher overall rotational velocity compared to stars with lower masses regardless of generation, the difference between the generations themselves is erased more efficiently for higher mass stars. The final panel of Fig. 4, for example, shows that at 20 t_rh,i the rotational velocity of the two higher mass groups is nearly identical for FG and SG stars, while the lowest mass group is still characterised by discernible differences.

To further illustrate the differences between the rotational kinematics of FG and SG stars and their dependence on the stellar mass, we show, in Fig. 5, the time evolution of the difference between the FG and SG peak rotational velocities. As discussed in Section 2, we utilised best-fit lines for the rotational velocity profiles to give a clearer view of the trends. This figure clearly shows the gradual decrease of the differences between the rotation of the FG and the SG stars and how the initial differences are erased more rapidly for high-mass stars. The dependence of the strength of the rotational differences on the stellar mass implies that their detection in observational studies focused on high-mass stars (see e.g. Cordoni et al. 2020a,b; Martens et al. 2023; Dalessandro et al. 2024; Leitinger et al. 2025) may not reveal significant rotational differences in clusters where such differences could still be present in low-mass stars. We interpret this trend as the result of the interplay between the effects of mass segregation causing massive stars to migrate towards the cluster’s inner regions, where shorter values of the local relaxation timescale lead to a more rapid dynamical mixing of the two populations.

It is also interesting to note the rapid change in the kinematic differences at the beginning of the simulation. This suggests that caution is necessary when making claims about the initial kinematic properties of clusters based on present-day kinematics, as even in dynamically young clusters the strength of the presentday kinematic differences may be significantly smaller than the initial differences imprinted at the time of cluster formation (see also Tiongco et al. 2019 for a discussion and cautionary remarks concerning the observational detection of differences between the rotation of the FG and the SG populations). This is true for each mass regime, even in low-mass stars where the kinematic differences persist for the longest amount of time.

Finally, in Fig. 6 we show the radial profile of the velocity anisotropy for different mass bins. For all the mass bins, the SG star is characterised by a stronger anisotropy than the FG star, and within each stellar population we do not find significant dependence of the degree of anisotropy on stellar mass.

Fig. 3 Radial profiles of velocity anisotropy are shown at representative times in Delta0 model. We define anisotropy to be the ratio of tangential velocity dispersion over the radial velocity dispersion, σT/σR, with a completely isotropic system having a value of σT/σR = 1. Time and radius normalisations are the same as described in Fig. 2.

Fig. 4 Radial profiles of rotational velocity are shown at representative times for Delta0 model for different mass groups. The line of sight is the same as described in Sect. 2. The rotational velocity is normalised by the central velocity dispersion (inner 1% of stars) of all stars (FG + SG). Time and radius normalisations are the same as described in Fig. 2. Mass groups are split into low- (0.1-0.3 M⊙), intermediate-(0.3-0.7 M⊙), and high- (0.7-1 M⊙) mass bins.

Fig. 5 Time evolution of differences between maximum rotational velocities of FG and SG stars is shown for different mass groups. Velocities are normalised by the peak rotational velocity for all stars (FG + SG), computed at each time. Dots correspond to the peak velocity values taken directly from the rotation profiles. Lines represent the peak velocity values taken from best-fit rotation curves to the rotation profiles (as described in Eq. (2)). The time was limited to when simulations had approximately 10% of the initial stars remaining.

Fig. 6 Radial profiles of velocity anisotropy are shown at representative times in the Delta0 model for different mass groups. Anisotropy is defined as in Fig. 3. The line of sight is the same as described in Sect. 2. Time and radius normalisations are the same as described in Fig. 2.

3.3 The role of the orientation of cluster’s internal angular momentum relative to the orbital plane

We continue our analysis by discussing how the initial orientation of the cluster’s internal angular momentum relative to the plane of the cluster’s orbit around the centre of the host galaxy affects the evolution of its internal kinematic properties. In particular, we explore how the differences in the initial configuration may affect the interplay between internal rotation and the torque due to the tidal field of the host galaxy. This can possibly lead to a more complex radial variation of the rotation pattern as the system evolves towards (partial) tidal synchronisation (see e.g. Tiongco et al. 2018).

Figure 7 shows that the difference between the FG and SG peak rotational velocity is generally similar for the four models. However, the two models with initial rotation axis not parallel or antiparallel to the orbital angular momentum (Delta45, Delta90) are characterised by an extended phase (between about 15 t_rh,i and 40 t_rh,i), during which the difference between the FG and SG rotational velocity is approximately constant or decreases at a slower rate than that of the Delta0 and Delta180 models.

Figure 8 shows the peak rotational velocity difference between generations for individual mass groups in each model. The dependence on mass is similar regardless of our choice of model. The lowest mass bin most closely resembles the behaviour seen in Fig. 7. This is simply due to the fact that the majority of stars are those in the low-mass group, and they thus have a dominant effect on the time evolution shown in Fig. 7 for the entire system.

It is important to note, however, that despite the large FG-SG difference in the low-mass stars, these are the stars that are observationally most difficult to detect. Most observational studies of the internal kinematics of globular clusters are instead focused on upper-main-sequence and red-giant stars for which more accurate measurements of line-of-sight and/or propermotion velocities are possible (see e.g. Dalessandro et al. 2024). As shown by our results, the high-mass group is the one characterised by the smallest differences between the FG and SG rotation; larger differences between the rotational kinematics of multiple populations may thus be revealed by studies of the kinematics of low-mass stars.

Figure 8 also shows that all models are characterised by the rapid decrease of the initial differences between the rotational kinematics of FG and SG stars. As pointed out in Section 3.2, this implies that - regardless of the initial orientation of the internal angular momentum - even in relatively dynamically young clusters the present-day differences between the FG and SG rotational velocities are significantly smaller than the primordial ones.

Fig. 7 Time evolution of difference between maximum rotational velocities of FG and SG is shown for different models including all stellar masses. Velocities are normalised by the peak rotational velocity for all stars (FG + SG). Dots correspond to the peak velocity values taken directly from the rotation profiles. Lines represent the peak velocity values taken from best-fit rotation curves to the rotation profiles (as described in Eq. (2)). The time was limited to when simulations had approximately 10% of the initial stars remaining.

Fig. 8 Time evolution of difference between maximum rotational velocities of FG and SG is shown for different models and relevant mass groups. The low-mass bin (0.1-0.3 M⊙) is shown in blue, the intermediate-mass bin (0.3-0.7 M⊙) in orange, and the high-mass bin (0.7-1 M⊙) in green. Delta0 is denoted by a solid line, Delta45 by a dashed line, Delta90 by a dotted line, and Delta180 by a dash-dotted line. For the clarity of the figure, we only show best-fit lines to the data. The largest difference between generations can be seen in the lowest mass bin.

3.4 Angular momentum

We now consider the investigation of the time evolution of the internal angular momentum of FG and SG stars and how it is affected by the external tidal field of the host galaxy. We also explore how the initial internal rotation axis relative to the orbital angular momentum affects the cluster’s internal angular momentum over time. In our analysis, we explore the direction of the internal angular momentum through the use of δ, the angle measured from the model’s orbital angular momentum vector.

Figure 9 shows the time evolution of the δ angle for the total internal angular momentum of stars within the tidal radius for each simulation. The internal rotation of the Delta0 model begins aligned with the model’s orbital angular momentum vector and, therefore, the orientation of its internal rotation remains approximately constant during the entire simulation. The internal rotation of the Delta180 model is counter-aligned with the orbital rotation vector and also does not vary significantly for much of its lifetime. The Delta45 and Delta90 models, however, have a large change in their orientation over time due to interactions with the tidal field. The tidal field acts to align the internal angular momentum with the orbital angular momentum (see also Tiongco et al. 2018, 2019 for studies of the evolution of the orientation of a cluster’s internal rotation axis in single-population and multiple-population clusters). Hereafter, we investigate the Delta45 model more closely as a reference case.

Figure 10 shows the time evolution of the total internal angular momentum orientation for FG and SG stars within the tidal radius of the Delta45 model. The stars in the cluster interact with the tidal field, causing the total internal angular momentum of the system to align with the orbital angular momentum by the end of the simulation. After a small initial change in the direction of the FG stars’ total internal angular momentum, likely due to the effect of the tidal field on the more extended FG population, the two populations are characterised by similar orientations and evolution towards alignment. Figures 9 and 10 also show small nutations in the total internal angular momentum (see e.g. Tiongco et al. 2018).

The effect of the tidal field can be seen more clearly in Fig. 11, which shows the evolution of total internal angular momentum orientation for different radial regions in the Delta45 model. In all cases (FG, SG, and FG + SG), the stars in the outer regions begin to align their internal angular momentum with the orbital angular momentum earlier than stars in more central areas of the cluster. The same behaviour is present in both the Delta45 and Delta90 models.

As the system evolves and the two populations dynamically mix, the difference between the direction of their total internal angular momentum is also erased. Overall, the differences between the orientation of the FG and the SG systems are always small and likely difficult to detect observationally. Our results are consistent with those of Tiongco et al. (2019) who also found that the total internal angular momentum of the FG and SG initially evolve at different rates, but eventually align with the orbital angular momentum and each other.

We also explored the angular momentum of stars of different masses. Fig. 12 shows the orientation of the total internal angular momentum for stars in different mass groups within the tidal radius of the Delta45 model. Low-mass stars begin to change their δ faster than higher mass stars. We interpret this as the result of the preferential migration of low-mass stars towards the outer regions of the cluster, where the tidal field can more easily affect their orbital properties.

We find small differences between the direction of internal angular momentum of FG and SG stars in the different stellar mass groups. These differences are present in all cases that do not begin with their internal and orbital angular momentum vectors aligned, but the differences are strongest in the Delta90 and Delta180 models. Fig. 13 displays the orientation of the total internal angular momentum for different stellar mass groups within the tidal radius of the Delta90 and Delta180 models. There is a small but consistent difference between the direction of total internal angular momentum of the two generations in the Delta90 model for most of the system’s evolution. Delta180 shows the largest difference in the direction of total internal angular momentum between the generations; however, this is only for low-mass stars and only at later dynamical times. Similarly to the behaviour seen in Fig. 12, Fig. 13 also shows that low-mass stars begin to align with the direction of the orbital angular momentum faster than stars of higher masses.

It is interesting to note that, in each simulation, the generations are spatially mixed by the time the model’s direction of total internal angular momentum undergoes significant change (~40 t_rh,i). This suggests that the difference in rotational velocities between generations may be a driving factor behind the delay in alignment of the different generations.

The SG stars’ higher rotational velocity may cause the SG to take longer to align with the direction of orbital angular momentum. We stress again that while for many of the N-body models presented in this study there is a difference between the direction of total internal angular momentum of the two generations, the difference is small and would likely be difficult to detect observationally.

The results of the analysis presented thus far focus on the total internal angular momentum of the stars in the models. In order to gain further insight into the models’ kinematics, we examined the distribution of the orientation of internal angular momentum for the individual stars. Figures 14, 15, 16, and 17 utilise a method described in Irodotou & Thomas (2021) to create maps of the internal angular momentum of stars of different generations. Each map shows the distribution of the orientation of the individual stars’ internal angular momentum at a given time. The maps present a 2D projection of the histogram of the internal angular momentum distribution with brighter yellow areas corresponding to angles of angular momentum orientation with a larger fraction of stars. The orbital angular momentum vector is aligned with the bottom of the maps.

The centre of the map is aligned with the direction of the galactic centre. The first row of these figures demonstrates that the density of the direction of internal angular momentum when considering all stars is initially dictated by the highly rotating SG. The density of the direction of the internal angular momentum becomes more diffuse as the system evolves and the generations become spatially and kinematically mixed. The density of orientations for the two generations appears nearly identical after approximately 40 t_rh,i.

Figure 14 shows the difference in the spread of orientation of the internal angular momentum between generations. Nonetheless, the location of the peak in the distribution of the orientations stays constant due to the internal angular momentum already being aligned with the direction of the orbital angular momentum. Figures 15, 16, and 17 show the evolution of the distribution of the orientation of the internal angular momentum within the system as well as the difference between the distribution of the internal angular momentum of the two generations. Each model shows that there is an initially large difference in the distribution of the orientation of internal angular momentum between the stars of the two generations. The disc-like, high-rotation structure of the SG stars causes the distribution of orientations to be more peaked around the initial-angularmomentum orientation, while the more diffuse, lower rotation FG stars has a larger spread in the direction of internal angular momentum for stars within the generation. That difference in the density of orientations is gradually erased as the FG and SG subsystems mix. The Delta45 and Delta90 models clearly display the effects of the tidal field on the orientation of internal rotation as the peak density and overall distribution of internal angular momentum orientation moves to align with the direction of the orbital angular momentum at the bottom of the maps. In this way, these maps present another view of what is shown in Fig. 9.

While the Delta0, Delta45, and Delta90 models all end with the orientation of their internal angular momentum aligned with the orbital angular momentum, the same is not true of the Delta180 model. At the end of the simulation, the orientation of the internal angular momentum is still counter-aligned with the direction of the orbital angular momentum. The final time step of Fig. 17, however, shows that there is some diffusion away from this orientation. High-density regions on these maps at 80 t_rh,i show that stars have begun to change their orientations and move towards aligning with the direction of orbital angular momentum.

This connects to what is shown in Fig. 13, where the total internal angular momentum vector of low-mass stars within the tidal radius has largely aligned with the orbital angular momentum. While for many of the stars in the Delta180 model the direction of the internal angular momentum is still counteraligned with the direction of the orbital angular momentum, low-mass stars largely have the direction of their internal angular momentum aligned with the orbital angular momentum by the end of the simulation.

We can further explore this in Fig. 18. Here, we investigate the distribution of the direction of internal angular momentum for different mass bins in the Delta180 model. It can clearly be seen that the low-mass stars are aligned more closely with the direction of the orbital angular momentum, while the intermediate- and high-mass stars maintain an orientation similar to that of the initial conditions.

Finally, we utilise these maps to show an alternative view of what is displayed in Fig. 11. Fig. 19 shows that the peak density of the orientation of internal angular momentum, as well as the overall distribution of orientations, becomes increasingly aligned with the direction of the orbital angular momentum towards the outer regions of the cluster over time. Again, this is due to the effects of the interaction with the tidal field being stronger towards the outer regions of the system.

Fig. 9 Time evolution of orientation of total internal angular momentum vector of stars within cluster’s tidal radius for different models. The total (FG + SG) internal angular momentum is shown for each model. Time is normalised by the initial half-mass relaxation time. The Delta0, Delta45, and Delta90 models each end with the orientation of their total internal angular momentum aligned largely with the orbital angular momentum at zero degrees. The Delta180 model does not align with the orbital angular momentum by the end of the simulation.

Fig. 10 Time evolution of orientation of total internal angular momentum vectors of different subsystems including stars within the cluster’s tidal radius of the Delta45 model. Time is normalised by the initial halfmass relaxation time.

Fig. 11 Time evolution of orientation of total internal angular momentum vectors for different radial regions is shown for the different subsystems in the Delta45 model. Time is normalised by the initial halfmass relaxation time. Radial regions are set by the Lagrangian radii of all stars (FG + SG), regardless of subsystem. In each subsystem, the total internal angular momentum of stars in the outer regions begins to align with the orbital angular momentum faster than stars in more central regions.

Fig. 12 Time evolution of orientation of total internal angular momentum vectors for all stars (FG + SG) in Delta45 model is shown for different mass bins. Time is normalised to the initial half-mass relaxation time. Low-mass stars align with the orbital angular momentum more rapidly than higher mass stars.

Fig. 13 Time evolution of orientation of total internal angular momentum vectors for Delta90 (top) and Delta180 (bottom) models is shown for different mass bins and generations. Time is normalised to the initial half-mass relaxation time. Low-mass stars align with the orbital angular momentum more rapidly than higher mass stars in both models.

Fig. 14 Internal angular momentum maps are shown for FG + SG (left), FG (middle), and SG (right) for the Delta0 model. The colour map corresponds to the concentration of stars whose internal angular momentum points in a given direction. Higher concentrations are represented by a yellow colouration, while lower concentrations are represented by a black colouration. A top-hat smoothing was applied to the density maps. Rows show representative time steps throughout the simulation. The orbital angular momentum vector is aligned with the bottom of the maps.

Fig. 15 Similar to Fig. 14, but for the Delta45 model.

Fig. 16 Similar to Fig. 14, but for the Delta90 model.

Fig. 17 Similar to Fig. 14, but for the Delta180 model. Unlike the other models, the peak direction of internal angular momentum does not fully align with the direction of the orbital angular momentum by the end of the simulation for the FG, SG, or for FG + SG.

Fig. 18 Similar to Fig. 17, but for individual mass bins of all stars (FG + SG) in the Delta180 model. The density of the orientation of internal angular momentum for low-mass stars is much closer to aligning with the direction of the orbital angular momentum than that of the intermediate- or high-mass stars.

Fig. 19 Angular momentum maps are shown for FG + SG (top), FG (middle), and SG (bottom) for different radial regions in the Delta45 model at 60 trh,i. The internal angular momentum of stars in the outer regions are shown to align more closely with the direction of the orbital angular momentum.

4 Conclusions

In this work, we investigated the dynamical evolution of rotating multiple-population globular clusters. Our study is based on a suite of N-body simulations and addresses a number of questions concerning the long-term evolution of the differences between the dynamical properties of FG and SG stars imprinted at the time of their formation. Informed by a number of studies of the formation of multiple stellar populations in rotating clusters (see e.g. Bekki 2010, 2011; Lacchin et al. 2022), our simulations start with clusters characterised by a rapidly rotating and flattened SG subsystem concentrated in the inner regions of a slowly rotating, spherical FG cluster. In order to investigate the interplay between the effects of the external tidal field on the clusters’ internal rotation, we studied the evolution of systems starting with different initial inclinations of the internal angular momentum relative to the cluster’s orbital plane.

The results of our study are as follows.

• We explored the models’ dynamical evolution through the evolution of their mass density, rotational velocity, and velocity anisotropy profiles. By examining the density profiles, we see that the spatial differences between generations are greatly diminished within a few initial half-mass relaxation times (Fig. 1). We find that, even after several initial half-mass relaxation times, there is still a difference in the rotational velocities of the two generations (Fig. 2). However, it is important to emphasise that initial differences between the rotation of FG and SG stars rapidly decrease, and, even in dynamically young clusters with ages equal to just a few halfmass relaxation times, these differences may be significantly smaller than those imprinted at the time of the clusters’ formation. This is the case for all models (Fig. 7) and for stars with different masses (Figures 5 and 8). This suggests that caution is necessary in drawing conclusions concerning the primordial rotation differences and their extent based on present-day kinematics, even when considering relatively dynamically young clusters (see also Tiongco et al. 2019 and Dalessandro et al. 2024 for further discussion of this point.)

  In agreement with previous investigations, our analysis also reveals differences between the anisotropy of the velocity distribution of the FG and SG populations. SG stars are characterised by a stronger velocity anisotropy than that of FG stars. Although these differences also decrease during the clusters’ long-term evolution, they tend to persist longer than those between FG and SG rotational kinematics. Observational studies with wide fields spanning a broad range of clustercentric differences are more likely to detect differences between the FG and SG spatial and kinematic differences. Studies limited to the inner regions may fail to reveal dynamical differences or may reveal much smaller differences than those present in the intermediate and outer regions.

• By studying the time evolution of the rotation curves for different stellar masses and the resulting peak velocities, we find a dependence of rotational velocity on stellar mass. Higher mass stars maintain a higher peak rotational velocity than low-mass stars. We also see that the kinematic differences between generations are erased more efficiently for high-mass stars as compared to low-mass stars (see Figures 4 and 5). We therefore find that the inclusion of low-mass stars would improve the chances of detecting kinematic differences observationally.

• By exploring the interaction of the models’ internal rotation with the galactic tidal field, we find that models whose internal rotation is not initially aligned or counter-aligned with the direction of their orbital angular momentum may maintain differences between the rotational kinematics of FG and SG stars for a longer time (Fig. 7). These differences between generations are present predominantly in low-mass stars (Fig. 8). Our analysis reveals, however, that these differences may be small and difficult to detect observationally. The observational detection may be hindered by various factors such as errors in velocity measurements, limited statistical samples, possible misclassification of FG and SG stars (see also Tiongco et al. 2018 for possible problems due to projection effects and observations of clusters along different lines of sight).

• We explored the interaction of the model’s internal angular momentum and the galactic tidal field. We find that, in agreement with previous studies, the tidal field can cause the internal angular momentum of globular clusters to align with the direction of their orbital angular momentum over time. Models Delta45 and Delta90 largely align their internal angular momentum with the direction of their orbital angular momentum by the end of our simulations. The direction of internal angular momentum for Delta0, which begins with its internal and orbital angular momentum aligned, shows little variation throughout the simulation (see Figures 9, 10, 14, 15, and 16). The direction of internal angular momentum for the Delta180 model, however, remains largely counteraligned with the direction of the orbital angular momentum even at the end of the simulation (Fig. 17).

• We also explored the mass dependence of the interaction between the cluster’s internal angular momentum and the galactic tidal field. In our analysis, we determine that low-mass stars align their internal angular momentum with the direction of the cluster’s orbital angular momentum more rapidly than higher-mass stars. In addition, while the internal angular momentum of the Delta180 model does not align with the direction of the orbital angular momentum by the end of the simulation when considering stars of all masses, the internal angular momentum of the low-mass stars of that model does (see Figures 12, 13, and 18). We interpret this as the result of the preferential migration of low-mass stars towards the model’s outer regions, where the effects of the tidal field are stronger.

• By exploring the orientation of the models’ internal angular momenta in different radial regions, we find - in agreement with previous studies - that there is a radial dependence of the direction of internal angular momentum within a cluster. The galactic tidal field more strongly affects stars at larger radii, leading those stars to align the direction of their internal angular momentum with the direction of the orbital angular momentum faster than stars closer to the centre of the cluster (Figures 11 and 19).

Acknowledgements

E.V. acknowledges support from NSF grant AST-2009193. E.V. acknowledges also support from the John and A-Lan Reynolds Faculty Research Fund. E.D. acknowledges financial support from the INAF Data Analysis Research Grant (ref. E. Dalessandro) of the “Bando AstrofisicaFondamentale 2024”. A.L.V. acknowledges support from a UKRI Future Leaders Fellowship (MR/X011097/1). We thank the referee for a thorough review of this paper.

References

All Tables

All Figures

Fig. 1 Surface-mass-density profiles as a function of projected distance from cluster’s centre normalised by the projected half-mass radius (Rh) of the total cluster (see Sect. 2 for a complete definition). The profiles are shown at a few representative times in the Delta0 model. All densities are normalised by the central bin of the total surface-mass density (FG + SG). Time is normalised by the initial half-mass relaxation time of the cluster, trh,i. The difference between the profiles of the two generations is gradually erased as the system evolves.

In the text

Fig. 2 Radial profiles of rotational velocity are shown at representative times for the Delta0 model. The rotational velocity is normalised by the central velocity dispersion (inner 1% of stars) of all stars (FG + SG) (see Sect. 2 for details about the rotational velocity calculation). The radius is normalised by the projected half-mass radius (Rh) of all stars (FG +SG). Time is normalised by the initial half-mass relaxation time of the cluster, trh,i.

In the text

	Fig. 3 Radial profiles of velocity anisotropy are shown at representative times in Delta0 model. We define anisotropy to be the ratio of tangential velocity dispersion over the radial velocity dispersion, σT/σR, with a completely isotropic system having a value of σT/σR = 1. Time and radius normalisations are the same as described in Fig. 2.
In the text

Fig. 4 Radial profiles of rotational velocity are shown at representative times for Delta0 model for different mass groups. The line of sight is the same as described in Sect. 2. The rotational velocity is normalised by the central velocity dispersion (inner 1% of stars) of all stars (FG + SG). Time and radius normalisations are the same as described in Fig. 2. Mass groups are split into low- (0.1-0.3 M⊙), intermediate-(0.3-0.7 M⊙), and high- (0.7-1 M⊙) mass bins.

In the text

Fig. 5 Time evolution of differences between maximum rotational velocities of FG and SG stars is shown for different mass groups. Velocities are normalised by the peak rotational velocity for all stars (FG + SG), computed at each time. Dots correspond to the peak velocity values taken directly from the rotation profiles. Lines represent the peak velocity values taken from best-fit rotation curves to the rotation profiles (as described in Eq. (2)). The time was limited to when simulations had approximately 10% of the initial stars remaining.

In the text

	Fig. 6 Radial profiles of velocity anisotropy are shown at representative times in the Delta0 model for different mass groups. Anisotropy is defined as in Fig. 3. The line of sight is the same as described in Sect. 2. Time and radius normalisations are the same as described in Fig. 2.
In the text

Fig. 7 Time evolution of difference between maximum rotational velocities of FG and SG is shown for different models including all stellar masses. Velocities are normalised by the peak rotational velocity for all stars (FG + SG). Dots correspond to the peak velocity values taken directly from the rotation profiles. Lines represent the peak velocity values taken from best-fit rotation curves to the rotation profiles (as described in Eq. (2)). The time was limited to when simulations had approximately 10% of the initial stars remaining.

In the text

Fig. 8 Time evolution of difference between maximum rotational velocities of FG and SG is shown for different models and relevant mass groups. The low-mass bin (0.1-0.3 M⊙) is shown in blue, the intermediate-mass bin (0.3-0.7 M⊙) in orange, and the high-mass bin (0.7-1 M⊙) in green. Delta0 is denoted by a solid line, Delta45 by a dashed line, Delta90 by a dotted line, and Delta180 by a dash-dotted line. For the clarity of the figure, we only show best-fit lines to the data. The largest difference between generations can be seen in the lowest mass bin.

In the text

Fig. 9 Time evolution of orientation of total internal angular momentum vector of stars within cluster’s tidal radius for different models. The total (FG + SG) internal angular momentum is shown for each model. Time is normalised by the initial half-mass relaxation time. The Delta0, Delta45, and Delta90 models each end with the orientation of their total internal angular momentum aligned largely with the orbital angular momentum at zero degrees. The Delta180 model does not align with the orbital angular momentum by the end of the simulation.

In the text

	Fig. 10 Time evolution of orientation of total internal angular momentum vectors of different subsystems including stars within the cluster’s tidal radius of the Delta45 model. Time is normalised by the initial halfmass relaxation time.
In the text

Fig. 11 Time evolution of orientation of total internal angular momentum vectors for different radial regions is shown for the different subsystems in the Delta45 model. Time is normalised by the initial halfmass relaxation time. Radial regions are set by the Lagrangian radii of all stars (FG + SG), regardless of subsystem. In each subsystem, the total internal angular momentum of stars in the outer regions begins to align with the orbital angular momentum faster than stars in more central regions.

In the text

	Fig. 12 Time evolution of orientation of total internal angular momentum vectors for all stars (FG + SG) in Delta45 model is shown for different mass bins. Time is normalised to the initial half-mass relaxation time. Low-mass stars align with the orbital angular momentum more rapidly than higher mass stars.
In the text

	Fig. 13 Time evolution of orientation of total internal angular momentum vectors for Delta90 (top) and Delta180 (bottom) models is shown for different mass bins and generations. Time is normalised to the initial half-mass relaxation time. Low-mass stars align with the orbital angular momentum more rapidly than higher mass stars in both models.
In the text

Fig. 14 Internal angular momentum maps are shown for FG + SG (left), FG (middle), and SG (right) for the Delta0 model. The colour map corresponds to the concentration of stars whose internal angular momentum points in a given direction. Higher concentrations are represented by a yellow colouration, while lower concentrations are represented by a black colouration. A top-hat smoothing was applied to the density maps. Rows show representative time steps throughout the simulation. The orbital angular momentum vector is aligned with the bottom of the maps.

In the text

	Fig. 15 Similar to Fig. 14, but for the Delta45 model.
In the text

	Fig. 16 Similar to Fig. 14, but for the Delta90 model.
In the text

	Fig. 17 Similar to Fig. 14, but for the Delta180 model. Unlike the other models, the peak direction of internal angular momentum does not fully align with the direction of the orbital angular momentum by the end of the simulation for the FG, SG, or for FG + SG.
In the text

	Fig. 18 Similar to Fig. 17, but for individual mass bins of all stars (FG + SG) in the Delta180 model. The density of the orientation of internal angular momentum for low-mass stars is much closer to aligning with the direction of the orbital angular momentum than that of the intermediate- or high-mass stars.
In the text

	Fig. 19 Angular momentum maps are shown for FG + SG (top), FG (middle), and SG (bottom) for different radial regions in the Delta45 model at 60 trh,i. The internal angular momentum of stars in the outer regions are shown to align more closely with the direction of the orbital angular momentum.
In the text