© The Authors 2025

1 Introduction

A well-established feature of globular clusters (GCs) is that the vast majority host multiple stellar populations (MPs), that are identifiable due to differences in the light-element abundances between stars, shown both in spectroscopy and photometry (see reviews by Bastian & Lardo 2018 and Milone & Marino 2022). In general, the MPs consist of a ‘primordial’ population of stars (‘P1’ stars) with light-element abundance patterns similar to surrounding field stars of the same metallicity. There is also at least one ‘enriched’ population of stars (‘P2’ stars), that are commonly enriched in elements, such as N, Na, and Al, but depleted in C, O, and sometimes Mg, in comparison to P1 stars. It is currently unknown which initial conditions are responsible for the birth of these MPs, that occurs during the early stages of GC formation. Many formation theories suggest a generational gap between the populations, in which P1 stars were born first, while P2 stars were created from enriched material at some stage afterwards – that is, the Asymptotic Giant Branch (AGB) scenario (Cottrell & Da Costa 1981; Dantona et al. 1983; Ventura et al. 2001), fast rotating massive stars (FRMS) (Decressin et al. 2007b,a), super-massive stars (SMS) (Denissenkov & Hartwick 2014; Gieles et al. 2018), or massive interacting binaries (de Mink et al. 2009; Renzini et al. 2022). Different formation mechanisms imply different kinematic imprints on the MPs, but as most of the GCs that exist today within our Milky Way are very old (>11 Gyr), long-term dynamical evolution most likely has erased traces of the initial conditions of their birth. However, GCs that have not undergone significant dynamical evolution during their lifetime can be described as ‘dynamically young’ and are expected to retain their initial conditions, particularly in the outer regions (Vesperini et al. 2013). In the absence of direct observations of MPs as they form in high redshift galaxies, young massive clusters and dynamically young GCs serve as the next best proxies for studying the initial conditions of the formation of MPs.

The dynamics of GCs can be separated into two aspects: a wider view of the overall motion of GCs within the host galaxy and a closer view of the internal dynamics as determined by the interactions of individual stars within the cluster. Within the host galaxy, the orbital information of GCs can provide insight into whether specific GCs formed in situ or were accreted into the host galaxy (Massari et al. 2019; Belokurov & Kravtsov 2024), with dynamical modelling providing insight into the effect of tidal interactions on the clusters (Sollima & Baumgardt 2017). Within the cluster itself, gravitational interactions influence the motions of the stars and two-body relaxation causes the cluster to become isotropic over time. Each cluster can also exhibit varying levels of rotation, in which the stars exhibit ordered motions. This present-day rotation can also be indicative of the rotation of the cluster at birth, depending on the degree of dynamical evolution the cluster has undergone within its lifetime (Sollima et al. 2020). Long-term dynamical evolution within GCs causes the initial orbital properties of the multiple stellar populations to mix in phase space (Decressin et al. 2008; Vesperini et al. 2013; Hénault-Brunet et al. 2015). The two-body relaxation timescale is inversely proportional to the local density, so since the density of a cluster increases towards the centre (Spitzer 1987), the relaxation timescale is longer in the outer regions. This means that the initial configurations in phase-space ultimately become mixed on a faster timescale in the central regions, whereas the outer regions are more likely to be representative of the initial configurations. It is therefore crucial to extend any kinematic analysis into the outer regions of GCs in order to obtain the best chance of observing any remaining initial conditions.

Some formation theories imply that MPs within GCs have different ages, initial radial configurations, or kinematics. Studying the rotation of stars within a GC, and in particular any differences in the rotational signatures of the multiple populations, can help to distinguish between the possible formation scenarios. For formation theories that assume (1) that P1 stars are located primarily in the outer regions and (2) that ~90% of P1 stars are removed from a cluster within its lifetime (i.e. as expected in the AGB scenario; Decressin et al. 2007b; D’Ercole et al. 2008; Schaerer & Charbonnel 2011; Cabrera-Ziri et al. 2015), the result is that a large amount of angular momentum would be removed with them. In this case, one would expect P1 stars (if they are predominantly in the outer regions) to have lower angular momentum (Hénault-Brunet et al. 2015) than the enriched stars. If P2 stars form with high angular momentum in the core of a cluster with low angular momentum, these P2 stars are predicted to have higher initial rotational amplitudes than the P1 stars located further out, due to the conservation of angular momentum. Observationally, Dalessandro et al. (2024) found this to be the case after analysing the 3D rotation profiles of 16 Galactic GCs, concluding that P2 stars rotate faster than P1 stars for the majority of clusters in their sample. Cordero et al. (2017) also discovered the same result for NGC 6205 after analysing the spectra of 111 giant stars. They found that the P2 stars in their sample have higher rotational amplitudes than the P1 stars, and found differences of 0–45 degrees between the spin axes of the populations. In addition, P2 stars were found to rotate faster than P1 stars in NGC 5272 (Szigeti et al. 2021), NGC 6093 (Kamann et al. 2020) and NGC 6362 (Dalessandro et al. 2021).

Similarly, Martens et al. (2023) analysed the rotational differences between populations in 25 Milky Way GCs using a combination of HST photometry and MUSE spectroscopy and found that only two GCs in their sample had P2 stars rotating faster than P1 stars; however, they also discovered one GC in their sample that showed the opposite. This was also observed by Bellini et al. (2018) using HST photometry of ω Centauri, where they discovered P1 stars with higher rotational amplitudes than P2 stars. A possible explanation for this is that if P2 stars – regardless of where they formed in the cluster – initially passed through the cluster core, they would be on radial orbits with low angular momentum and exhibit lower rotational amplitudes (Hénault-Brunet 2015). In many GCs, no differences in the rotational amplitudes of their MPs have been observed, such as the clusters NGC 104 (Milone et al. 2018; Cordoni et al. 2020b), NGC 6121 and NGC 6838 (Cordoni et al. 2020b), NGC 6352 (Libralato et al. 2019), as well as NGC 7078 (Szigeti et al. 2021). Many formation theories predict differences between the rotational amplitudes of the MPs, but the absence of a difference has not been explicitly predicted by any formation theories besides the FRMS scenario (Krause et al. 2013), that predicts that P1 and P2 stars should share the same kinematics, since P2 stars are formed within the decretion discs of massive P1 stars.

Theoretically, anisotropy profiles should also provide insight into the possible initial conditions and evolution of GCs (Pavlík & Vesperini 2021, 2022; Pavlík et al. 2024). It is predicted that if a GC is born with tangential anisotropy, the relaxation processes are accelerated in comparison to GCs that were born isotropically or with radial anisotropy. It is also expected that the outer regions of a GC have a higher chance of retaining any initial anisotropy of a GC, with radial anisotropy expected for GCs that are tidally underfilled and tangential anisotropy expected for tidally filled GCs (Zocchi et al. 2016; Pavlík et al. 2024). Observationally, the outer regions of dynamically young GCs should therefore be ideal for determining clues about the initial configuration.

The anisotropy of MPs are also expected to provide insight into the kinematic evolution and potential initial conditions of their GCs according to N-body simulations (Tiongco et al. 2019; Vesperini et al. 2021). Theoretically, if P2 stars are more centrally concentrated and diffuse in the outer regions, then they should exhibit radial anisotropy in the outskirts, which has been observed in a handful of MW GCs (Richer et al. 2013; Bellini et al. 2015; Milone et al. 2018; Libralato et al. 2018). This is especially expected for clusters that begin compact in comparison to their tidal radii, as the expansion of the cluster forces stars into strongly radial orbits in the outer regions, until many of these highly radial stars escape the cluster (Tiongco et al. 2016, 2019). On the other hand, tangential anisotropy, or isotropy, is expected for clusters that have already filled their Roche lobe during formation and therefore have not expanded significantly in their lifetime.

The present-day signature of radial anisotropy in the outer regions is not specific to any particular formation theory, so it cannot help to narrow down the possible options. All generational scenarios predict that P2 stars are centrally concentrated during their formation, leading them to be scattered on wider radial orbits as the cluster evolves (Hénault-Brunet et al. 2015). In fact, in any scenario where the gas cools and collects in the centre of the cluster and P2 stars are formed, or the P2 stars are the first to traverse the cluster core on radial orbits, the resulting anisotropy profiles of both populations will end up being the same in the present day according to simulations by Hénault-Brunet et al. (2015).

The rotational orientation of GCs and their MPs can also be useful for comparisons between observations and formation theories. For the SMS scenario described in Gieles et al. (2018), it is theorised that P2 stars are formed from the winds of an SMS and therefore are not expected to show coherent motion relative to the P1 stars, nor are their rotational orientations expected to match those of the P1 stars. This is because the SMS is formed by stellar collisions, so its angular momentum increases through a random walk, causing the SMS to exhibit a random spin direction relative to the pre-existing P1 stars. Although the rotational velocity of the SMS can reach the order of ~100 km/s, this velocity is not inherited by the P2 stars, which instead formed from the stellar winds of the SMS, with the rotational velocity of the SMS becoming negligible once the wind reaches ~1 pc. Using 3D hydrodynamical simulations that follow the formation process of the AGB scenario, McKenzie & Bekki (2021) found that after formation, P2 stars will rotate faster than P1 stars and exhibit highly coherent rotation in a disc-like structure parallel to the plane of the parent galaxy. Similarly, the 3D hydrodynamical simulations of Lacchin et al. (2022), also based on the AGB formation scenario, result in P2 stars rotating faster than P1 stars, with both sets of simulations showing that the P1 and P2 stars rotate in phase with each other (i.e. their rotation peaks are aligned). However, these simulations differ in terms of the spatial concentration of the MPs after formation. At the end of the 370 Myr timeframe of the McKenzie & Bekki (2021) simulation, P1 stars are found to be centrally concentrated within the cluster, while Lacchin et al. (2022) found the opposite. Current simulations have yet to accurately match observations of MPs in terms of both the present day spatial distributions and kinematic signatures.

In this paper, we have homogeneously analysed 30 Galactic GCs and explored the differences in the dynamics of each cluster, as well as their MPs. Unveiling clues about the formation history of these GCs is a task best approached by investigating a large, diverse sample of resolved GCs, as the observable trends will have higher significance than if the focus is to thoroughly investigate one or two clusters. The GCs in our sample are diverse in terms of the range of masses, dynamical ages, metallicities and locations within the halo of the Milky Way, with every GC in our sample found to contain multiple stellar populations.

2 Observational data

We investigated 30 Galactic globular clusters with distances from the Sun between 1.85 ± 0.02 ≤ R_⊙ (kpc) ≤ 18.5 ± 0.18 (Baumgardt & Vasiliev 2021), masses between (3.78 ± 0.14) ⋅ 10⁴ ≤ M_⊙ ≤ (8.5 3 ± 0.05) ⋅ 10⁵ (Baumgardt & Hilker 2018)¹, and metallicities between –2.37 ≤ [Fe/H] ≤ –0.59 (Harris 2010). Our GC sample is listed in Table A.1 and contains the 28 GCs for which we identified multiple stellar populations in Leitinger et al. (2023), plus the addition of NGC 104 and NGC 6656 as discussed in 2.1. The positions of our clusters imposed on an artist’s conception image of the Milky Way are shown in Figure 1. The cluster parameters such as position, distance from the Galactic centre, proper motion of the cluster, etc. used throughout this paper were taken from the Galactic Globular Cluster Database Baumgardt et al. (2023b), unless specified otherwise in the text.

Fig. 1 Sample of 30 Galactic globular clusters in galactocentric coordinates plotted over an artist’s conception image of the Milky Way (R. Hurt: NASA/JPL-Caltech/SSC) using mw-plot. The size of the points are proportional to the mass of the clusters.

2.1 Photometric catalogues

In our first paper (Leitinger et al. 2023), we used space-based HST (Piotto et al. 2015; Nardiello et al. 2018) and groundbased (Stetson et al. 2019) photometric catalogues to obtain a wide-field view of 28 clusters and classified their multiple stellar populations using RGB stars. The normalised, cumulative radial distributions of the P1 and P2 stars were then calculated in order to determine which population was more centrally concentrated in the cluster. Quantitatively, this was expressed by the A⁺ parameter (Alessandrini et al. 2016; Dalessandro et al. 2019), with negative (positive) values indicating the P2 (P1) stars were more centrally concentrated, or a value close to zero indicating the populations were spatially mixed. Using the same method, we added NGC 104 and NGC 6656 to the sample and determined the A⁺ values and enriched star fractions P₂/P_tot, which are shown in Table A.1. Figures 15 and 16 of Leitinger et al. (2023) showed the combined A⁺ values as a function of both the dynamical age (age/relaxation time) and mass-loss ratio (M_current/M_initiai), so an updated version of these plots is described and included in Appendix B, in order to reflect the addition of NGC 104 and NGC 6656 and any changes to the global structural parameters of our sample.

2.2 Proper motions

In order to determine the kinematics of the globular clusters in the plane of the sky, we used proper motions from Gaia DR3 (Gaia Collaboration 2016, 2023) that cover mainly the outer regions of each cluster up to the tidal radius, combined with proper motions from the Hubble Space Telescope Atlases of Cluster Kinematics (HACKS) (Libralato et al. 2022) that covers the central regions of the clusters.

For each globular cluster, we searched for stars in Gaia DR3 out to the tidal radius (r_t). The list of candidate members found in Gaia DR3 was then cleaned of non-members following the method described in Section 2 of Vasiliev & Baumgardt (2021), as well as proper motion cleaning using the χ² test described by Equation 3 from Leitinger et al. (2023). By using the Gaia BP- RP and G magnitudes, we isolated stars in the resulting CMD that were one magnitude fainter than the turn-off (G_RGBturn-off) and performed N – σ clipping to keep only stars that belonged to the RGB, AGB, HB, and MS-TO of each cluster. This created a cleaned catalogue of Gaia DR3 proper motions, which we used for the remainder of the analysis.

Vasiliev & Baumgardt (2021) found that the proper motion uncertainties given in the Gaia DR3 catalogue (Gaia Collaboration 2016, 2023) are underestimated by 10–20%, particularly in the dense central regions of globular clusters. To correct for this, we first calculated the number density (Σ) of all stars and followed the method described in Vasiliev & Baumgardt (2021), applying their Equation (3) to the original Gaia proper motion errors in order to correct this underestimation. We selected the value of ${ϵ}_{\overline{\mu },\text{sys}}\simeq 0-0.02$ mas for each cluster based on Figure 2 of Vasiliev & Baumgardt (2021). Throughout this paper we have used these externally calibrated proper motion errors, ${ϵ}_{\overline{\mu },\text{ext}}$, rather than the errors given in the Gaia catalogue.

In the HACKS astro-photometric catalogue (Libralato et al. 2022), we removed stars with unreliable photometry and proper motions before using them in our kinematic analysis. To do this, we selected well-measured objects in the HACKS photometric catalogues for each cluster in our sample, following the method described in points (i) to (vi) in Section 4 of Libralato et al. (2022). We then selected reliable stars in the corresponding HACKS proper motion catalogues by removing all objects with $N_u^{\text{PM}}/N_f^{\text{PM}}$ < 0.8 to 0.9 (with the value dependent on each cluster), ${\chi }_{{\mu }_{\alpha }\cos \delta }^2,{\chi }_{{\mu }_{\delta }}^2$ > 1.25 to 1.5 and proper motion errors >0.5 mas yr^–1, as also detailed by Libralato et al. (2022).

For each cluster, stars with reliable photometry were matched with the corresponding cleaned proper motion catalogue in order to create CMDs using the combinations of available F606W and F814W filters. The resulting F 606W – F814W versus F 606W CMDs served as a means to isolate the RGB, AGB, HB, and MS-TO stars for the kinematic analysis. We implemented a cut to remove stars fainter than one magnitude below F606W_RGBtum-off, as well as removing blue straggler stars from the final sample, so as not to affect the velocity dispersions of the remaining stars. We also removed stars located too far from the RGB, AGB, HB, and MS-TO using iterative N – σ clipping along the main distribution of the stars.

The proper motions included in the HACKS catalogue are reported in a relative reference system in which the motion of a star is measured relative to the surrounding stars, unlike the absolute reference system used by the Gaia DR3 catalogue. The rotation of the clusters can therefore not be observed in HACKS, and it was not possible to perform a rotational analysis using the HACKS data. However, we did use the HACKS proper motions to calculate the velocity dispersion and anisotropy profiles of each cluster in Sections 3.3 and 3.4.

2.3 Line-of-sight velocities

We obtained line-of-sight (LOS) velocities from the homogeneous catalogue created by Baumgardt (2017) and Baumgardt & Hilker (2018), which included velocities derived from ESO FLAMES/UVES, Keck HIRES/DEIMOS spectra, and published literature values. We also used the additional line-of-sight velocities of globular cluster stars which Baumgardt et al. (2023a) compiled from the APOGEE DR17 (Abdurro’uf et al. 2022), Gaia DR3 (Gaia Collaboration 2016, 2023), Galah DR3 (Buder et al. 2021), LAMOST DR8 (Cui et al. 2012), MIKIS (Ferraro et al. 2018), and RAVE DR5 (Kunder et al. 2017) surveys. In total, this data set contains about 37 000 LOS velocity measurements for about 22 000 unique cluster stars in the 30 studied clusters, with an average of 1.7 LOS measurements per star. Additionally, MUSE LOS velocities from Martens et al. (2023) were used for stars in the inner regions of the clusters. We removed stars with line-of-sight velocity errors >5 km/s, and averaged the remaining velocities and associated errors for each star. The distribution of LOS velocities as a function of the projected distance from the centre of each cluster was then fitted with a linear function and cleaned of >3σ outliers using an iterative process.

We began with the same sample of GCs as included in Leitinger et al. (2023), with the addition of NGC 104 and NGC 6656. However, from this combined sample of GCs, three GCs did not have a sufficient (>100) number of stars in the LOS velocity catalogues to perform the rotation analysis described in Section 3.1 and were therefore not included in the multiple stellar population analysis of Section 3.2. These clusters were NGC 5053 (N_LOS = 70), NGC 6934 (N_LOS = 77), and NGC 6981 (N_LOS = 65), where the quoted N_LOS of each cluster refers to the full sample of stars with available LOS velocities before any cleaning was applied. Despite having insufficient LOS velocities to determine rotation, the number of proper motion measurements allowed us to include these clusters in the velocity dispersion and anisotropy analysis of Sections 3.3 through to 3.5.

We added archival MUSE observations of NGC 6101 from proposal ID 099.D-0824 (PI: Peuten), covering the inner ~3 × 3 arcmin of the cluster with eight pointings of the MUSE Wide Field Mode. The data were reduced with the standard MUSE pipeline (Weilbacher et al. 2020), resulting in a single reduced data cube for each of the pointings. Each cube contained the signal from four exposures of 600 s each, and small spatial dithers and derotator offsets of 90 degrees between them. Individual stellar spectra were extracted from the cubes using PAMPELMUSE (Kamann et al. 2013). The code uses an astrometric reference catalogue in order to determine the positions of resolvable stars and the point-spread function in the integral field data as a function of wavelength. This information is then used to optimally extract the spectra of the resolvable stars while accounting for blends between adjacent stars. The photometric reference catalogue required by PAMPELMUSE was taken from the ACS survey of Galactic globular clusters (Sarajedini et al. 2007; Anderson et al. 2008). The extracted spectra were processed with SPEXXY (Husser et al. 2016), which performed a full-spectrum fit of each input spectrum against the GLIB (Husser et al. 2013) synthetic spectral library. The required initial guesses for the surface gravity and effective temperature of each star were obtained by comparing the aforementioned photometry to an isochrone from the database of Bressan et al. (2012), selected to match the metallicity of NGC 6101. Initial guesses for the line-of-sight velocity of each spectrum were obtained by cross-correlating it with a GLIB template of matching stellar parameters. In total, the SPEXXY fits yielded 427 velocity measurements that passed our quality cuts.

We included new observations of the central part of the cluster NGC 5024. We used a mosaic of eight pointings with the MEGARA@GTC IFU. The data were reduced using the MEGARA Data Reduction Pipeline (MDRP, Pascual et al. 2022). The process began by subtracting the bias level present in the images. To achieve this, we utilised the bias calibration and the MegaraBiasImage task. The bias level varies slightly between the top and bottom of the image (by approximately 100 counts) due to the MEGARA CCD (an E2V CCD231-84) being read out through two diagonally opposite amplifiers. At this stage, defective pixels were automatically masked using the master_bpm.fits file provided by the MDRP. Next, the spectra were processed using the Trace and ModelMap tasks, that trace the fibre spectra from the halogen lamp images. The MDRP was then used to perform wavelength calibration using ThNe lamp images. Flat-field correction was carried out using the halogen lamp images. For each calibrated exposure, we computed the median of the images within each observing block. Additionally, the MDRP subtracts the mean sky spectrum, derived from the dedicated sky fibres, to produce the fully reduced spectra. The reduced MEGARA data of NGC 5024 consisted of raw-stacked spectra (RSS) files, with each one-dimensional spectrum containing the flux recorded by a single hexagonal fibre. Again, we extracted individual stellar spectra from the data using PAMPELMUSE (Kamann et al. 2013). To avoid having to resample the MEGARA data to a rectangular grid, PAMPELMUSE was modified so that it accepts the native (hexagonal) MEGARA sampling as an input format. As basis for the reference catalogue, we used the F814W photometry of NGC 5024 from the HACKS survey (Libralato et al. 2022). As many bright stars are missing in the catalogue, likely because they were saturated in the underlying HST images, we complemented the catalogue with Gaia DR3.

After extraction of the spectra, we determined stellar radial velocities with the help of the IRAF² fxcor task, using as templates the spectra of cool giant stars of a metallicity that is comparable to the cluster metallicities. We created the template spectra with the stellar synthesis programme SPECTRUM (Gray & Corbally 1994) using ATLAS9 stellar model atmospheres (Castelli & Kurucz 2003). For the later analysis, we only use stars with radial velocity errors <3 km/sec, magnitude errors <0.5, and signal-to-noise ratio S/N > 3.

3 Method

Our approach was to first calculate the rotation and anisotropy profiles of the 30 Galactic globular clusters as a whole, using the full sample of RGB, AGB, HB, and MS-TO stars available in the observational data described in Section 2. In order to then investigate the differences in 3D kinematics between the multiple stellar populations of each cluster, we cross-match the cleaned, observational data for each cluster with the photometrically classified multiple stellar populations of RGB stars presented in Leitinger et al. (2023), and completed the same analysis again on this smaller subset.

We analysed the 3D rotation of stars in each cluster by combining the cleaned proper motions from Gaia DR3 with the LOS velocity catalogues, determining the rotation axes and 3D total rotational amplitudes for all RGB, AGB, HB, and MS-TO stars available. Section 3.1 describes the analysis of this large sample of stars, for each cluster. Section 3.2 then describes the analysis of the RGB subset of this sample that were classified into multiple stellar populations.

Similarly, in order to analyse the anisotropy of each cluster, we combined the HACKS proper motions with Gaia DR3 proper motions for all RGB, AGB, HB, and MS-TO stars available.

Section 3.3 describes the velocity dispersion analysis and Section 3.4 describes the anisotropy analysis for this large sample. Then, Section 3.5 describes the analysis of the subset that was cross-matched with the photometric RGB sample of stars classified into multiple stellar populations.

Fig. 2 Number of stars with measured proper motions (indigo) and LOS velocities (orange) as a function of the projected distance from the cluster centre for NGC 6809. The median values are shown by dashed lines.

3.1 Rotation axes and total rotational amplitudes

We began with all cleaned RGB, AGB, HB, and MS-TO stars available in the Gaia DR3 and LOS catalogues for the clusters in our sample, except NGC 5053, NGC 6934, and NGC 6981, due to the low number of LOS velocities for those clusters.

In order to account for perspective rotation that can add an artificial rotation signal to the data, we corrected the proper motions and LOS velocities. A quick estimate shows that for our sample, perspective rotation is only important for NGC 3201 due to its large LOS velocity of R_V = 495.39 ± 0.06 km/s and close proximity to the observer at R_⊙ = 4.74 ± 0.04 kpc. We therefore only correct the motion of stars in NGC 3201. The method for correcting the velocities is described in van de Ven et al. (2006) and Wan et al. (2021) using Equations (1) and (2) for the components of the proper motion and Equation (3) for the LOS velocities. In these equations, ${\mu }_{x\prime }^{sys},{\mu }_{y\prime }^{sys}\text{\hspace{0.17em}and\hspace{0.17em}}{v}_{ɀ\prime }^{sys}$ are the systematic velocities, x′ and y′ are the projected coordinates as calculated by Equation (2) of van de Ven et al. (2006) and D is the distance to the cluster $${\mu }_{x\prime }^{\text{pr}}=-6.1363\times {10}^{-5}x\prime v_{ɀ\prime }^{sys}/D\text{\hspace{0.17em}mas}/\text{yr}$$(1) $${\mu }_{y\prime }^{\text{pr}}=-6.1363\times {10}^{-5}y\prime v_{ɀ\prime }^{sys}/D\text{\hspace{0.17em}mas}/\text{yr}$$(2) $$v_{ɀ\prime }^{\text{pr}}=1.3790\times {10}^{-3}\left(x\prime {\mu }_{x^{\prime }}^{sys}+y\prime {\mu }_{y\prime }^{sys}\right)D\text{\hspace{0.17em}km}/\text{s}.$$(3)

We compared the projected distances of stars with measured proper motions and radial velocities for all clusters, as shown in Figure 2 for NGC 6809, to ensure the median distances are within 100'' of each other for all clusters. This guarantees that both samples are taken at similar radii on average so that the proper motions and LOS velocities can directly be compared with each other. In this way, we avoided the requirement that each star be matched in both the proper motion and LOS catalogue, which would decrease the number of stars available for analysis. If the median distances of the stars with LOS measurements disagreed by more than 100″ from that of the stars with PM measurements, we cross-matched stars from the MUSE catalogue with Gaia DR3, keeping only the stars that could be identified in Gaia DR3. This cross-match allowed the medians to align for all clusters except NGC 104, which instead required a cross-match between the proper motions and LOS catalogues in order to ensure we were observing rotation at comparable radii. Due to the large number of observations for stars in NGC 104, we were still left with 1409 stars after cross-matching in this way.

When using Gaia data for globular clusters, it is good practice to adjust the celestial coordinates and proper motions to account for projection effects due to the curvature of the sky, which is especially necessary for large and nearby clusters. We used the standard orthographic projection shown in Equation 2 of Gaia Collaboration (2018) to transform the celestial coordinates and corresponding proper motions into Cartesian coordinates on a tangent plane. We then subtracted the average cluster proper motions (µ_α*,cluster, µ_δ,cluster) and LOS velocity (v_LOS,cluster) using the values determined in Baumgardt & Hilker (2018). Using the transformed coordinates and proper motions, we calculated the radial (v_rad) and tangential (v_tan) components of the proper motions, such that $${\text{v}}_{\text{rad}}=\frac{x.{\mu }_x+y.{\mu }_y}{r}$$(4) $${\text{v}}_{\text{tan}}=\frac{-y.{\mu }_x+x.{\mu }_y}{r},$$(5)

for stars with positions (x, y) relative to the cluster centre, proper motions (µ_x, µ_y), and the projected distance of each star to the centre of the cluster (r), as described by van Leeuwen et al. (2000).

We calculated the position angle of each star (θ_i) using the projected distance of the star from the centre of the cluster (x, y), with values varying between 0 < θ_i < 360. In general, the radial component (v_rad) as a function of the position angle is expected to be close to zero, as a non-zero value would indicate the cluster is experiencing radial fluctuations, implying the cluster is ‘breathing’. Negative v_rad values indicate stars approaching the cluster centre, while positive values indicate stars moving away from the centre. The tangential component (v_tan) can be negative, zero, or positive, with negative values indicating clockwise rotation of the stars and positive values indicating counter-clockwise rotation. In a rotating cluster, the LOS component (v_LOS) can display a sinusoidal rotation signal as a function of the position angle, with the amplitude of this sinusoidal signal dependent on both the strength of the cluster rotation and the orientation of the cluster with regard to the observer.

The rotation axis of stars within a rotating cluster can be expressed using the two angles shown in Figure 3 – the position angle of the rotation axis (0° ≤ θ₀ < 360°) and inclination angle (0° ≤ i ≤ 180°). The position angle describes the rotation axis of stars in the plane of the sky and was calculated for each cluster by plotting ∆v_LOS as a function of the position angle θ_i, fitting a sin curve to the rotation signal and determining the angle in which the sin function amplitude reached zero before increasing. The inclination angle i = arctan(v_tan/A_LOS), describes the inclination of the rotation axis in relation to the observer. An inclination of i = 0° corresponds to a face-on view of the cluster and therefore a strong v_tan signal, whereas i = 90° provides an edge-on view and a strong v_LOS signal. In order to clarify that these angles were recoverable from the combination of v_rad, v_tan and vLOS as functions of position angle θ_i, we created mock clusters using 100 000 stars with position and velocity components. Using different combinations of θ₀ and i input values and using only the projected positions and 3D velocities, we were able to reliably recover the input parameters within reasonable (<5°) uncertainties.

To recover the rotation axis angles of all stars in the clusters within the sample, we determined ∆v_rad, ∆v_tan and ∆v_LOS as functions of position angle θ_i, as shown in Figure 4. We used the Markov Chain Monte Carlo method (emcee in Python; Foreman-Mackey et al. 2013) to simultaneously fit the radial, tangential, and LOS velocity components of the stars in each cluster, using 32 walkers, 10 000 steps, and a burn-in of 2000 samples. For the radial and tangential components, we fit for the means (∆v_rad, ∆v_tan) and dispersions (${\sigma }_{{\text{v}}_{\text{rad}}},{\sigma }_{{\text{v}}_{\text{tan}}}$), while for the LOS component, we fit for the amplitude (∆A_LOS), phase shift (θ₀), and dispersion (${\sigma }_{{\text{v}}_{\text{L}}{}_{\text{OS}}}$). Figure 4 shows the result of this analysis for NGC 104: the top and middle panels show the distribution of stars (black) in each position angle bin for the radial and tangential velocities, respectively, with the averaged mean velocity ($\mathrm{\Delta }\overline{\text{v}}$) shown in the legends. The bottom panel shows the sinusoidal rotation signal recovered from the LOS velocities, with the position angle of the rotation axis (θ₀) shown as a white star with error bars.

The inclination angle i of each cluster was calculated using Equation (6) from the amplitude of the rotation signal in the LOS velocity and the mean tangential velocity. The uncertainty was calculated analytically using the associated errors of the LOS amplitude (∆A_LOS,err) and mean tangential velocity (∆v_tan,err), shown in Equation (7) $$i=\arctan \left(\frac{{\text{ΔA}}_{\text{LOS}}}{\text{Δ}{\overline{\text{v}}}_{\text{tan}}}\right)$$(6) $$\begin{array}{l}\text{Δ}i=\sqrt{{\left(\frac{\text{Δ}{\overline{\text{v}}}_{\text{tan}}}{{\text{ΔA}}_{\text{LOS}}^2+\text{Δ}{\overline{\text{v}}}_{\text{tan}}^2}\cdot {\text{ΔA}}_{\text{LOS},\text{err}}\right)}^2+}\\ \overline{{\left(\frac{-{\text{ΔA}}_{\text{LOS}}}{{\text{ΔA}}_{\text{LOS}}^2+\text{Δ}{\overline{\text{v}}}_{\text{tan\hspace{0.17em}}}^2}\cdot \text{Δ}{\overline{\text{v}}}_{\text{tan,err\hspace{0.17em}}}\right)}^2}.\end{array}$$(7)

We calculated the total rotational amplitude using the mean tangential velocity and amplitude of the LOS velocity signal ${\text{A}}_{\text{total\hspace{0.17em}}}=\sqrt{\mathrm{\Delta }{\text{A}}_{\text{LOS}}{}^2+{\overline{\text{v}}}_{\text{tan\hspace{0.17em}}}^2}$, as well as the total velocity dispersion $\sigma =\sqrt{{\sigma }_{\text{LOS}}{}^2+{\sigma }_{{\tan }^2}^2}$, in order to determine the total rotation over dispersion ratio (A_total/σ).

In order to calculate the photometric semi-major axis orientation of each cluster, we fit the spatial distribution of its members with a bivariate Gaussian distribution, using only stars brighter than the turn-off. To obtain a robust estimate of the position angle (PA) of the semi-major axis, we took the median and the standard deviation of 100 bootstrap experiments. In Table D.1, we report the estimated PA of the photometric semimajor axis and eccentricity of each cluster. These values were compared with those calculated in White & Shawl (1987), as well as Kamann et al. (2018) for clusters overlapping with our sample. We found consistency between the orientations from all three sources for clusters NGC 104, NGC 1851, and NGC 2808, but although we found good agreement in general with the values calculated by Kamann et al. (2018), we did not find our values aligned with those of White & Shawl (1987) and therefore we opted to only use our values for this analysis.

In the left panel of Figure 5, we show the position angle of the rotation axis of NGC 104 (θ₀ = 134 ± 4°) in cyan, with the corresponding photometric major axis orientation of the cluster shown in grey and the uncertainty represented by the width of the line. We found that for NGC 104, the position angle in the plane of the sky is perpendicular to the photometric major axis orientation, which is to be expected for a cluster flattened by initial rotation. The inclination angle (i = 33° ± 2°) for NGC 104 shown in magenta in the middle panel of Figure 5 indicates the orientation of the cluster is approximately mid-way between face-on and edge-on, which is also reflected in the strong signal of both the observed tangential and LOS velocities of Figure 4. The length of the inclination angle vector corresponds to the total rotational amplitude (A_total) of the cluster as displayed on the x- axis, while the width of the inclination angle vector indicates the uncertainty in the calculated value. Finally, the right panel of Figure 5 demonstrates the probability distribution of the total rotational amplitudes (A_total) calculated from the samples of the MCMC analysis for all stars.

Fig. 3 Schematic of the rotation axis angles for a globular cluster with an inclined plane. The light grey plane shows a face-on cluster with position vectors in the plane of the sky: RA, Dec, and the line-of- sight (LOS) vector towards or away from the observer. In dark grey is an inclined plane, with projected position vectors and a rotation axis described by two angles with respect to the face-on cluster: position angle of the rotation axis θ0 and inclination angle i.

Fig. 4 3D velocity components of NGC 104 with respect to the mean cluster velocity, as a function of the position angle θi . The combined sample of RGB, AGB, HB, and MSTO stars in each velocity catalogue are shown as black points with errorbars. Top panel: the radial component of the proper motion, with the mean value shown in green. Middle panel: the same as the top panel, but for the tangential component of the proper motion. Bottom panel: the LOS velocity as a function of θi, with the sinusoidal fit to the distribution shown in green and the position angle of the rotation axis θ0 shown as a white star with errorbars.

3.2 Rotation of the multiple stellar populations

The homogeneous classifications of the multiple stellar populations was completed for 30 Galactic globular clusters using a combination of ground-based (Stetson et al. 2019) and spacebased HST (Piotto et al. 2015; Nardiello et al. 2018) photometry, as described in Leitinger et al. (2023). The multiple stellar populations that were determined in Leitinger et al. (2023) were classified into: ‘P1’, ‘P2’, and in some clusters, ‘P3’. However, for the sake of the spatial distribution analysis described in Leitinger et al. (2023), as well as the kinematic analysis described in this paper, the ‘enriched stars’ (P2 and P3) were combined into one population (P2) for comparisons against P1, as our main focus is on the kinematic differences between stars that show ‘primordial’ chemical abundance patterns (P1), versus any stars that are enriched in comparison to P1 stars.

We added NGC 104 and NGC 6656 to the original sample of 28 GCs in Leitinger et al. (2023) by using the same classification method to determine the multiple stellar populations of the RGB stars, but removed NGC 5053, NGC 6981, and NGC 6934 due to the low number of LOS velocities for each cluster. We cross-matched the photometrically classified stars of each cluster with the cleaned LOS and proper motion catalogues described in Section 3.1, in order to assign a population to the stars in common with the velocity catalogues. We then performed the same analysis as described in Section 3.1 on the combined catalogues and determined the position angle of the rotation axis, inclination angle, and total rotational amplitude of each population in the 27 GCs.

The rotation signals as a function of the position angle for the multiple stellar populations in NGC 104 are shown in Figure 6, with P1 (P2) stars shown in the left (right) panels. The 3D velocity components exhibit a difference of 3.5σ significance in the mean velocities and rotation amplitudes between the P1 and P2 stars, with P2 stars rotating faster than P1 stars. Despite this, the left panel of Figure 7 shows similar position angles for P1 (blue) and P2 (red) stars, with uncertainties represented as solid lines. The middle panel shows only a slight difference in the inclination angles of the P1 and P2 stars, with the length of the vectors representing the total rotational amplitude and the width representing the uncertainties. The right panel shows the probability distribution of the calculated total rotational amplitudes A_total of each population, with solid lines indicating the median values and dashed lines indicating the ±1σ values.

NGC 104 has a dynamical age (age/relaxation time) of age/T_rh = 2.5 ± 0.1, meaning it has experienced just over two relaxation times, so while we classify it as one of the ‘dynamically younger’ clusters in our sample, the present day conditions may not be entirely representative of its initial conditions. With this in mind, we cannot confidently conclude the cluster was born this way – with P2 stars rotating faster than P1 stars and similar rotation axes for both populations. N -body simulations may help to explore how much the cluster has changed from its initial conditions and investigate scenarios that can produce the same kinematic results we have observed.

The rotation axes plots shown in Figures 5 and 7 were combined for the 27 clusters with adequate LOS velocities, to show the rotation axes and total rotational amplitudes of all RGB, AGB, HB, and MS-TO stars (green), as well as the RGB stars of the P1 (blue) and P2 (red) populations. As the data set for all stars tends to extend to a larger projected distance than the data set with population classifications, the position angles for each data set are shown at their corresponding average distances. These figures are publicly available online for the dynamically young clusters in our sample. Additionally, Table D.2 in Appendix D lists the calculated rotation parameters of the MPs, similar to Table A.1 for all stars.

Fig. 5 Calculated rotation parameters of NGC 104. Left panel: the position angle of the rotation axis (θ0) of NGC 104 is shown as a cyan circle with associated uncertainty, at the average projected distance of the velocity catalogues (r = 751 arcsec). In grey is the photometric major axis of the cluster derived in this work, with the uncertainty shown by the spread of the line. Middle panel: the inclination angle is shown in magenta, with associated uncertainty shown by the width of the line. An inclination angle of i = 0° indicates the cluster is face-on with respect to the observer, while i = 90° indicates an edge-on orientation. The length of the inclination angle line corresponds to the total rotational amplitude of the cluster (Atotal), shown on the x-axis. Right panel: the probability distribution of the total rotational amplitudes calculated for each MCMC sample.

Fig. 6 3D velocity components of NGC 104 as shown in Figure 4, but for the multiple stellar populations with P1 (blue) shown in the left panels and P2 (red) shown in the right panels.

3.3 Velocity dispersions

The velocity dispersions and the anisotropy for the full sample of 30 GCs were determined from the Gaia DR3 and HACKS proper motions. The number of stars with available proper motions are shown in Column 2 of Table A.1 for each cluster.

Since HACKS covers a square field of view, stars located close to the corners were removed by creating annuli from the centre of the cluster outwards, stopping at the radius in which there is 85% completeness. This method is described in greater detail in Section 2.1 of Leitinger et al. (2023), as it was also used for combining the HST and ground-based photometric catalogues. Stars in the Gaia DR3 catalogue located roughly within the HACKS field of view in the cluster centres suffered from stellar crowding, resulting in large proper motion uncertainties. For this reason, stars within the aforementioned radius limit were also removed, creating a smooth transition between the HACKS stars in the centre and Gaia stars in the outer regions. This was performed only for the velocity dispersion and anisotropy analysis, as the HACKS proper motions were not used to determine the rotation of the clusters.

The radial (v_rad) and tangential (v_tan) components of the proper motions were grouped into radial bins with equal numbers of stars, in which both the number of bins and number of stars in each bin depended on the sample available for each cluster.

A caveat on the calculation of anisotropy: there are two ways to calculate the velocity dispersions and resulting anisotropy of stars in a GC. In the first method, the mean velocity of each component (µ_rad, µ_tan) is subtracted from each radial bin, removing the bulk motion due to rotation and creating a reference frame in which only the intrinsic random motion between stars is observed. This reference frame is useful for observing the internal dynamical differences between stars, but does not necessarily reflect the general orbits of the stars – especially for fast rotating clusters. In the second method, the contribution of rotation is not subtracted, so that the bulk motion is also retained. This method can be useful for identifying the anisotropy of stars under the influence of rotation, that is then a combination of the bulk motion and intrinsic motion. In the case of a strongly rotating cluster, that greatly influences the orbits of the stars, this version of the anisotropy provides a more general view of these orbits. So for the purpose of identifying the internal dynamical processes, rotation-subtracted velocity dispersions are used, but for the purpose of observing the general orbits of stars in the context of energy, the rotation is included. The relative reference frame used to measure the HACKS proper motions does not include the bulk rotation of the cluster. Therefore, for the sake of consistency between the HACKS and Gaia proper motions, we used the method in which we remove the rotation to calculate the velocity dispersions and anisotropy, such that the bulk motion does not influence the anisotropy of the stars. However, we also present a version of the anisotropy results without the rotation removed and without including the HACKS proper motions (Figure C.1 in Appendix C).

We calculated the velocity dispersions of the radial and tangential components of the proper motions, using a modified version of the method outlined by Pryor & Meylan (1993). Briefly, we initially assume a normal distribution for the individual velocity components, described by the probability density function $$f\left(v_i\right)=\frac{1}{\sqrt{2\pi \left({\sigma }_c^2+{\sigma }_{e,i}^2\right)}}\exp \left(-\frac{{\left({\text{v}}_{\text{i}}-\overline{\text{v}}\right)}^2}{2\left({\sigma }_{\text{c}}^2+{\sigma }_{\text{e},\text{i}}^2\right)}\right),$$(8)

where $\overline{v}$ is the average cluster velocity, σ_c is the intrinsic cluster dispersion, v_i are the individual stellar velocities, and σ_e,i are the associated errors for the individual velocities. The likelihood function of all N stars was then calculated with the following set of equations: $${\displaystyle \sum _{i=1}^N\frac{v_i}{\left({\sigma }_c^2+{\sigma }_{e,i}^2\right)}}-\overline{v}{\displaystyle \sum _{i=1}^N\frac{1}{\left({\sigma }_c^2+{\sigma }_{e,i}^2\right)}}=0$$(9) $${\displaystyle \sum _{i=1}^N\frac{{\left(v_i-\overline{v}\right)}^2}{{\left({\sigma }_c^2+{\sigma }_{e,i}^2\right)}^2}}-{\displaystyle \sum _{i=1}^N\frac{1}{\left({\sigma }_c^2+{\sigma }_{e,i}^2\right)}}=0,$$(10)

that were numerically solved through iteration for both $\overline{v}$ and σ_c, where the initial guesses for these parameters assumed uniform uncertainties, leading to the simplified equations $$\overline{v}=\frac{1}{N}{\displaystyle \sum _{i=1}^N{v}_i}{\sigma }_c^2=\frac{1}{N}{\displaystyle \sum _{i=1}^N{\left(v_i-\overline{v}\right)}^2}-{\sigma }_{e,i}^2.$$

During this iteration process, stars deviating beyond a 2.5σ rejection tolerance from the median were systematically removed for each velocity dispersion bin until no further rejections were necessary, which will underestimate the true velocity dispersion by a few percent. The upper and lower bounds of the uncertainties were estimated by iteratively determining dispersion values that corresponded to a likelihood difference of 0.5, with respect to the maximum likelihood. Figure 8 shows the resulting radial and tangential velocity dispersion as a function of the projected distance for NGC 104, with dashed (solid) lines indicating the core (half-light) radius of the cluster.

Fig. 7 Position angles, inclination angles, and total rotational amplitude probability densities of NGC 104, as shown in Figure 5 for all stars, but now for the multiple stellar populations: P1 (blue) and P2 (red).

Fig. 8 Radial (indigo) and tangential (orange) velocity dispersion as a function of the projected distance for NGC 104. The dashed black line shows the core radius rc, while the solid black line shows the half-light radius rhi.

3.4 Anisotropy

Using the radial σ_rad and tangential σ_tan velocity dispersions, the anisotropy could be calculated using $\beta =1-\left(\frac{{\sigma }_{{\tan }^2}}{{\sigma }_{\text{rad}}{}^2}\right)$, providing limits of –∞ < β ≤ 1 for the possible anisotropy values. In order to remove the asymmetry in the possible values, we normalised this version of the anisotropy $$\tilde{\beta }=\frac{\beta }{(2-\beta )},$$(11)

in which the new limits are $-1\le \tilde{\beta }\le 1$, where $\tilde{\beta }=-1$ is fully tangentially anisotropic, $\tilde{\beta }=1$ is fully radially anisotropic, and $\tilde{\beta }=0$ is isotropic. We calculated this normalised anisotropy as a function of the projected distance from the centre of the cluster for the full sample of stars in each cluster, as demonstrated in Figure 9 for NGC 104, that shows isotropic behaviour for the inner ~100″ and radial anisotropy in the outer regions.

Since we used HACKS proper motions for the inner regions of each cluster, we compared and found consistency between our results and those of Libralato et al. (2022), when adopting the same definition of the anisotropy as Libralato et al. (2022): σ_tan/σ_rad. However, regardless of the initial conditions of cluster anisotropy, we would expect the inner regions of the cluster to become isotropic faster than the outer regions due to dynamical mixing. The outer regions of a cluster are expected to show the most significant deviations from isotropy. We especially found this to be the case for the dynamically youngest (age/T_rh < 4.5) globular clusters in the sample, that demonstrate a significant (>2σ) central concentration of either P1 or P2 stars (see Figure 15 of Leitinger et al. (2023)). This is discussed further in Section 4.3.

Fig. 9 Normalised anisotropy of NGC 104 as a function of the projected distance from the centre of the cluster in terms of arcseconds (bottom x-axis) and half-light radii (top x-axis). We removed the contribution due to rotation in this version of the anisotropy. The dashed black line indicates the core radius rc, while the solid black line indicates the halflight radius rhi.

3.5 Velocity dispersions and anisotropy of the multiple stellar populations

For the multiple stellar populations, we performed the same analysis as described in Sections 3.3 and 3.4, but for the cleaned proper motion data described in Section 2.2 matched with the combined HST and ground-based photometry outlined in Section 2.1. To ensure the two populations were comparable, we binned stars using multiple annuli from the centre of the cluster to the outermost regions. Since the anisotropy was investigated as a function of radius, it was important that the stars of each population occupy the same area in each bin. Otherwise, if we binned by equal number of stars in each population instead, we risked comparing P1 and P2 stars of different radii.

The median value of the projected distance for each bin varied slightly between each cluster, as shown in Figure 10 for NGC 104, but the bin edges were the same for each population. As NGC 104 had the largest sample of stars, the general radial anisotropy found for stars at ≳ 1 r_hl in Figure 9 was found to be driven by the P2 population more than P1. However, for many clusters, the anisotropy of the MPs was inconclusive due to the lower number of stars that matched with the photometric catalogue, leading to results such as NGC 3201 in Figure 11, for example. Although we found that NGC 3201 shows radial anisotropy in the outer regions, the errorbars of Figure 11 make it difficult to conclude whether one population is responsible for this, or both. Any significant anisotropy results for the MPs will be included and discussed in Section 5 for the relevant, dynamically young clusters.

Fig. 10 Normalised anisotropy of NGC 104 as a function of the projected distance from the cluster centre in terms of arcseconds (bottom x-axis) and half-light radii (top x-axis) for the P1 (blue) and P2 (red) populations. The dashed red line indicates isotropy. The dashed black line indicates the core radius rc, while the solid black line indicates the half-light radius rhi.

Fig. 11 Same as Figure 10 but for NGC 3201.

4 Results

This section focuses on the implications of the 3D kinematic analysis of the clusters in our sample, in terms of the rotation over dispersion ratio (A_total/σ), position angle of the rotation axis (θ₀), and inclination angle (i), for the full sample of stars in Section 4.1, and the subset that was classified into multiple stellar populations in Section 4.2.

Fig. 12 Rotation over dispersion ratio Atotal/σ as a function of: current mass (top left) (Baumgardt et al. 2023b), relaxation time (top right) (Baumgardt et al. 2023b), enriched star fraction (bottom left) (Leitinger et al. 2023) and concentration parameter (bottom right) (Harris 2010). Clusters are colour-coded by $A_4^{+}$ values (Leitinger et al. 2023) and each panel shows the Spearman and Pearson rank correlation coefficients for each parameter, with corresponding log p-values to show the significance of the correlations (log(p) ≤ −1.3 indicates a p-value ≤ 0.05, meaning the correlation is statistically significant and not a result of random chance).

4.1 3D rotation of the clusters

We calculated the 3D rotation over dispersion ratios (A_total/σ) for 27 Galactic GCs in our sample, as shown in Table A.1, in which we found significant (> 3σ) rotation for 21 of the clusters. We also investigated Atotal/σ as a function of every structural and orbital parameter available on the Galactic Globular Cluster Database (Baumgardt et al. 2023b), as well as parameters calculated in this work, our previous work (Leitinger et al. 2023), and the catalogue of parameters for Milky Way globular clusters (Harris 2010). Additionally, we divided the sample into in situ (seven GCs) versus accreted (23 GCs), using the classifications listed in Massari et al. (2019). We present the Pearson and Spearman rank correlation coefficients and p-values for the strongest correlations, as shown in Figure 12.

The strongest correlations with A_total/σ were found to be with the mass of the cluster (both current and initial), the fraction of enriched stars, the relaxation time, and the concentration of the cluster. We show these correlations in Figure 12, with in situ clusters marked as inverse triangles and accreted clusters marked as circles. We found a very strong correlation with the current mass of the cluster (log(Mass)), shown in the top left panel of Figure 12, that was also found when using the estimated initial mass of the cluster (log(M_i)). Each panel of Figure 12 is colour-coded by the $A_4^{+}$ parameters of Leitinger et al. (2023), indicating whether a cluster was found to be P1 centrally concentrated (blue), P2 centrally concentrated (red) or spatially mixed (black). Both the current and initial mass correlations produced Pearson and Spearman rank correlation coefficients of r_s = 0.85, r_p = 0.76 and r_s = 0.66, r_p = 0.58, respectively. Previously, Bianchini et al. (2018) analysed the rotation of 51 Galactic GCs using Gaia DR2 proper motions and compared V/σ against the ellipticity, relaxation time, metallicity, concentration, mass-to-light ratio, and total cluster mass. When isolating clusters with significant (>3σ) rotation, Bianchini et al. (2018) found moderate correlations with the mass, metallicity, and relaxation time, in that order of significance. However, they note that when considering all clusters, the correlation with metallicity is erased. Our results support the correlation with mass found by Bianchini et al. (2018), but when we also restrict our sample to only those clusters with significant (>3σ) rotation, it did not significantly increase our rank correlation coefficients for correlations with mass. We also did not find a strong correlation with metallicity in either case.

In the top right panel of Figure 12, we observe a strong correlation between rotation and relaxation time (log(T_rh[yr])), with the exception of three dynamically young clusters in particular – NGC 3201 (blue), NGC 5024 (red), and NGC 6101 (blue) – at log(T_rh) > 9.6 and A_total/σ ≈ 0.2. The only notable difference between these three clusters and the remaining clusters that follow the correlation, is that they exhibit some of the largest perigalactic distances of the sample (NGC 3201: R_peri = 8.34 kpc, NGC 5024: R_peri = 9.28 kpc, and NGC 6101: R_peri = 10.19 kpc), alongside NGC 4590 (R_peri = 8.88 kpc) that also exhibits low rotation with log(T_rh) = 9.48. This appears to simply be coincidental, as we found no correlations between rotation and eccentricity of the Galactic orbit, apo- or perigalac- tic distance. Regardless of these outliers, the correlation between rotation and relaxation time confirms the results of Kamann et al. (2018), Bianchini et al. (2018), and the N-body simulations of Tiongco et al. (2017), that conclude that if a cluster is born with an initial rotation, the ability of that cluster to retain its rotation signal is dependent on the relaxation time. As a cluster dynamically relaxes, angular momentum is redistributed towards the outer regions, effectively slowing the total rotation of the cluster.

In the lower right panel of Figure 12 a correlation between rotation and the concentration of a cluster can be seen, with concentration defined as the ratio between the tidal radius and core radius of a cluster (log(r_t/r_c)), taken from Harris (2010). This parameter is driven primarily by the density profile within a few half-light radii, indicating that core-collapsed clusters are rotating faster. In general, the P2 centrally concentrated clusters have higher concentration values (c > 1.7), while the P1 centrally concentrated clusters have lower concentrations (c < 1.3). The exception is NGC 6809, that has lost a significant fraction of its initial mass (M_cuɪrent/M_initial = 0.26; Leitinger et al. 2023). The two outlier clusters with low rotation (A_total/σ < 0.3) but high concentration (c = 2.5) are NGC 7099 and NGC 6752, both of which have undergone core collapse (Harris 2010), alongside NGC 7078 with c = 2.3, that is rotating significantly. It is possible that due to the redistribution of angular momentum during core collapse, the process of core collapse has affected each cluster differently, depending on its initial rotation or mass.

The strong correlation between rotation and the enriched star fraction (P₂/P_total) of the clusters in the lower left panel of Figure 12 was to be expected, as it traces the well established correlation between the mass of a cluster and the fraction of enriched stars it hosts (Milone et al. 2017; Leitinger et al. 2023). Similarly, more massive clusters tend to have longer relaxation times, while longer relaxation times allow clusters to retain their initial rotation more efficiently (Bianchini et al. 2018). We tested whether the strong correlation between rotation and relaxation time was also caused by a correlation with mass. To confirm this, we performed a residual analysis to remove the effect of mass from the rotation, finding that all three of the strong correlations we found (P₂/P_total , log(T_rh), and c) disappeared.

We also present the 3D rotation over dispersion ratios (A_total/σ) calculated in this work in Table A.2 alongside previous literature values for comparison. The fastest rotating clusters in our sample: NGC 104, NGC 5904, NGC 6656, and NGC 7078 are also shown to exhibit fast rotation in the vast majority of previous kinematic results in Table A.2. Overall there is only one rotating cluster in our sample: NGC 5986, that somewhat disagrees with previous literature, although we can only compare to Sollima et al. (2019) and Vasiliev (2019) who state the rotation is ambiguous, while Bianchini et al. (2018) found no rotation. In general, our results for the rotation of each cluster agree with at least one other previous literature result, where available (note that NGC 6101 was not included in any of the previous works listed).

Our sample of GCs has the largest overlap with the work of Martens et al. (2023), who performed a large-scale kinematic analysis of 25 Galactic GCs using LOS velocities derived from MUSE spectroscopy, 13 of which are present in our sample. Our works agree in finding unambiguous rotation for all 13 of these overlapping clusters, despite their analysis focusing only on the central regions of each cluster, without the use of proper motions. Additionally, our results agree for all except one cluster in common with the analysis of Petralia et al. (2024), who used APOGEE DR14 spectra to determine rotational signatures from LOS velocities. For the cluster in disagreement between our works: NGC 6218, we also see contention in whether it is rotating from the previous analyses in Table A.2.

Our rotation results also agree strongly with the results of Sollima et al. (2019), who analysed a large sample of 62 Galactic GCs using Gaia DR2 proper motions and the 2018 version of the radial velocity catalogue from Baumgardt & Hilker (2018). More specifically, we both found unambiguous rotation signals for seven clusters in common: NGC 104, NGC 2808, NGC 5904, NGC 6205, NGC 6656, NGC 7078, and NGC 7089. For the majority of the significantly rotating clusters in our sample such as NGC 5272, NGC 6218, NGC 6341, NGC 6752, and NGC 6809, Table 1 of Sollima et al. (2019) lists these clusters as having a P ≥ 99.9% chance of rotating, even if they classify the rotation as ‘uncertain’ (see Section 3.2 in Sollima et al. (2019) for more details). There is disagreement between our works for the clusters NGC 1851, NGC 5286, NGC 6121, and NGC 6254, however, we note that in each case it is because we found a rotation signal that was not apparent in the analysis of Sollima et al. (2019). We also note that both kinematic data sets in this work, Gaia DR3 and the LOS velocity catalogues first produced by Baumgardt & Hilker (2018), are updated versions of the data sets used by Sollima et al. (2019), and therefore allow for results with lower error bars, which may explain the discrepancies between our works.

We also found good agreement with Kamann et al. (2018), who completed an analysis on the rotation of 22 Galactic GCs using MUSE line-of-sight velocities, in which nine GCs show unambiguous rotation in both their work and ours, out of 12 GCs that overlap with our sample. For those other three overlapping clusters, we found rotation at a > 3σ significance, but Kamann et al. (2018) found rotation at only a 2σ significance, that may be due to MUSE targeting only the inner regions of each cluster, while we use velocity measurements in the outer regions as well. The position angles of the rotation axes calculated by Kamann et al. (2018) for the 12 overlapping clusters are in agreement with our results and show that the position angles are generally perpendicular to the photometric semi-major axis for most clusters. These results support the argument presented by Fabricius et al. (2014), in which an oblate, isotropic, rotating cluster is expected to have its rotation axis perpendicular to its flattening axis. Additionally, Kamann et al. (2018) discovered that five of the clusters in their sample showed evidence of a decoupled core, exhibited by a change in the position angle as a function of radius. This implies the centre of the cluster is rotating differently from the outer regions and is proposed to be explained by tidal interactions influencing the outer regions, while the central regions still retain traces of the rotational signatures created at birth. Of these five GCs with decoupled cores, three are part of our sample as well: NGC 5904, NGC 6254, and NGC 7078. All three of these clusters are fast rotators, but only NGC 5904 (age/T_rh = 3.6 ± 0.2) is considered dynamically young, albeit on the older end. It is expected that fast rotating GCs that have already gone through multiple relaxation times were rotating even faster at birth (Hénault-Brunet et al. 2015), which may also be the case for these clusters.

Fig. 13 Total rotational amplitudes divided by the velocity dispersion for the P1 (blue) and P2 (red) stars for each of the dynamically young clusters in our sample, as well as for all stars (green) as described in Section 3.1 and listed in Table A.1. Solid circles indicate that both populations contain enough stars with LOS velocities for a statistically robust result (NLOS,P1 + NLOS,P2 > 300), while open circles indicate the opposite. Cluster names are colour-coded based on whether the cluster has a significant P1 (blue) central concentration, P2 (red) central concentration or are spatially mixed (black). Clusters are ordered based on the Atotal/σ values derived using all available stars as listed in Table A.1, such that clusters on the left (right) have lower (higher) rotation.

4.2 3D rotation of the multiple stellar populations

We now focus on the dynamically youngest clusters in our sample (age/T_rh < 4.5), in terms of the 3D rotation over dispersion ratios (A_total/σ) of their multiple stellar populations, as shown in Figure 13 and displayed in Table D.2 in Appendix D. The x-axis of Figure 13 is ordered by the A_total/σ values using all stars, as described in Section 3.1, with NGC 288 showing the lowest degree of rotation (A_total/σ = 0.10 ± 0.04) and NGC 104 showing the highest (A_total/σ = 0.60 ± 0.02). The clusters in Figure 13 are colour-coded based on the values of $A_4^{+}$ calculated in Leitinger et al. (2023), indicating whether they exhibit a significant P1 central concentration (blue), P2 central concentration (red), or are spatially mixed (black). We found no consensus between whether one population rotates faster than the other, even if one population is more centrally concentrated. For example, in the P1 centrally concentrated clusters NGC 288, NGC 3201, and NGC 6101, we do not find that the P1 stars located primarily in the inner regions rotate faster than P2 stars. NGC 3201 and NGC 6101 could be argued as having P1 stars rotating faster than P2 stars, but NGC 288 shows the opposite (ignoring error bars). Similarly, the P2 centrally concentrated clusters NGC 104, NGC 5024, NGC 5272, and NGC 6809 show no evidence of the P2 stars consistently rotating faster. One may conclude the P2 stars of NGC 104 are rotating faster than the P1 stars, but in NGC 5024, NGC 5272, and 6809 (ignoring error bars), the opposite is observed. This general lack of evidence for one population preferentially rotating faster than the other is supported by Martens et al. (2023), but is at odds with the work of Dalessandro et al. (2024), who investigated the 3D kinematics of MPs in 16 Galactic GCs and observed that P2 stars preferentially rotate faster than P1 stars.

The largest difference between the A_total/σ values of the P1 and P2 stars in our sample of clusters is observed in NGC 104, with P2 rotating faster than P1 at a 3.5σ significance (see Table D.2), but with no differences in the position angles of the MPs. This result is supported by the work of Dalessandro et al. (2024), who found that P2 stars rotate faster than P1 stars in NGC 104, with no significant difference between the position angles. The only other cluster in which we found a difference of >2σ between the populations is NGC 7089, in which P1 is rotating faster than P2, but NGC 7089 was not part of the sample in Dalessandro et al. (2024). The results of Martens et al. (2023) for NGC 104 and NGC 7089 did not identify a difference in the rotation over dispersion, (υ/σ_HL, evaluated at the half-light radius) between the populations. However, Martens et al. (2023) did find differences between the rotation of the multiple stellar populations for NGC 2808, NGC 6093 (not in our sample) and NGC 7078, at around ~1–2σ significance. In NGC 2808 we found that P1 is rotating faster than P2 at a 1.8σ significance, that is supported by the result of Martens et al. (2023) with the same significance level. In NGC 7078 we found that P1 is rotating faster than P2 at a significance level of <1σ, while Martens et al. (2023) found the opposite: P2 stars rotating faster than P1 at a 2.2σ significance. As the work of Martens et al. (2023) focuses on the centre of the cluster with MUSE and NGC 7078 has been shown to exhibit a decoupled core (Kamann et al. 2018), itis expected that the inner and outer regions of the cluster are rotating differently, that may explain the disparity between our works for NGC 7078, while for NGC 2808 we produce compatible results. The calculated rotation values for the multiple populations of all clusters in our sample is shown in Table D.2 in Appendix D.

4.3 Anisotropy

In our anisotropy analysis, we isolated the dynamically youngest clusters in our sample (age/T_rh < 4.5) and normalised the projected radii of each cluster to its half-light radius. We then truncated every cluster to 5r_hl where possible and calculated the anisotropy for all cluster stars in our sample, shown together in Figure 14, with the contribution by rotation removed. All clusters have stars out to 5r_hl in Figure 14, except for NGC 5053. Here, we show that clusters that were found to have a P1 central concentration (A⁺ > 0.2) exhibit radial anisotropy in the outer regions, while clusters with a P2 central concentration (A⁺ < −0.2) tend to exhibit tangential anisotropy in the outer regions, with the exception of NGC 104 and NGC 5272. The dynamically old clusters in our sample (age/T_rh > 4.5) have mixed populations (A⁺ ~ 0) and are largely isotropic throughout, as shown in the bottom panel of Figure 14. But regardless of dynamical age, all clusters show approximately isotropic behaviour in their inner regions (<1r_hl in Figure 14), that are covered by the HACKS catalogue. In agreement with Watkins et al. (2015) and Libralato et al. (2022), isotropy is expected for the denser central regions of clusters, especially those that have undergone several relaxation times in their lifetime. However, this work has uncovered a previously undiscovered result: clusters with P1 central concentrations exhibit radial anisotropy in their outskirts, while those with P2 central concentrations show tangential anisotropy or isotropy in their outer regions.

An alternative version of Figure 14 is shown in Appendix C, in which rotation is not removed, therefore excluding stars with only HACKS proper motions due to the lack of rotation and instead only including the Gaia proper motions with rotation included. There are minimal differences in the anisotropy profiles for the majority of clusters whether the rotation is removed (Figure 14) or not (Figure C.1), except for NGC 104 that changed from radially anisotropic in the outer regions without rotation, to tangentially anisotropic with rotation. In Figure C.1, we find that all P2 centrally concentrated clusters except NGC 5272 are tangentially anisotropic in the outer regions, with NGC 5272 showing more isotropy in its outer regions than shown in Figure 14. All P1 centrally concentrated clusters still show radial anisotropy in the outer regions when rotation is included.

NGC 5024 is the dynamically youngest cluster in our sample to contain a P2 central concentration, followed by NGC 104 and then NGC 5272 in the upper panels of Figures 14 and C.1. NGC 6809 is another cluster with a P2 central concentration, but although it is still considered dynamically young, it has undergone a high degree of mass loss (Figure B.2) and is therefore not likely to be indicative of its initial conditions despite showing significant tangential anisotropy. Even so, it is difficult to make robust conclusions regarding a correlation between P2 central concentrations and tangential anisotropy, as two of the dynamically young clusters that contain spatially mixed populations according to Leitinger et al. (2023): NGC 4590 and NGC 5904, also exhibit clear tangential anisotropy in the outskirts, with NGC 4590 displaying the strongest tangential anisotropy in our sample. The remaining spatially mixed clusters: NGC 6205, NGC 6656, and NGC 7089 are approximately isotropic throughout. Clusters with the highest P1 central concentrations: NGC 288, NGC 3201, and NGC 6101 all exhibit radial anisotropy in the outskirts, which may be a result of dynamical evolution from an initially compact cluster (e.g. Tiongco et al. 2016), or a surviving signature from the initial conditions. We investigate these clusters in further detail in Section 5.2.

As the central regions of each cluster are approximately isotropic, the focus of the next section will be on the outer regions; specifically, we only use stars >2r_hl from the cluster centre and average the anisotropy values for these outer regions. In this way, we are focusing on the regions that have a greater chance of still being representative of the initial conditions at birth. We average the anisotropy bins in these outer regions – that corresponds to the bins with projected distance/half-light radius >2 in Figure C.1. For this section we only use the anisotropy profiles in which rotation is included, as the aim is to investigate trends with the general motion of the stars, as opposed to the intrinsic motion of the stars relative to one another.

We investigated the average rotation-included anisotropy in the outer regions covered only by Gaia ($\tilde{\beta }\left(r_{\text{hl}}\ge 2\right)$) as a function of every structural and orbital parameter available on the Galactic Globular Cluster Database (Baumgardt et al. 2023b), (Leitinger et al. 2023), (Harris 2010) and parameters calculated in this work. The clusters were classified as in situ (seven GCs) or accreted (23 GCs), using the classifications of Massari et al. (2019). We present the Pearson and Spearman rank correlation coefficients and p-values for each parameter that showed the strongest correlations, shown in Figure 15.

We found correlations with anisotropy between both the cumulative radial distribution parameter ($A_4^{+}$ from Leitinger et al. 2023) and the z component of the angular momentum (L_z). Accreted clusters are shown as circles, while in situ clusters are shown as inverted triangles. The clusters are colour-coded by their progenitor, as determined by Massari et al. (2019). Clusters with the highest P2 central concentrations ($A_4^{+}$ < −0.2) display either tangential anisotropy or isotropy. Of these clusters, NGC 5024 and NGC 5272 originate from the Helmi streams, NGC 104 formed in situ and NGC 6809 is part of an unclassified low- energy structure (see Massari et al. 2019 for details regarding the classifications). Clusters that formed in situ within the Milky Way are marked as red inverted triangles and are largely isotropic with spatially mixed populations. The exception to this is the in situ cluster NGC 104, that – when keeping the contribution of rotation in the anisotropy analysis – displays tangential anisotropy in its outer regions. The most P1 centrally concentrated clusters, NGC 3201 and NGC 6101 ($A_4^{+}$ > 0.3), are the most radially anisotropic clusters and are also the only two accreted clusters in our sample classified as originating from either Sequoia or Gaia-Enceladus, with Massari et al. (2019) stating that the progenitor is preferred to be Sequoia. It is possible that either 1) the specific star and cluster formation conditions within the progenitor system influenced both the P1 central concentration and radial anisotropy, 2) the accretion event of this progenitor influenced these parameters or 3) a combination of both.

The correlation between anisotropy in the outer regions and the z-component of the angular momentum further supports the idea that the physical conditions in the progenitor systems had an influence on the early stages of dynamical evolution in these clusters. However, there are three clusters in particular that seem to contribute the most to this correlation: NGC 4590 at L_z < −2kpc²/Myr, and NGC 3201 and NGC 6101 both at L_z > 2kpc²/Myr. In order to test the validity of this correlation, we experimented with removing certain clusters from the sample. The results of this experiment are shown in Table 1, that show the Pearson and Spearman rank correlation coefficients and associated log p-values for each test. We began by removing all in situ clusters (red inverted triangles in Figure 15) to focus only on accreted clusters, which strengthened the correlation. Removing the extremes in L_z: NGC 4590 or both NGC 3201 and NGC 6101 still produced strong correlations. Removing the in situ clusters, NGC 4590, NGC 3201, and NGC 6101 still produced a strong correlation for the Spearman coefficient, but the p-value for the Pearson coefficient suggested the result could be a product of random chance. Finally, keeping the in situ clusters, but removing the extreme L_z clusters: NGC 4590, NGC 3201, and NGC 6101, produced p-values for both coefficients that suggest the result is not statistically significant. Regardless, these extreme clusters are a part of our sample and the full sample of clusters does produce a strong correlation with L_z that, when coupled with the $A_4^{+}$ parameter correlation, suggests that the physical conditions of star cluster formation could be different in the former host galaxies, and thus influence the early stages of the dynamical evolution of their star clusters, leaving present-day clues in their outer regions. Future investigations would benefit from the inclusion of more clusters with extreme L_z values. For example, adding the anisotropy of the clusters NGC 7006, NGC 4499, NGC 1758, and Pal 13 – all of which exhibit L_z > 1 kpc²/Myr, along with Pal 10, Pal 12, NGC 5824, NGC 2419, Pal 4, Pal 3, and Pal 1 with L_z < −2 kpc²/Myr, would serve to fill in the gaps in Figure 15 and confirm whether the correlation is still significant.

Fig. 14 Normalised anisotropy of the dynamically youngest GCs in the sample (age/Trh < 4.5) (top) and the dynamically old clusters (age/Trh > 4.5) (bottom) as a function of the projected distance normalised by the half-light radius of each cluster. Each cluster is colour-coded by its A+ parameter, as calculated in Leitinger et al. (2023). Clusters in blue (red) are shown to have P1(P2) more centrally concentrated, while clusters in black have spatially mixed populations.

Fig. 15 Average normalised anisotropy in the outer regions (>2rhl) using only Gaia proper motions, without removing the contribution by rotation. Top panel: the correlation with the cumulative radial distribution parameter $A_4^{+}$ (Leitinger et al. 2023), and bottom panel: the correlation against the z-component of the angular momentum, where positive (negative) values indicate retrograde (prograde) orbits. Clusters are colour-coded by their progenitor structures, as classified by Massari et al. (2019). M-D: main disc of the Milky Way progenitor, L-E: unassociated low-energy group, H-E: high-energy unassociated group, H99: Helmi streams, G-E: Gaia-Enceladus, Seq: Sequoia. Each panel also shows the Spearman and Pearson rank correlation coefficients for each parameter, with corresponding log p-values to show the significance of the correlations (log(p) ≤ −1.3 indicates a p-value ≤0.05, meaning the correlation is statistically significant and not likely to be the result of random chance).

5 The dynamically youngest clusters

The focus in this section will be on the dynamically youngest clusters (age/T_rh < 4.5) in our sample, combining all available information from our spatial analysis in Leitinger et al. (2023) and kinematic analyses in this work, comparing with previous literature results where possible. The dynamically young clusters are separated into three sections: those that exhibit a P2 central concentration are discussed in Section 5.1, P1 centrally concentrated clusters are discussed in Section 5.2, and those with an approximately homogeneous mix of the populations are discussed in Section 5.3.

5.1 P2 centrally concentrated clusters

In this section we discuss the individual results of the clusters NGC 104, NGC 5024, NGC 5272, and NGC 6809, that contain a significant central concentration of the enriched stars. The majority of these clusters also exhibit tangential anisotropy in the outer regions, with the exception of the isotropic cluster NGC 5272 (See Figure 14) and NGC 104 when removing the effect of rotation (see the caveat in Section 3.3 and discussion in Appendix C). Additionally, each cluster was found to have higher concentration parameters as calculated by Harris (2010), with the exception of NGC 6809, that has undergone significant mass loss in comparison to the other clusters.

5.1.1 NGC 104 (47 Tucanae)

NGC 104 is one of the brightest, closest and most massive GCs in the Milky Way and has thus been the focus of many photometric and kinematic studies. In this work, we identified a P2 central concentration in NGC 104 that has also been well established in previous literature (Milone et al. 2012; Cordero et al. 2014; Jang et al. 2022), with a high enriched star fraction of P₂/P_total = 0.73 ± 0.02. It is also the fastest rotating cluster in our sample with A_total/σ = 0.60 ± 0.02, that is supported by many previous kinematic analyses (Lane et al. 2010; Bellazzini et al. 2012; Kimmig et al. 2015; Kamann et al. 2018; Gaia Collaboration 2018; Bianchini et al. 2018; Vasiliev 2019; Sollima et al. 2019; Martens et al. 2023). We found clear and high rotation of NGC 104, but with no differences between the P1 and P2 stars in terms of position angle of the rotation axis or inclination angle. However, we did find a difference between the rotation (A_total/σ shown in Table D.2 in Appendix D) of the MPs, in which P2 rotates faster than P1 at a 3.5σ significance, that supports the same result found by Dalessandro et al. (2024). However, this result is not supported by Martens et al. (2023) using MUSE LOS velocities, suggesting the contribution to rotation given by the outer regions and/or the tangential component of the proper motions contributes to this significant difference in rotation between the populations. A lack of significant rotational differences between the populations of NGC 104 was also reported by Milone et al. (2018) using Gaia DR2 proper motions, in which they observed P1 stars with higher tangential velocities than P2 stars, to a ~2σ level of significance – the opposite of our results. Similarly, Cordoni et al. (2020b) combined Gaia DR2 proper motions with VLT LOS velocities and found no significant differences between the rotation of the MPs. An analysis by Scalco et al. (2023) using Gaia DR3 proper motions and magnitudes converted into stellar masses through isochrone fitting, found a rotation-mass relation in which rotational velocity increases with stellar mass. They claim this is not the result of mass segregation for NGC 104, but is instead due to long term dynamical evolution causing more massive stars to rotate faster around the core of the cluster than low-mass stars. Scalco et al. (2023) did not assign population classifications to the stars in their sample, but since rotational velocity is a function of radius, the P2 central concentration in NGC 104 suggests that the P2 stars are on average located closer to the peak of the rotation curve than P1 stars.

In terms of anisotropy in NGC 104, we found an approximately isotropic central region, but two different results for the outer regions, depending on whether we remove the contribution by rotation or not (see the discussion in Appendix C). There is significant tangential rotation shown in the middle panel of Figure 4, that greatly affects the motion of the stars in the cluster, as the bulk motion of the stars is tangential motion. Removing this rotation when calculating the velocity dispersion reveals that the intrinsic random motion of the stars shows radial anisotropy in the outer regions. This result is consistent with the anisotropy analysis of Milone et al. (2018), Bianchini et al. (2018), Cordoni et al. (2020b), and Libralato et al. (2022), who also remove the contribution by rotation. Additionally, we found that the P2 stars are contributing to this radial anisotropy (Figure 10), with P1 stars exhibiting isotropy throughout the cluster. This approach allows us to observe the intrinsic motion of the stars in a nonrotating reference frame, such that the motion is only showing the intrinsic differences between stars, allowing assumptions to be made about the internal dynamical processes between stars. Because of this reference frame, Milone et al. (2018) concluded that the combination of rotation with the radial anisotropy of P2 stars is indicative of a scenario in which the P2 stars were born more centrally concentrated and are diffusing outwards due to the radial anisotropy, while P1 stars are isotropic throughout. However, we found that when considering the general orbits of stars in the context of the bulk motion due to rotation and the intrinsic motions (Appendix C), NGC 104 is tangentially anisotropic in the outer regions. Furthermore, we found this tangential anisotropy was driven by the P1 stars, whereas the faster rotating P2 stars are isotropic throughout. To understand the complex interplay between rotation and anisotropy observed in NGC 104, detailed dynamical modelling of the cluster will be needed, which is beyond the scope of this paper. Overall, NGC 104 shows high rotation with tangential anisotropy and no significant differences between the rotation axes of its MPs, but with P2 stars rotating faster than P1 stars.

5.1.2 NGC 5024 (M 53)

NGC 5024 is one of the dynamically youngest clusters in our sample, exhibiting a relatively low rotation over dispersion ratio with A_total/σ = 0.19 ± 0.05. Boberg et al. (2017) studied the rotation of NGC 5024 using spectra from the Hydra sample on the WIYN 3.5 m telescope for 245 stars, finding it has a peak rotation amplitude of 1.4 ± 0.1 km/s. These stars are a subset of the LOS catalogue we have analysed in this work, with our full sample of 344 stars exhibiting a lower rotation of ΔA = 0.7 ± 0.3 km/s for the LOS amplitude alone, but a total rotational amplitude of A_total = 1.4 ± 0.3 km/s. Additionally, Boberg et al. (2017) state that their discovery of a radial gradient of the position angle of the rotation axis suggests that the inner region of NGC 5024 is still reflective of the cluster formation and early dynamics, while the outer region is more affected by the Milky Way tidal field and interactions with the nearby cluster NGC 5053.

We found that the inner regions (<30 arcseconds) of NGC 5024 are slightly tangential, before shifting to radial anisotropy, with a tangential anisotropy again developing beyond ~2r_hl. Separating the sample into MPs shows that the two populations follow the same degree of anisotropy until the outer region >1 r_hl, where the P1 stars become radially anisotropic, while P2 stars show tangential anisotropy, as shown in Figure 16. It is unclear whether the radial anisotropy of P1 stars in the outer regions supports a scenario in which P1 stars were preferentially lost by the cluster. Given that NGC 5024 is considered to be associated with the cluster NGC 5053 due to an angular separation of only 0.96 degrees and a physical separation of ∼1 kpc (Jordi & Grebel 2010), the claim that interactions are altering the outer regions of NGC 5024 could imply that the outer regions are no longer representative of their initial conditions. Using SDSS photometry, Chun et al. (2010) claimed to identify a tidal bridge between the two clusters, but the analysis of Jordi & Grebel (2010), also using SDSS photometry, was unable to support this. If the outer regions have indeed been affected by interactions, there should also be evidence of stripped stars in the region surrounding the cluster; however, only seven stars were identified as potential members of NGC 5024 amongst the extratidal field star population by Hanke et al. (2020). A caveat included in the work of Hanke et al. (2020) is that the field between NGC 5024 and NGC 5053 was only covered by ∼ 10 fibres, making it unlikely to unveil evidence of a potential tidal bridge.

Nevertheless, the overall tangential anisotropy in the outer regions of NGC 5024 found in this work could be the result of: 1) the cluster being born with tangential anisotropy that has been preserved until the present day, or 2) the tidal field and interactions with NGC 5053 preferentially causing stars on radial orbits in the outer regions to be stripped, leaving only those stars that are on tangential orbits (possibly preferentially stripping P1 stars). Considering that NGC 5024 is a relatively compact cluster, (c = 1.72; Harris 2010), the stellar density of the cluster drops rapidly from the centre outwards, so it may be less likely for the outer regions to lose a significant enough fraction of stars to produce an observable tidal bridge within ~ 1.49 relaxation times. Future work could focus on identifying a tidal bridge between NGC 5024 and NGC 5053, or extratidal structures around NGC 5024 and, if identified, determining which population the stripped stars belong to.

Fig. 16 Normalised anisotropy of the P1 stars (blue) and P2 stars (red) in NGC 5024. The dashed vertical line indicates the core radius, while the solid vertical line indicates the half-light radius.

5.1.3 NGC 5272 (M 3)

We found significant rotation for NGC 5272, that has been established by many previous works to be rotating at a significance level of ≥2σ (Fabricius et al. 2014; Ferraro et al. 2018; Gaia Collaboration 2018; Bianchini et al. 2018; Sollima et al. 2019). We found that P1 stars appear to be rotating slightly faster (~1σ significance) than P2 stars, which is an unexpected result for a P2 centrally concentrated cluster, considering most formation theories predict the population in the centre to have a higher rotational amplitude. However, our results are at odds with Dalessandro et al. (2024), who found the opposite: P2 stars rotate slightly faster than P1 stars. There are many peculiarities surrounding this cluster, not only within our analysis, but also from previous literature. For example, NGC 5272 is the only dynamically young, P2 centrally concentrated cluster we have found to be approximately isotropic throughout, with slight radial anisotropy developing beyond ∼3r_hl. It has been proposed that NGC 5272 was either accreted from the Helmi stream (Massari et al. 2019) or the result of a merger between two GCs in a dwarf galaxy environment and then accreted (Lee 2021). For a cluster that is proposed to have had a tumultuous introduction to the Milky Way, there is no evidence of a large-scale tidal structure surrounding the cluster (Jordi & Grebel 2010; Hanke et al. 2020), which also agrees with the low mass loss fraction we reported in Leitinger et al. (2023) (M_c/M_i = 0.48). Another peculiarity is that the P1 stars are extended in the HST pseudocolour index used for chromosome maps: ΔC_F275W,F814W (Marino et al. 2019; Jang et al. 2022; Leitinger et al. 2023), which may be driven by small variations in Fe/H (Lardo et al. 2022). The combination of these elements creates great difficulty in tracing back the formation history of this cluster, despite its young dynamical age.

5.1.4 NGC 6809 (M 55)

In Leitinger et al. (2023) we discussed the structural properties of NGC 6809, in which the main focus was that for a dynamically young cluster (age/T_rh = 4.1 ± 0.2), it currently only has ~1/4 of its initial mass. Within this remaining mass, we identified a P2 central concentration in Leitinger et al. (2023), with P₂/P_total = 0.56 ± 0.04. In this paper we have found NGC 6809 to be tangentially anisotropic in the outer regions with a low overall rotation A_total/σ = 0.13 ± 0.02. The combination of these observables may suggest that the outer regions of NGC 6809 exhibit such strong tangential anisotropy because all stars on radial orbits were stripped away, leaving only those on tangential orbits.

5.2 P1 centrally concentrated clusters

In this section we focus on the P1 centrally concentrated clusters: NGC 288, NGC 3201, and NGC 6101. We note that in the case of NGC 3201, the $A_4^{+}$ parameter in Leitinger et al. (2023) suggests the cluster is P1 centrally concentrated, but the enriched star fraction P₂/P_total in both Leitinger et al. (2023) (Figure A3) and Cadelano et al. (2024) (Figure 5) suggests there are primarily P2 stars within ~100″ of the cluster centre, primarily P1 stars between ~100″ and ~300″ and finally, primarily P2 stars again beyond ~300″ . This result could suggest that NGC 3201 formed P2 centrally concentrated, with many of the inner P2 stars then distributed to the outer regions on radial orbits. But regardless of whether P2 stars are preferentially in the inner ~100″ , we still observe a greater fraction of P2 stars in the outer regions than P1 stars, which is why we still label the current state of this cluster as P1 centrally concentrated. In both NGC 288 and NGC 6101, the P₂/P_total fraction does not exhibit this trend (Leitinger et al. 2023), so we still consider them as fully P1 centrally concentrated clusters. The focus of this section is therefore to investigate any common features between these three clusters, that all primarily contain P2 stars in their outer regions.

NGC 288, NGC 3201, and NGC 6101 show some of the highest radial anisotropy in the outer regions when considering the intrinsic random motion of stars relative to each other (Figure 14), as well as the anisotropy including rotation (see Figure C.1). All three clusters have low rotation (A_total/σ < 0.25), and lower concentration values than the P2 centrally concentrated clusters according to Harris (2010) (with the exception of NGC 6809), which could be explained by Tiongco et al. (2016), who used N -body simulations to conclude that initially isotropic and tidally underfilling clusters (in terms of the halfmass radius to Jacobi radius ratio) can produce a build-up of radially anisotropic stars in the outer regions. In these cases, complete isotropy is not expected to be observable until they have undergone significant mass loss and dynamical evolution, that has not yet occurred for these dynamically young clusters. Interestingly, all three clusters were also found by Ibata et al. (2021) to have tidal tails, whereas no tidal tails were identified for any of the P2 centrally concentrated clusters.

5.2.1 NGC 288

NGC 288 has retained only a fraction (M_c/M_i ~ 0.24) of its initial mass, with an enriched star fraction of P₂/P_total = 0.37 ± 0.04 (Leitinger et al. 2023). From our kinematic analysis, we found that the P2 stars rotate slightly faster (<1σ) than the P1 stars, which is consistent with the result of Dalessandro et al. (2024). The outer regions of the cluster exhibit strong radial anisotropy, that can be expected for a cluster with such high mass loss. However, NGC 6809 has also experienced a similar degree of mass loss and instead exhibits tangential anisotropy in its outer regions. Significant tidal tails have been identified surrounding NGC 288, but not for NGC 6809 (Ibata et al. 2021), and the cluster has likely recently experienced a tidal shock from the disc and bulge (Leon et al. 2000), supporting a history in which NGC 288 has lost a large amount of its initial mass. As there were too few stars in the cross-match between our kinematic catalogues and our photometric catalogue, both P1 and P2 stars appeared approximately isotropic throughout, due to the size of the errorbars, so we cannot identify whether a specific population is responsible for the radial anisotropy.

5.2.2 NGC 3201

We found no significant differences between the rotation parameters of the MPs in NGC 3201, in terms of position angle, inclination angle, or total rotational amplitude, that is consistent with the results of both Martens et al. (2023) and Dalessandro et al. (2024). In our analysis of the anisotropy, we discovered NGC 3201 is radially anisotropic in the outskirts, both when the anisotropy is defined in terms of the intrinsic motion and, to a greater degree, when the anisotropy also includes the rotation (see Appendix C for details).

Bianchini et al. (2019) and Wan et al. (2021) investigated the peculiar kinematics in the outskirts of NGC 3201, that contains tidal tails and exhibits flattened velocity dispersions. The fact that P2 stars are primarily located in the outer regions of NGC 3201 could indicate that the peculiar kinematics found by Bianchini et al. (2019) and Wan et al. (2021) is driven by the P2 stars. By that same logic, the radial anisotropy we found in the outer regions of NGC 3201 has a higher chance of being driven by the P2 stars as well. Although our sample has too few stars to determine which population drives the radial anisotropy, the recent kinematics analysis of Cadelano et al. (2024) supports this idea of enriched stars driving the radial anisotropy. They classified the multiple stellar populations using the same method and photometric data as we did in Leitinger et al. (2023), obtaining the same results in terms of the cumulative radial distribution and enriched star fraction. They concluded that the bimodal distribution of the enriched star fraction is a result of an irregular 2D distribution of the populations, rather than a symmetrical distribution. Cadelano et al. (2024) then used Gaia DR3 proper motions to investigate the velocity dispersions and anisotropy of 297 P1 and 325 P2 stars. They discovered the P1 stars in their sample are isotropic throughout the cluster, while P2 stars are isotropic in the centre and become radially anisotropic beyond the half-mass radius. Such a result, coupled with the U-shaped P₂ /P_total distribution (Figure A3 of Leitinger et al. (2023) and Figure 5 of Cadelano et al. (2024)), shows evidence for a formation scenario in which the P2 stars were initially more concentrated before the expansion of the cluster ejected stars on radial orbits. If this is the case, one would expect to find a significant amount of P2 stars in the observed tidal tails of NGC 3201.

NGC 3201 has also been investigated in terms of hosting a black hole population, that could contribute to the peculiar structural and dynamical features observed. Giesers et al. (2019) concluded in their analysis of NGC 3201 that the cluster shows evidence of an extensive subpopulation of black holes, with Askar et al. (2018), Kremer et al. (2019), and Weatherford et al. (2020) estimating the number of black holes in NGC 3201 is between 100 < N_BH < 200, although Weatherford et al. (2020) also concludes that this number is significantly reduced to $N_{BH}={41}_{-34}^{+40}$ if the effects of mass segregation are accounted for, making it less likely that the black holes in NGC 3201 have significantly contributed to the peculiar current day observations.

5.2.3 NGC 6101

The rotation of NGC 6101 is mainly representative of the rotation of the inner regions covered by MUSE, as the majority of LOS velocities were taken from the ESO programme(s) 099.D- 0824(A), while only 26 stars were taken from the LOS catalogue compiled by Baumgardt & Hilker (2018). We found low rotation for NGC 6101, with no significant differences between the rotation amplitudes or position angles between the multiple populations. Our anisotropy analysis was performed using 2109 proper motions for all stars, showing the strongest radial anisotropy in the outskirts, of any cluster in our sample. Combining these proper motion measurements with the photometric sample of 481 stars in Leitinger et al. (2023) meant that the anisotropy of the MPs included only three bins, as shown in Figure 17. Changing the number of bins to two and four produced the same result – the P2 stars within the core radius were preferentially on radial orbits, while both populations were contributing to the radial anisotropy in the outer regions.

NGC 6101 has lost a very low amount of mass within its lifetime, but nevertheless, Ibata et al. (2021) discovered the presence of tidal tails for this cluster. Peuten et al. (2016) and Dalessandro et al. (2015) found no evidence of mass segregation within the cluster and, using N-body simulations, Peuten et al. (2016) concluded that a population of stellar mass black holes may be responsible for this. The analysis of Askar et al. (2018) performing simulations of GCs using the MOCCA code concluded that NGC 6101 is likely to be host to a large subsystem of 100 to 250 black holes, with Weatherford et al. (2020) supporting this idea with their confident estimation of 75–236 black holes. It may be possible that these stellar mass black holes have influenced the kinematics of the cluster, but additional observations and more data are required to further investigate this.

Fig. 17 Normalised anisotropy of the P1 stars (blue) and P2 stars (red) in NGC 6101. The dashed vertical line indicates the core radius, while the solid vertical line indicates the half-light radius.

5.3 Spatially mixed populations

This section focuses on the dynamically young clusters that have retained most of their initial conditions, but are nevertheless spatially mixed in terms of their populations. These clusters include NGC 4590, NGC 5053, NGC 5904, NGC 6205, and NGC 7089. These spatially mixed clusters show varying degrees of mass loss and rotation over dispersion ratios. We also observe variations between either isotropy or tangential anisotropy in their outer regions.

5.3.1 NGC 4590 (M 68)

We found NGC 4590 to have very low rotation, that is supported by the results of Lane et al. (2010), Kimmig et al. (2015), Bianchini et al. (2018), Vasiliev (2019), and Sollima et al. (2019), who also did not detect significant rotation. Interestingly, it exhibits the highest degree of tangential anisotropy in the outer regions, of any cluster in our sample. The 361 stars classified into MPs in Leitinger et al. (2023) were matched with the proper motion catalogues included in this work, and the resulting anisotropy of the MPs showed that the tangential anisotropy in the outer regions is driven mainly by the P2 stars. Changing the number of bins for this anisotropy analysis did not change the result of P2 stars showing tangential anisotropy, but the P1 stars varied between showing isotropy or slight radial anisotropy in the outer regions. This was also the case for the anisotropy of P1 and P2 stars in NGC 5024. Overall, the tangential anisotropy in the outer regions suggests that NGC 4590 is a dynamically stable cluster despite its young dynamical age.

From previous literature on NGC 4590, Baumgardt et al. (2019) found that NGC 4590 has large perigalactic (8.95 ± 0.06 kpc) and apogalactic (29.51 ± 0.42 kpc) distances, Massari et al. (2019) predicts that the Helmi streams is the progenitor of NGC 4590, and it has been reported to have very long tidal tails (Ibata et al. 2021). The large peri- and apogalactic distances suggest tidal stripping was unlikely to have removed a significant fraction of stars, but it is currently unclear how the accretion of the progenitor system of NGC 4590 into the Milky Way affected the kinematics of this cluster. NGC 5024, NGC 5272, and NGC 6981 were accreted from the same progenitor, but none of these clusters have tidal tails like NGC 4590. However, NGC 5024 and NGC 6981 exhibit tangential anisotropy in the outer regions similar to NGC 4590, while NGC 5272 exhibits slight radial anisotropy.

5.3.2 NGC 5053

We found NGC 5053 to have consistent radial anisotropy from the centre of the cluster to ~1r_hl, but there were too few stars in the photometry and LOS velocity catalogues for NGC 5053 in order for the anisotropy of MPs or the rotational analysis to be performed. Previous work has found NGC 5053 to be dynamically complicated as it contains significant tidal tails (Lauchner et al. 2006; Jordi & Grebel 2010) and it is postulated that NGC 5053 and NGC 5024 were accreted into the Milky Way from the same progenitor, possibly undergoing interactions with one another. Despite this common history, the two clusters show no similarities to one another. NGC 5053 shows slight radial anisotropy and has a concentration of c = 0.7, making it far less compact compared to the tangentially anisotropic NGC 5024 which has c = 1.7. Additionally, there is a P2 central concentration within NGC 5024, while NGC 5053 is spatially mixed, although both clusters have approximately similar enriched star fractions (Leitinger et al. 2023).

5.3.3 NGC 5904 (M 5)

NGC 5904 is a very interesting cluster in terms of kinematics and dynamics. We found NGC 5904 to be the second-fastest rotator in our sample (A_total/σ = 0.50 ± 0.03), consistent with previous kinematic analyses (Bellazzini et al. 2012; Fabricius et al. 2014; Kimmig et al. 2015; Kamann et al. 2018; Gaia Collaboration 2018; Bianchini et al. 2018; Vasiliev 2019; Sollima et al. 2019; Cordoni et al. 2020b; Martens et al. 2023). However, we found no differences in the rotation amplitudes, nor rotation axes, of the spatially mixed populations within NGC 5904, that is supported by the results of Martens et al. (2023), but partially at odds with the work of Dalessandro et al. (2024) who found P2 rotates faster than P1 (~2σ significance), with no differences between the rotation axes of the MPs. Cordoni et al. (2020b) also analysed the MPs of NGC 5904, finding the two populations exhibit the same rotational amplitude, but different rotation axes in terms of only the position angles. NGC 5904 has been found by Kamann et al. (2018) to exhibit a decoupled core, in which the position angle changes from the central region outwards, so this may account for the differences between works. The rotation-mass relation found for NGC 104 by Scalco et al. (2023) was also discovered for NGC 5904, meaning higher mass stars rotate faster around the cluster than low-mass stars (unrelated to mass segregation), but the spatially mixed populations within this cluster suggest that this may not necessarily be attributed to a specific population. There is also evidence of tidal tails surrounding NGC 5904 (Jordi & Grebel 2010; Grillmair 2019; Ibata et al. 2021), which is interesting as NGC 4590 is the only other cluster in our sample that exhibits significant tangential anisotropy in its outer regions and also has tidal tails. The progenitor system of NGC 5904 may be the Helmi streams, similar to NGC 4590, but the classification by Massari et al. (2019) is unclear and suggests NGC 5904 could belong to either the Helmi streams or Gaia-Enceladus.

5.3.4 NGC 6205 (M 13) and NGC 7089 (M 2)

We have combined the results of NGC 6205 and NGC 7089 into the same section due to the similarities between their masses, mass-loss ratios, dynamical ages, progenitor systems, concentration values, rotation, and anisotropy profiles. We found that NGC 6205 and NGC 7089 have the same anisotropy profiles, with both clusters exhibiting approximately isotropic behaviour in the inner regions, with slight radial anisotropy in the outer regions. We also found that the MPs of both clusters are approximately isotropic throughout. Both clusters share very similar concentration values (c = 1.5 for NGC 6205 and c = 1.6 for NGC 7089) and are classified as originating from Gaia-Enceladus (Massari et al. 2019). However, only NGC 7089 was observed to have tidal tails (Jordi & Grebel 2010; Hanke et al. 2020; Ibata et al. 2021), with one extratidal giant found to be a P2 star, that was argued to suggest there is a large fraction of P2 stars that have escaped the cluster (see Section 4.4 of Hanke et al. 2020 for caveats).

In terms of rotation, both NGC 6205 and NGC 7089 exhibit a high degree of rotation in our analysis. The same result for NGC 6205 has been found by Fabricius et al. (2014) and Sollima et al. (2019), as well as Bianchini et al. (2018) who discovered rotation in NGC 6205 to a 2σ level, while Vasiliev (2019) found no evidence of rotation in this cluster. For NGC 7089, our results agree with every literature result in common in Table A.2. Previously, Cordero et al. (2017) separated NGC 6205 into three populations, discovering the most enriched population rotates faster than the intermediate population, which they see as a signature of the initial formation. On the other hand, our results show the populations rotating to the same degree. We compared our sample with that of Cordero et al. (2017) and found agreement in the classification of the multiple stellar populations. We then found that we obtain the same result as Cordero et al. (2017) if we perform our rotation analysis on the same stars common between our data sets (113 stars), but that those differences disappear when using our full sample (855 stars), suggesting their analysis included too few stars. Using a larger sample of stars with 313 LOS velocities and 1201 proper motions, Dalessandro et al. (2024) found P2 stars rotating faster than P1 stars at ~1σ significance, consistent with our results for the MPs. In NGC 7089 we have found ~2σ differences in the rotation of the MPs, with P1 stars rotating faster than P2 stars despite the cluster being spatially mixed, but this cluster was not included in the analysis of Dalessandro et al. (2024) for comparison.

5.3.5 NGC 6656 (M 22)

We discovered significant rotation in NGC 6656, that is supported by previous results (Lane et al. 2010; Bellazzini et al. 2012; Kamann et al. 2018; Gaia Collaboration 2018; Bianchini et al. 2018; Vasiliev 2019; Sollima et al. 2019; Martens et al. 2023). However, we did not detect any significant differences in the rotation of the MPs, in terms of rotational amplitude or rotation axes. We cannot compare with the results of Martens et al. (2023), as they could not find differences due to a low number of stars in each population with associated MUSE velocities. NGC 6656 exhibits the most consistent isotropy amongst the dynamically young clusters in our sample, that is also supported by Cordoni et al. (2020a) using Gaia DR2 and HST data; however, they classify the populations as Fe-rich and Fe- poor, finding the populations share similar spatial distributions, with both populations isotropic. This classification of the populations is not directly comparable with our classifications, but considering we found isotropy for all stars, P1 stars, and P2 stars, while they found isotropy for Fe-rich and Fe-poor stars, it is reasonable to assume the populations are indeed isotropic in both analyses.

6 Conclusions

We have combined Hubble Space Telescope and Gaia DR3 proper motions together with a comprehensive set of line-of- sight velocities to determine the 3D rotational amplitudes, rotation axes, and anisotropy profiles of 30 Galactic globular clusters and their MPs. The focus of this paper was on the dynamically young clusters in our sample: NGC 104, NGC 288, NGC 3201, NGC 4590, NGC 5053, NGC 5024, NGC 5272, NGC 5904, NGC 6101, NGC 6205, NGC 6656, NGC 6809, and NGC 7089, as these clusters have the highest chance of retaining the initial conditions of formation.

We discovered significant (>3σ) rotation in 21 GCs, consistent with previous rotational analyses (e.g. Kamann et al. 2018; Bianchini et al. 2018; Sollima et al. 2019; Martens et al. 2023). We found no significant differences between the total rotational amplitudes of the MPs for the clusters in our sample, with the exception of NGC 104, in which P2 rotates faster than P1 at a ~3.5σ significance. We found no significant differences in the rotation axes of the MPs, in terms of position angle or inclination angle. For the clusters in our sample that demonstrated differences in the 3D rotational amplitudes of the MPs and exhibited a P1 or P2 central concentration, it was not necessarily the centrally concentrated population that displayed higher rotation. There has been increasing evidence in previous literature that there are very few clusters that exhibit significant differences in the rotation of the multiple populations (Milone et al. 2012; Libralato et al. 2019; Cordoni et al. 2020b; Szigeti et al. 2021; Martens et al. 2023). While for clusters in which differences are found, there has not been consistency between which population is rotating faster (i.e. P2 always rotating faster than P1) (Bellini et al. 2018; Kamann et al. 2020; Dalessandro et al. 2021; Szigeti et al. 2021; Martens et al. 2023). The exception to this is the recent work of Dalessandro et al. (2024), who concluded that P2 stars preferentially rotate faster than P1 stars for the 16 clusters in their sample. Although our results are mostly consistent with Dalessandro et al. (2024) for the nine clusters in common between our works, we cannot draw the same conclusion based on the clusters in our sample and the significance of the differences we find between the rotation of the MPs. A general lack of differences between the rotation of the MPs in our sample contradicts the expectations of some formation theories that place P2 stars as more centrally concentrated and with higher angular momentum compared to P1, therefore expecting that P2 stars will rotate faster (e.g. Hénault-Brunet et al. 2015).

We found that 3D rotation strongly correlates with the current and initial mass, relaxation time, enriched star fraction, and concentration of the clusters in our sample, but that these correlations are mainly attributed to correlations with mass. This implies that either the degree of initial rotation or the ability of a cluster to retain its initial rotation depends strongly on the mass.

We determined the anisotropy of each cluster, as well as the anisotropy of the MPs where possible, investigating correlations with the structural parameters, orbital parameters, and accretion history of the clusters from their progenitor systems. We discovered that the dynamically young clusters with the highest central concentration of primordial stars showed significant radial anisotropy in the outer regions (>2 half-light radii) of the cluster. Of these clusters, all were previously confirmed to exhibit tidal tails, while two were accreted from the same progenitor system. The dynamically young clusters with a central concentration of enriched stars displayed isotropy or tangential anisotropy in their outer regions, when also considering the rotation in the anisotropy analysis. These clusters generally displayed a higher mass and concentration (tidal radius/core radius) than the primordially concentrated clusters, with no confirmed tidal tails. We also found a correlation between the anisotropy in the outer regions of the clusters, with the z-component of the angular momentum of their orbits, suggesting the ancient host galaxy of a cluster not only ‘donors’ the orbit of the cluster around the Milky Way, but may also determine the kinematics of stars within the GC. This could also imply that the present day observed kinematics of a GC were determined by the physical conditions of star and star cluster formation within their progenitor systems, but further investigation into this correlation is necessary.

Data availability

The rotation plots similar to Figure 7 for all dynamically young clusters in our sample are publicly available on Zenodo.

Acknowledgements

This study was supported by the Klaus Tschira Foundation. Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme(s) 099.D- 0824(A). Based on observations made with the Gran Telescopio Canarias (GTC), installed in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, in the island of La Palma. This work is partly based on data obtained with MEGARA instrument, funded by European Regional Development Funds (ERDF), through Pro- grama Operativo Canarias FEDER 2014–2020. MG acknowledges support from the Ministry of Science and Innovation (PID2021-125485NB-C22, CEX2019- 000918-M funded by MCIN/AEI/10.13039/501100011033) and from AGAUR (SGR-2021-01069). EL acknowledges support from the ERC Consolidator Grant funding scheme (project ASTEROCHRONOMETRY, https://www.asterochronometry.eu, G.A. no. 772293).

Appendix A Tables of rotation parameters

Appendix B Updated cumulative radial distributions of the 30 Galactic GCs as function of the dynamical age and mass-loss ratio

Since Leitinger et al. (2023) was published, some structural parameters of the globular clusters in our sample have been updated on the Galactic Globular Cluster Database (Baumgardt et al. 2023b) from which they were taken. In particular, the relaxation times and current cluster masses have been updated, due to the inclusion of additional HST mass functions.

Following this update, the cumulative radial distribution plots presented in Leitinger et al. (2023) (Figure 15 as a function of the dynamical age and Figure 16 as a function of the mass-loss ratio) also required an update in order to reflect the clusters that are still considered ‘dynamically young’, matching the criteria of age/T_rh < 4.5.

The most notable changes in the dynamical ages occur for NGC 2808 (previously: age/T_rh = 3.6 ± 0.1 and currently: age/T_rh = 5.3 ± 0.2), NGC 3201 (previously: age/T_rh = 3.5 ± 0.2 and currently: age/T_rh = 2.5 ± 0.2), and NGC 7078 (previously: age/T_rh = 4.1 ± 0.1 and currently: age/T_rh = 8.0 ± 0.2). There are no substantial changes to the mass-loss ratios (M_current/M_initial) of each cluster.

Fig. B.1 Updated version of the normalised, cumulative radial distributions as a function of the dynamical age (Figure 15 from Leitinger et al. (2023)), including the two additional GCs, NGC 104 and NGC 6656. We have used updated dynamical ages from Baumgardt et al. (2023b), that change the sample of the ‘dynamically young’ clusters of our sample by removing NGC 2808 and NGC 7078, as the updated dynamical ages of these clusters now exceed age/Trh > 4.5. Please note that the colour-coding of the clusters in this figure do not correspond to the same colour-coding in Figure 15 of Leitinger et al. (2023).

Fig. B.2 Updated version of the normalised, cumulative radial distributions as a function of the mass-loss ratio (Figure 16 from Leitinger et al. (2023)), including the two additional GCs, NGC 104 and NGC 6656, as well as updated mass-loss ratios from Baumgardt et al. (2023b).

Appendix C The combined anisotropy

In Section 3.3, we calculated the velocity dispersion using the combination of HACKS and Gaia DR3 proper motions. As the HACKS proper motions were measured in a relative reference system that does not allow rotation to be observed, removing the rotational contribution from the Gaia proper motions allowed the two data sets to be compatible, as in Figure 14.

If the contribution by rotation is not first removed from the Gaia radial bins, there can be a discontinuity between the velocity dispersions for the HACKS radial range and the Gaia range. Due to this incompatibility, the HACKS proper motions cannot be used for the anisotropy calculation, which includes the contribution of rotation. We therefore included a second version of Figure 14 that shows the anisotropy using only the Gaia proper motions, in which rotation is not removed. Removing this bulk tangential motion as in Figure 14 reveals the intrinsic random motion of the stars, but this result complicates the physical meaning of the orbits of the stars we wished to use to investigate the correlations in Section 4.3.

In the Gaia only version of the anisotropy vs radius plot (Figure C.1) the most notable differences are in the clusters NGC 104 and NGC 3201. NGC 104 exhibits the fastest rotation of all clusters in our sample, with the tangential component of the rotation especially contributing to the orbits of the stars. Under the influence of rotation, the stars move preferentially on tangential orbits as expected. There are no other clusters in our sample that switch between radial and tangential anisotropy after the rotation is removed. The next most significant change is NGC 3201, that becomes less radially anisotropic after the rotation is removed. However, in general the anisotropy of the clusters do not significantly change when including or removing the rotation, nor do the other fast rotating clusters - NGC 5904 and NGC 6656 - show a significant difference between Figure 14 and Figure C.1.

Fig. C.1 Same as Figure 14, but using only Gaia proper motions, with no HACKS proper motions. We used the mean square velocity to calculate the anisotropy, such that the bulk motion and intrinsic random motion are both considered. For more details on this method, see the discussion above in Appendix C.

Appendix D Tables of the photometric semi-major axis orientation and rotation parameters of the multiple stellar populations

References

All Tables

All Figures

	Fig. 1 Sample of 30 Galactic globular clusters in galactocentric coordinates plotted over an artist’s conception image of the Milky Way (R. Hurt: NASA/JPL-Caltech/SSC) using mw-plot. The size of the points are proportional to the mass of the clusters.
In the text

	Fig. 2 Number of stars with measured proper motions (indigo) and LOS velocities (orange) as a function of the projected distance from the cluster centre for NGC 6809. The median values are shown by dashed lines.
In the text

Fig. 3 Schematic of the rotation axis angles for a globular cluster with an inclined plane. The light grey plane shows a face-on cluster with position vectors in the plane of the sky: RA, Dec, and the line-of- sight (LOS) vector towards or away from the observer. In dark grey is an inclined plane, with projected position vectors and a rotation axis described by two angles with respect to the face-on cluster: position angle of the rotation axis θ0 and inclination angle i.

In the text

Fig. 4 3D velocity components of NGC 104 with respect to the mean cluster velocity, as a function of the position angle θi . The combined sample of RGB, AGB, HB, and MSTO stars in each velocity catalogue are shown as black points with errorbars. Top panel: the radial component of the proper motion, with the mean value shown in green. Middle panel: the same as the top panel, but for the tangential component of the proper motion. Bottom panel: the LOS velocity as a function of θi, with the sinusoidal fit to the distribution shown in green and the position angle of the rotation axis θ0 shown as a white star with errorbars.

In the text

Fig. 5 Calculated rotation parameters of NGC 104. Left panel: the position angle of the rotation axis (θ0) of NGC 104 is shown as a cyan circle with associated uncertainty, at the average projected distance of the velocity catalogues (r = 751 arcsec). In grey is the photometric major axis of the cluster derived in this work, with the uncertainty shown by the spread of the line. Middle panel: the inclination angle is shown in magenta, with associated uncertainty shown by the width of the line. An inclination angle of i = 0° indicates the cluster is face-on with respect to the observer, while i = 90° indicates an edge-on orientation. The length of the inclination angle line corresponds to the total rotational amplitude of the cluster (Atotal), shown on the x-axis. Right panel: the probability distribution of the total rotational amplitudes calculated for each MCMC sample.

In the text

	Fig. 6 3D velocity components of NGC 104 as shown in Figure 4, but for the multiple stellar populations with P1 (blue) shown in the left panels and P2 (red) shown in the right panels.
In the text

	Fig. 7 Position angles, inclination angles, and total rotational amplitude probability densities of NGC 104, as shown in Figure 5 for all stars, but now for the multiple stellar populations: P1 (blue) and P2 (red).
In the text

	Fig. 8 Radial (indigo) and tangential (orange) velocity dispersion as a function of the projected distance for NGC 104. The dashed black line shows the core radius rc, while the solid black line shows the half-light radius rhi.
In the text

	Fig. 9 Normalised anisotropy of NGC 104 as a function of the projected distance from the centre of the cluster in terms of arcseconds (bottom x-axis) and half-light radii (top x-axis). We removed the contribution due to rotation in this version of the anisotropy. The dashed black line indicates the core radius rc, while the solid black line indicates the halflight radius rhi.
In the text

	Fig. 10 Normalised anisotropy of NGC 104 as a function of the projected distance from the cluster centre in terms of arcseconds (bottom x-axis) and half-light radii (top x-axis) for the P1 (blue) and P2 (red) populations. The dashed red line indicates isotropy. The dashed black line indicates the core radius rc, while the solid black line indicates the half-light radius rhi.
In the text

	Fig. 11 Same as Figure 10 but for NGC 3201.
In the text

Fig. 12 Rotation over dispersion ratio Atotal/σ as a function of: current mass (top left) (Baumgardt et al. 2023b), relaxation time (top right) (Baumgardt et al. 2023b), enriched star fraction (bottom left) (Leitinger et al. 2023) and concentration parameter (bottom right) (Harris 2010). Clusters are colour-coded by $A_4^{+}$ values (Leitinger et al. 2023) and each panel shows the Spearman and Pearson rank correlation coefficients for each parameter, with corresponding log p-values to show the significance of the correlations (log(p) ≤ −1.3 indicates a p-value ≤ 0.05, meaning the correlation is statistically significant and not a result of random chance).

In the text

Fig. 13 Total rotational amplitudes divided by the velocity dispersion for the P1 (blue) and P2 (red) stars for each of the dynamically young clusters in our sample, as well as for all stars (green) as described in Section 3.1 and listed in Table A.1. Solid circles indicate that both populations contain enough stars with LOS velocities for a statistically robust result (NLOS,P1 + NLOS,P2 > 300), while open circles indicate the opposite. Cluster names are colour-coded based on whether the cluster has a significant P1 (blue) central concentration, P2 (red) central concentration or are spatially mixed (black). Clusters are ordered based on the Atotal/σ values derived using all available stars as listed in Table A.1, such that clusters on the left (right) have lower (higher) rotation.

In the text

Fig. 14 Normalised anisotropy of the dynamically youngest GCs in the sample (age/Trh < 4.5) (top) and the dynamically old clusters (age/Trh > 4.5) (bottom) as a function of the projected distance normalised by the half-light radius of each cluster. Each cluster is colour-coded by its A+ parameter, as calculated in Leitinger et al. (2023). Clusters in blue (red) are shown to have P1(P2) more centrally concentrated, while clusters in black have spatially mixed populations.

In the text

Fig. 15 Average normalised anisotropy in the outer regions (>2rhl) using only Gaia proper motions, without removing the contribution by rotation. Top panel: the correlation with the cumulative radial distribution parameter $A_4^{+}$ (Leitinger et al. 2023), and bottom panel: the correlation against the z-component of the angular momentum, where positive (negative) values indicate retrograde (prograde) orbits. Clusters are colour-coded by their progenitor structures, as classified by Massari et al. (2019). M-D: main disc of the Milky Way progenitor, L-E: unassociated low-energy group, H-E: high-energy unassociated group, H99: Helmi streams, G-E: Gaia-Enceladus, Seq: Sequoia. Each panel also shows the Spearman and Pearson rank correlation coefficients for each parameter, with corresponding log p-values to show the significance of the correlations (log(p) ≤ −1.3 indicates a p-value ≤0.05, meaning the correlation is statistically significant and not likely to be the result of random chance).

In the text

	Fig. 16 Normalised anisotropy of the P1 stars (blue) and P2 stars (red) in NGC 5024. The dashed vertical line indicates the core radius, while the solid vertical line indicates the half-light radius.
In the text

	Fig. 17 Normalised anisotropy of the P1 stars (blue) and P2 stars (red) in NGC 6101. The dashed vertical line indicates the core radius, while the solid vertical line indicates the half-light radius.
In the text

Fig. B.1 Updated version of the normalised, cumulative radial distributions as a function of the dynamical age (Figure 15 from Leitinger et al. (2023)), including the two additional GCs, NGC 104 and NGC 6656. We have used updated dynamical ages from Baumgardt et al. (2023b), that change the sample of the ‘dynamically young’ clusters of our sample by removing NGC 2808 and NGC 7078, as the updated dynamical ages of these clusters now exceed age/Trh > 4.5. Please note that the colour-coding of the clusters in this figure do not correspond to the same colour-coding in Figure 15 of Leitinger et al. (2023).

In the text

	Fig. B.2 Updated version of the normalised, cumulative radial distributions as a function of the mass-loss ratio (Figure 16 from Leitinger et al. (2023)), including the two additional GCs, NGC 104 and NGC 6656, as well as updated mass-loss ratios from Baumgardt et al. (2023b).
In the text

	Fig. C.1 Same as Figure 14, but using only Gaia proper motions, with no HACKS proper motions. We used the mean square velocity to calculate the anisotropy, such that the bulk motion and intrinsic random motion are both considered. For more details on this method, see the discussion above in Appendix C.
In the text