We present the first 3D kinematic analysis of multiple stellar populations MPs in a representative sample of 16 Galactic globular clusters GCs. For each GC in the sample we studied the MP lineofsight planeofthesky and 3D rotation and velocity distribution anisotropy. The differences between firstpopulation FP and secondpopulation SP kinematic patterns were constrained by means of parameters specifically defined to provide a global measure of the relevant physical quantities and to enable a meaningful comparison among different clusters. Our analysis provides the first observational description of the MP kinematic properties and of the path they follow during their longterm dynamical evolution. In particular we find evidence of differences between the rotation of MPs along all velocity components with the SP preferentially rotating faster than the FP. The difference between the rotation strength of MPs is anticorrelated with the cluster dynamical age. We also observe that FPs are characterized by isotropic velocity distributions at any dynamical age probed by our sample. On the contrary the velocity distribution of SP stars is found to be radially anisotropic in dynamically young clusters and isotropic at later evolutionary stages. The comparison with a set of numerical simulations shows that these observational results are consistent with the longterm evolution of clusters forming with an initially more centrally concentrated and more rapidly rotating SP subsystem. We discuss the possible implications these findings have on our understanding of MP formation and early evolution.
['techniques: photometric', 'techniques: radial velocities', 'stars: abundances', 'Hertzsprung-Russell and C-M diagrams', 'stars: kinematics and dynamics', 'globular clusters: general', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Solar and Stellar Astrophysics']
A 3D view of multiple populations kinematics in Galactic globular clusters
{'A 3D view of multiple populations kinematics in Galactic globular clusters': 'E. Dalessandro 1 , M. Cadelano 1 , 2 , A. Della Croce 1 , 2 , F. I. Aros 3 , E. B. White 3 , E. Vesperini 3 , C. Fanelli 1 , F. R. 1 2 1 2 1 2 1  \nFerraro , , B. Lanzoni , , S. Leanza , , and L. Origlia  \n- 3 I-40129 Bologna, Italy\n- 1 INAF - Astrophysics and Space Science Observatory Bologna, Via Gobetti 93 / e-mail: emanuele.dalessandro@inaf.it\n- 2 Dipartimento di Fisica e Astronomia, Via Gobetti 93 / 2 I-40129 Bologna, Italy\n- 3 Department of Astronomy, Indiana University, Swain West, 727 E. 3rd Street, IN 47405 Bloomington, USA  \nReceived June 10, 2024; accepted September 02, 2024', 'ABSTRACT': 'We present the first 3D kinematic analysis of multiple stellar populations (MPs) in a representative sample of 16 Galactic globular clusters (GCs). For each GC in the sample we studied the MP line-of-sight, plane-of-the-sky and 3D rotation as well as the velocity distribution anisotropy. The di ff erences between first- (FP) and second-population (SP) kinematic patterns were constrained by means of parameters specifically defined to provide a global measure of the relevant physical quantities and to enable a meaningful comparison among di ff erent clusters. Our analysis provides the first observational description of the MP kinematic properties and of the path they follow during the long-term dynamical evolution. In particular, we find evidence of di ff erences between the rotation of MPs along all velocity components with the SP preferentially rotating faster than the FP. The di ff erence between the rotation strength of MPs is anti-correlated with the cluster dynamical age. We observe also that FPs are characterized by isotropic velocity distributions at any dynamical age probed by our sample. On the contrary, the velocity distribution of SP stars is found to be radially anisotropic in dynamically young clusters and isotropic at later evolutionary stages. The comparison with a set of numerical simulations shows that these observational results are consistent with the long-term evolution of clusters forming with an initially more centrally concentrated and more rapidly rotating SP subsystem. We discuss the possible implications these findings have on our understanding of MP formation and early evolution.  \nKey words. globular clusters: general - Stars: Hertzsprung-Russell and C-M diagrams - kinematics and dynamics - abundances techniques: photometry - spectroscopy - astrometry', '1. Introduction': "Globular clusters (GCs) exhibit intrinsic star-to-star variations in their light element content. In fact, while some GC stars have the same light-element abundances as field stars with the same metallicity (first population / generation-FP), others show enhanced N and Na along with depleted C and O abundances (second population / generation-SP). The manifestation of such light-element inhomogeneities is referred to as multiple populations (MPs - see Bastian & Lardo 2018; Gratton et al. 2019 for a review of the subject).  \nLight-element abundance variations can have an impact on both the stellar structures (as in the case of He for example) and atmospheres (as for Na, O, C and N) and they can therefore produce a broadening or splitting of di ff erent evolutionary sequences in color-magnitude-diagrams (CMDs) when appropriate filter combinations are used (Piotto et al. 2007; Sbordone et al. 2011; Dalessandro et al. 2011; Monelli et al. 2013; Piotto et al. 2015; Niederhofer et al. 2017; Cadelano et al. 2023).  \nIt is well established now that the MP phenomenon is (almost) ubiquitous among massive stellar clusters. In fact, it has been shown that nearly all massive ( > 10 4 M ⊙ ; e.g., Dalessandro et al. 2014; Piotto et al. 2015; Milone et al. 2017; Bragaglia et al. 2017) and relatively old ( > 1.5-2 Gyr; Martocchia et al. 2018a; Cadelano et al. 2022) GCs host MPs. In addition, MPs are observed in GC systems in any environment where they have been  \nproperly searched for. They are regularly found in the Magellanic Clouds stellar clusters (Mucciarelli et al. 2009; Dalessandro et al. 2016), in GCs in dwarf galaxies such as Fornax (Larsen et al. 2012, 2018) and Sagittarius (e.g., Sills et al. 2019), in the M31 GC system (Schiavon et al. 2013; Nardiello et al. 2018) and there are strong indications (though indirect) that they are ubiquitous in stellar clusters in massive elliptical galaxies (e.g., Chung et al. 2011).  \nMPs are believed to form during the very early epochs of GC formation and evolution ( ∼ 10 -100 Myr; see Martocchia et al. 2018b; Nardiello et al. 2015; Saracino et al. 2020 for direct observational constraints) and a number of scenarios have been proposed over the years to describe the sequence of physical events and mechanisms involved in their formation. We can schematically group them in two main categories. One includes multi-epoch formation models, which predict that MPs form during multiple (at least two) events of star formation and they typically invoke self-enrichment processes, in which the SP forms out of the ejecta of relatively massive FP stars (e.g. Decressin et al. 2007; D'Ercole et al. 2008; de Mink et al. 2009; D'Antona et al. 2016). The second category includes models where MPs form simultaneously and SPs are able to accrete gas eventually during their pre-main sequence phases (e.g., Bastian et al. 2013; Gieles et al. 2018). However, independently on the specific differences, all models proposed so far have their own caveats and face serious problems to reproduce some or all the available ob-  \nrvations. As a matter of fact, we still lack a self-consistent explanation of the physical processes at the basis of MP formation (e.g., see Bastian & Lardo 2018; Gratton et al. 2019).  \nUnderstanding the kinematical and structural properties of MPs can provide new insights into the early epochs of GC formation and evolution. In fact, most formation models suggest that MPs form with di ff erent structural and kinematic properties. Di ff erences between the FP and the SP kinematics can be either imprinted at the time of SP formation (see e.g., Bekki 2010; Lacchin et al. 2022) or emerge during a cluster's evolution as a consequence of the initial di ff erences between the FP and SP spatial distributions (see e.g., Tiongco et al. 2019; Vesperini et al. 2021; Sollima 2021). Although the primordial structural and kinematic di ff erences between FP and SP stars are expected to be gradually erased during GC long-term dynamical evolution (e.g., Vesperini et al. 2013; Hénault-Brunet et al. 2015; Tiongco et al. 2019; Vesperini et al. 2021; Sollima 2021), some clusters are expected to still retain some memory of these initial di ff erences. Indeed, Dalessandro et al. (2019) measured the di ff erence in the spatial distributions of FP and SP stars for a large sample of massive Galactic and Magellanic Clouds clusters homogeneously observed with HST. The authors found that the differences between the FP and SP spatial distributions generally follow the evolutionary sequence expected for the long-term dynamical evolution of clusters forming with an initially more centrally concentrated SP subsystem (see also Leitinger et al. 2023; Onorato et al. 2023; Cadelano et al. 2024).  \nSpatial distributions alone can provide only a partial picture of the dynamical properties of MPs and further key constraints on the possible formation and dynamical paths of MPs are expected to be hidden in their kinematic properties. Because of the technical limitations to derive kinematic information for large and significant samples of resolved stars in dense environments, most of the available information so far have been obtained by using HST proper motions (PMs) and ESO / VLT MUSE line-ofsight ( LOS ) velocities sampling relatively small portions of the cluster and focusing typically on the innermost regions. In a few particularly well studied systems, MPs have been found to show di ff erent degrees of orbital anisotropy (e.g., Richer et al. 2013; Bellini et al. 2015; Libralato et al. 2023) and possibly di ff erent rotation amplitudes (e.g., Cordero et al. 2017; Kamann et al. 2020; Cordoni et al. 2020; Dalessandro et al. 2021a; Martens et al. 2023). In other cases however, no significant di ff erences between the MP kinematic properties have been observed (see for example Milone et al. 2018; Cordoni et al. 2020; Libralato et al. 2019; Szigeti et al. 2021; Martens et al. 2023) and as a consequence, the emerging internal kinematic results describe a pretty heterogeneous picture.  \nTo move a leap forward in our understanding of MP kinematic properties and their possible implications on GC formation, it is fundamental to perform a systematic and homogeneous study of clusters sampling a wide range of dynamical ages and including the analysis of the clusters' outer regions, which are expected to retain some memory of the primordial structural and kinematic di ff erences for longer timescales. As a first step in this direction, in this paper, we perform for the first time a self-consistent study of the 3D kinematics of MPs in a representative sample of Galactic GCs for which it is possible to sample virtually their entire radial extension. This study has the additional advantage to overcome the typical limitations connected with projection e ff ects, which typically arise when LOS or planeof-the-sky velocities (i.e., PMs) are used independently possibly hampering the detection of the actual di ff erences between  \n<!-- image -->  \nFig. 1. On the left, the vector-point-diagram obtained by using Gaia DR3 PMs for the GC 47 Tuc is shown along with the distribution of stars along the µ ∗ α and µδ velocity components. The red circle represents the 2 σ selection described in Section 2.3. The panel on the right shows the (U, U-I) and (U, C U , B . I ) CMDs for 47 Tuc obtained by using groundbased photometric catalogs published by Stetson et al. (2019). Likely member stars based on the PM selection shown in the left panel are highlighted in black, while in grey likely field interlopers are shown.  \n<!-- image -->  \nthe MP kinematic properties (see e.g. the discussion concerning these issues in Tiongco et al. 2019).  \nThe paper is structured as follows. In Section 2 the adopted sample and the observational database are presented. In Section 3 we describe the kinematic analysis along the three velocity components, while in Section 4 we present the approach adopted for the study of the morphological properties of MPs. In Sections 5 and 6 we present the main observational results along with a detailed comparison with dynamical simulations, and in Section 7 we compare them with the literature. In Section 8 we summarize our findings and discuss their possible implications in the context of massive clusters formation and early evolution.", '2. Sample definition and observational data-sets': 'The analysis presented in this paper targets 16 Galactic GCs. In detail, the sample includes all clusters analyzed by Ferraro et al. (2018) and Lanzoni et al. (2018a,b) in the context of the ESO / VLT Multi Instrument Kinematic Survey of Galactic GCs (MIKiS) but NGC 362 as it lacks near-UV photometric data needed for the study of MPs (see Section 2.2 for details). We added to the target list NGC 104 (47 Tucanae) as it is a massive, relatively close and well studied GC, which can be useful for comparative analysis. We included also NGC 6362 for which we secured a large kinematic data-set in Dalessandro et al. (2018b, 2021a), NGC 6089 (M 80) and NGC 6205 (M 13) as they have been found to show interesting kinematic properties in previous analysis by Cordero et al. (2017) and Kamann et al. (2020). Table 1 summarizes some useful properties of the targets, such as distances, structural properties, age, relaxation times and kinematic sample sizes. The selected clusters are representative of the overall Galactic GC population as they properly encompass the cluster dynamically-sensitive parameter space, spanning a large range of central densities and concentrations, di ff erent stages of dynamical evolution and di ff erent environmental conditions. They are also more massive than M > 10 4 M ⊙ and relatively close to the Earth (within ∼ 16 kpc), thus providing data with good signal-to-noise ratios (S / Ns) for a large sample of stars.  \nTable 1. Properties of the 16 clusters fully analyzed in the present work.  \nNotes. Distances are from Baumgardt et al. (2019), structural parameters from Harris (1996) and ages come from the compilation by Dotter et al. (2010) with the exception of NGC 1904 for which we used the age derivation by Dalessandro et al. (2013). N LOS and N PM represent the number of LOS velocities and PMs used for the kinematic analysis.', '2.1. Kinematic database': "The analysis performed in this paper is based on two main kinematic data-sets securing LOS velocities (RVs) and PMs for hundreds (or thousands in a few cases - Table 1) of red giant branch stars (RGBs) in each GC. For 12 out of 16 GCs, most of the adopted LOS RVs were obtained by using ESO / VLT KMOS and FLAMESdataobtained as part of the MIKiS survey. We refer the reader to Ferraro et al. (2018) for details about the overall observational strategy and data-analysis. MIKiS LOS RVs were then complemented with those from the publicly available catalog by Baumgardt et al. (2019) to improve the sampling of the external regions of the target clusters. All LOS RVs used for 47 Tuc come from the Baumgardt et al. (2019) catalog. For NGC 6362, M 80 and M 13 LOS RVs were obtained by using ESO / VLTMUSEand FLAMES data (Cordero et al. 2017; Dalessandro et al. 2018b, 2021a; Kamann et al. 2020).  \nFor each of the investigated GCs, astrometric information, namely absolute PMs ( µ ∗ α , µδ ) and relative errors, were retrieved from the ESA / Gaia DR3(Gaia Collaboration et al. 2023) archive out to the clusters' tidal radius. Only stars with ruwe < 1 . 3 1 were then used for the kinematic analysis.", '2.2. Photometric dataset, membership selection and differential reddening correction': 'We used the wide-field catalogs published by Stetson et al. (2019) and including U, B, V, R and I bands to identify MPs in the target GCs (see Section 2.3) While these data are seeinglimited and can su ff er of incompleteness in the crowded central regions, they have similar spatial resolution as the kinematic LOS RVs dataset. In addition, they are typically not a ff ected by saturation problems and therefore they maximize the number of bright stars in common with the kinematic samples. The photometric catalogs were cross-matched with the kinematic ones based on their absolute coordinates ( α , δ ) and by using the cross-  \nFig. 2. LOS , RAD and TAN velocity distributions of likely member stars of the GC 47 Tuc as a function of the cluster-centric distance. All velocities are shown with respect to the cluster systemic velocity along the corresponding component.  \n<!-- image -->  \nrrelation tool CataXcorr 2 . For each cluster in the sample, the final catalog includes all stars in common between Gaia and the photometric catalogs. A fraction (typically larger than ∼ 60% along the RGB) of these stars has also LOS RVs and is therefore suited for a full 3D analysis (see Table 1).  \nTo separate likely cluster members from field interlopers, we selected stars whose PMs are within 2 σ from the cluster systemic velocity (adopted from Vasiliev & Baumgardt 2021), where σ is the standard deviation of the observed ( µ ∗ α , µδ ) distribution of RGB stars. We have verified that for the clusters in our sample, reasonable variations of the adopted cluster membership selection criteria do not have a significant impact on the main results of the kinematic analysis. As an example, Figure 1 shows  \n<!-- image -->  \nFig. 3. Left panel: (U, C U , B , I ) CMD of 47 Tuc likely member stars. RGB stars adopted for the kinematic analysis are highlighted in black. The blue and red curves are the fiducial lines adopted to verticalize the color distribution. Right-hand panels: the top panel displays the verticalized color distribution of RGB stars, while the bottom panel shows the corresponding histograms. The red and blue curves represent the two best-fit Gaussians for the FP and SP, respectively, while the solid black curve is their sum.  \n<!-- image -->  \nthe PM distribution of 47 Tuc stars along with the (U, U-I) and (U, CU , B , I - where CU , B , I = ( U -B ) -( B -I ); Monelli et al. 2013) CMDs for both selected cluster stars and field interlopers. Figure 2 shows instead the distributions of the velocities along the LOS , the PM radial ( RAD ) and PM tangential ( TAN ) components as a function of the cluster-centric distance for likely members RGB stars of the same cluster. It is worth mentioning here that the kinematic catalogs obtained from the MIKiS survey already rely on the cluster membership selection performed by Ferraro et al. (2018) and Lanzoni et al. (2018a,b), which is based on both the LOS RVand [Fe / H] distributions (we refer to those papers for further details).  \nAvailable magnitudes were then corrected for di ff erential reddening by using the approach described in Dalessandro et al. (2018c) (see also Cadelano et al. 2020). In short, di ff erential reddening was estimated by using likely member stars selected in a magnitude range typically going from the RGB-bump level down to about one magnitude below the cluster turn-o ff . By using these stars a mean ridge-line (MRL) was defined in the (B, B-I) CMD. Then for all stars within 3 σ (where σ is the color spread around the MRL) the geometric distance from the MRL ( ∆ X ) was computed. For each star in the catalog, di ff erential reddening was then obtained as the mean of the ∆ X values of the 30 nearest (in space) selected stars. ∆ X was then transformed into di ff erential reddening δ E ( B -V ) by using Eq. 1 from Dalessandro et al. (2018c), properly modified to account for the di ff erent extinction coe ffi cients for the adopted filters. Di ff erential reddening corrections turn out to be relatively small ( < 0 . 1mag) for most GCs in the sample, with the most critical cases being NGC 3201 and NGC 5927 for which we find δ E ( B -V ) values larger than ∼ 0 . 2 mag.', '2.3. MP classification': 'Starting from the sample of likely member stars and di ff erential reddening corrected magnitudes, we identified MPs along the RGB by using their distribution in the ( U , CU , B , I ) CMD (Figure 3). It has been shown that this color combination is very e ff ective to identify MPs along the RGB with di ff erent C and N (and possibly He) abundances (Sbordone et al. 2011; Monelli et al. 2013). RGB stars were verticalized in the ( U , CU , B , I ) CMD with respect to two fiducial lines on the blue and red edges of the RGB, calculated as the 5th and 95th percentile of the color distribution in di ff erent magnitude bins (Figure 3 - see, e.g., Dalessandro et al. 2018a,c; Onorato et al. 2023; Cadelano et al. 2023, for similar implementations of the same technique). In the resulting verticalized color distribution ( ∆ CU , B , I ; right panel in Figure 3), stars on the red (blue) side are expected to be N-poor (-rich), i.e., FP (SP), respectively. We run a two components Gaussian Mixture Modeling 3 (GMM) analysis to the resulting ∆ CU , B , I distribution thus assigning to each star a probability of belonging to the FP and SP sub-populations. Stars with a probability larger than 50% to belong to one of the two sub-populations were then flagged as FPs or SPs. Figure 3 shows the result of the MP identification and separation for the GC 47 Tuc. While this approach may introduce a few uncertainties and over-simplifications in the MP classification as we are not directly deriving light-element chemical abundances, it secures statistically large samples of stars with MP tagging that are hard to obtain by using only spectroscopic data.', '3. Kinematic analysis': "For each cluster in the sample we first analyzed the kinematic properties in terms of velocity dispersion, rotation and anisotropy profiles for the LOS and plane-of-the-sky components separately (Sections 3.1 and 3.2) and then, for the fraction of stars for which all velocity components are available, we performed a full 3D study (Section 3.3). We adopted as clusters' centers those reported by Goldsbury et al. (2010). All velocities were corrected for perspective e ff ects induced by the clusters' systemic motions by using the equations reported in van Leeuwen (2009) and following the approach already adopted in Dalessandro et al. (2021b); Della Croce et al. (2023).", '3.1. 1D velocity dispersion and rotation profiles': 'To characterize the kinematic properties of the clusters in the sample and of their sub-populations, we adopted the Bayesian approach described in Cordero et al. (2017) (see also Dalessandro et al. 2018b), which is based on the use of a discrete fitting technique to compare simple kinematic models (including a radial dependence of the rotational amplitude and velocity dispersion of the cluster) with individual radial velocities. We stress that this is a purely kinematic approach aimed at searching for relative di ff erences among di ff erent clusters and sub-populations and it is not aimed at providing a self-consistent dynamical description of each system.  \nThe likelihood function for the radial velocities of individual stars depends on our assumptions about the formal descriptions of the rotation and velocity dispersion radial variations. For the velocity dispersion profile we assumed the functional form of the Plummer model (Plummer 1911), which is simply defined by its central velocity dispersion σ 0 and its scale radius a :  \nFig. 4. Bottom panels : observed velocity dispersion profiles along the LOS , radial and tangential velocity components for MPs in 47 Tuc obtained by using a maximum-likelihood approach on binned data. Upper panels: Best-fit velocity dispersion profiles as obtained by using the bayesian analysis on discrete velocities described in Section 3.1.  \n<!-- image -->  \nσ 2 ( R ) = σ 2 0 p 1 + R 2 / a 2 , (1)  \nwhere R is the projected distance from the centre of the cluster. We adopted the same formal description for all velocity components. For the rotation curve, we assumed cylindrical rotation and adopted the functional form expected for stellar systems undergoing violent relaxation during their early phases of evolution (Lynden-Bell 1967):  \nV rotsin i ( X PA0) = 2 A rot R peak X PA0 1 + ( X PA0 / R peak) 2 (2)  \nµ TAN = 2 V peak R peak R 1 + ( R / R peak) 2 (3)  \nfor the LOS (Eq. 2) and TAN (Eq. 3) velocity components, respectively. In Eq. 2 Vrotsin i represents the projection of the rotational amplitude along the LOS velocity component at a projected distance X PA0 from the rotation axis. A rot is the peak rotational amplitude occurring at the projected distance R peak from  \nthe cluster center. We defined the rotation axis position angle ( PA ) as increasing anti-clockwise in the plane of the sky from north (PA = 0 · ) towards east (PA = 90 · ). Since the inclination of the rotation axis is unknown, Vrotsin i represents a lower limit to the actual rotational amplitude. As an extreme case, if the rotation axis is aligned with the line of sight, the rotation would be in the plane-of-the-sky. In Eq. 3, Vpeak represents the maximum (in an absolute sense) of the mean motion in the TAN component.  \nThe fit of the kinematic quantities was performed by using the emcee 4 (Foreman-Mackey et al. 2013) implementation of the Markov chain Monte Carlo (MCMC) sampler, which provides the posterior probability distribution function (PDFs) for σ 0, a , A rot and PA 0. For each quantity, the 50th-, 16th- and 84thpercentile of the PDF distributions were adopted as the best-fit values and relative errors, respectively. We assumed a Gaussian likelihood and flat priors on each of the investigated parameter within reasonably large range of values. It is important to note that in general, since the analysis is based on the conditional probability of a velocity measurement, given the position of a star, our fitting procedure is not biased by the spatial sampling of the stars in the di ff erent clusters and sub-samples. However, the kinematic properties are better constrained in regions that are  \nFig. 5. Same as 4, but now for the rotation profiles along the LOS (left panel) and tangential (right panel) components for MPs in 47 Tuc.  \n<!-- image -->  \nbetter sampled (i.e. larger number of stars with available kinematic information).  \nFor sanity check and comparison, we also derived the velocity dispersion and rotation profiles by splitting the surveyed areas in a set of concentric annuli, whose width was chosen as a compromise between a good radial sampling and a statistically significant number of stars. In this case the analysis was limited radially within a maximum distance from the center of the clusters to guarantee a symmetric coverage of the field of view. The adopted limiting distance varies from one cluster to the other depending on the photometric and kinematic dataset field of view limits (Section 2.3). While this approach requires the splitting of the sample in concentric radial bins, whose number and width is at least partially arbitrary and can potentially have an impact on the final results, it has the advantage of avoiding any assumption on the model description of the velocity dispersion and rotation profiles.  \nIn each radial bin the velocity dispersion was computed by following the maximum-likelihood approach described by Pryor & Meylan (1993). The method is based on the assumption that  \nthe probability of finding a star with a velocity of vi and error ϵ i at a projected distance from the cluster center Ri can be approximated as:  \np ( vi , ϵ i , Ri ) = 1 2 π q σ 2 + ϵ 2 i exp ( vi -v 0) 2 -2( σ 2 + ϵ 2 i ) (4)  \nwhere v 0 and σ are the systemic velocity and the intrinsic dispersion profile of the cluster along the three components (i.e., LOS , RAD and TAN ), respectively.  \nAs for the rotation along the LOS component, we used the method fully described in Bellazzini et al. (2012) and adopted in the previous papers of our group (e.g., Ferraro et al. 2018; Lanzoni et al. 2018a; Dalessandro et al. 2021a; Leanza et al. 2022). In brief, we considered a line passing through the cluster center with position angle varying from -90 · to 90 · by steps of 10 · . For each value of PA , such a line splits the observed sample in two. If the cluster is rotating along the line of sight, we expect to find a value of PA that maximizes the di ff erence between the median  \nRVs of the two sub-samples, since one component is mostly approaching and the other is receding with respect to the observer. Moving PA from this value has the e ff ect of gradually decreasing the di ff erence in the median RV. Hence, the appearance of a coherent sinusoidal behavior as a function of PA is a signature of rotation and its best-fit sine function provides an estimate of the rotation amplitude ( A rot) and the position angle of the cluster rotation axis (PA0). For the plane-of-the-sky rotation we used instead the variation of the mean values within each radial bin of the tangential velocity component with respect to the systematic motion.  \nExamples of the results obtained with both the Bayesian and maximum-likelihood analyses are shown in Figures 4 and 5 for the MPs of the GC 47 Tuc. For all clusters in the sample we find a good agreement between the discrete and binned analysis, however in the following we will adopt the best-fit results (and errors) obtained with the Bayesian approach. Table A.1 reports the best-fit values and relative errors for the most relevant quantities along both the LOS and plane-of-the-sky for both the FP and SP.', '3.2. Anisotropy profiles': "We adopted a modified version of the Osipkov-Merrit (Osipkov 1979; Merritt 1985) model to provide a formal description of the velocity anisotropy profile for each cluster and subpopulation (see Aros et al., in prep. for further details). This model is isotropic in the cluster centre and it becomes increasingly anisotropic for R > ra (where ra is the anisotropy radius). Our anisotropy description includes two additional parameters, which are the 'outer' anisotropy value ( β ∞ ) and the truncation radius ( rt ):  \nβ ( R ) = β ∞ R 2 ( R 2 + ra ) 2 1 -R rt ! (5)  \nβ ∞ defines how radial or tangentially biased the velocity anisotropy profile is for R >> ra , and rt is the radius at which the velocity anisotropy becomes isotropic again. In practice, we adopt rt as the Jacobi radius of the cluster.  \nFollowing the same approach adopted for the velocity dispersion and rotation, we performed a MCMC fit of the anisotropy profiles (see results in Table A.1) by using the individual stellar velocities along the radial and tangential components ( µ RAD and µ TAN, respectively) and by assuming a Gaussian likelihood in which the best-fit velocity dispersion along the radial direction σ RAD can be described through a Plummer profile (see Equation 1) and the velocity dispersion along the tangential direction is simply given by σ 2 TAN = [1 -β ( R )] σ 2 RAD . It is important to stress that this assumption is only adopted for the best-fit anisotropy parameters derivation, while we will generally consider as best-fit velocity dispersion profiles those derived through the independent fit of σ RAD and σ TAN discussed in Section 3.1. For the binned analysis, the anisotropy profiles were obtained directly from the RAD and TAN velocity dispersion values obtained in Section 3.1. Figure 6 shows an example of the best-fit anisotropy profile for the MPs of 47 Tuc.", '3.3. Full 3D analysis': 'To perform a full 3D analysis we used the kinematic sample of member stars having both LOS RVs and Gaia PMs after quality selection (see Section 2.1). We also limited the analysis to  \nthe same radial extension adopted for the binned maximumlikelihood analysis (Section 3.1).  \nWe followed the approach described in Sollima et al. (2019), which has the advantage of constraining a cluster full rotation pattern by estimating the inclination angle of the rotation axis ( i ) with respect to the line of sight, the position angle of the rotation axis ( θ 0) and the rotation velocity amplitude ( A ), by means of a model-independent analysis. We refer the reader to that paper for a full description of the method. Here we only report the main ingredients and assumptions.  \nIn a real cluster the angular velocity is expected to be a function of the distance from the rotation axis (see for example Eqs. 1 and 2). To account for such a dependence in a rigorous way, a rotating model should be fitted to the data. However, to perform a model-independent analysis we considered an average projected rotation velocity with amplitude A 3D, which has been assumed to be independent on the distance from the cluster center. While of course, this represents a crude approximation of the rotation patterns expected in GCs and provides a rough average of the actual rotation amplitude, it is important to stress that it does not introduce any bias in the estimate of θ 0 and i .  \nA 3D, i and θ 0 were derived by solving the equations describing the rotation projection along the LOS ( V LOS) and those perpendicular ( V ⊥ ) and parallel ( V ∥ ) to the rotation axes (see Eq. 2 in Sollima et al. 2019). While the velocity component perpendicular to the rotation axis has a dependence on stellar positions within the cluster along the line of sight, we neglected it in the present analysis. In fact, we note that the dependence on the stellar distance along the line of sight does not a ff ect the mean trend of the perpendicular velocity component, but it can only introduce an additional spread on its distribution. We assumed that i varies in the range 0 · < i < 90 · with respect to the line of sight and the position angle in the range 0 · < θ 0 < 360 · . θ 0 grows anticlockwise from North to East and A 3D is positive for clockwise rotation in the plane of the sky. Following the approach already adopted for the 1D analysis, we derived the best-fit rotation amplitudes, position and inclination angles and relative errors by maximizing the likelihood function reported in Eq. 3 of Sollima et al. (2019) by using the MCMC algorithm emcee . Best-fit results are reported in Table A.1. Figure 7 shows the result of the best-fit analysis along the three velocity components for the FP and SP sub-populations of 47 Tuc.', '4. Morphological analysis: MP ellipticity': "We inferred the ellipticity of the MPs by constructing the socalled shape tensor (Zemp et al. 2011), which is defined by the following equation  \nS ij ≡ P k m k( R k)i( R k)j P k m k , (6)  \nwhere R k and m k are the projected distance from the cluster center and the mass of the k -th star, within the i -th and j -th element of the shape tensor grid. For the purpose of this study, we assumed that all stars have the same mass. This assumption is well justified as the adopted sample consists of RGB stars (the same used for the kinematic analysis), whose mass variations for a given cluster is expected to be ∼ 0 . 02 M ⊙ , depending on the metallicity.  \nThe shape tensor was computed starting from a spherical grid, whose nodes are the stellar radial distribution's 10th, 50th, and 90th percentiles. We adopted an iterative procedure for each bin during which the shape tensor is initially constructed from  \nFig. 6. Same as 4, but now for the anisotropy profile (see Section 3.2).  \n<!-- image -->  \nspherical distances. After that, being ( w 0; w 1) the eigenvalues (with w 0 > w 1) and ( v 0; v 1) the respective eigenvectors of the shape tensor, it follows (Zemp et al. 2011) that  \nϵ = 1 -w 1 w 0 and PA = arctan v 0 , y v 0 , x . (7)  \nThe particle coordinates were then rotated by the angle -PA, and distances to the center were defined by means of the circularized distance  \nR ell ≡ s x ' 2 + y ' 2 ϵ 2 , (8)  \nwhere ( x ' ; y ' ) are the rotated locally-Cartesian coordinates of the stars.  \nArticle number, page 8 of 26  \nFinally, stars were binned in the new coordinate system according to the initial grid and the shape tensor has been computed using R ell instead of R . Such procedure was then iterated until a relative precision of 5% on the axis ratio was reached. Best-fit ellipticity values and relative errors have been obtained by means of a bootstrapping analysis. In detail, the shape tensor was computed 1000 times for each cluster by randomly selecting at each time a sub-sample including only 90% of stars. The values corresponding to the 50th percentile of the distribution of all the ϵ values has been then adopted as best-fit, while the errors correspond to the 16th and 84th percentiles. As an example, Figure 8 shows the result of the analysis for the MPs in 47 Tuc, while Table A.1 reports the best-fit values for each GC and subpopulation.  \nFig. 7. Distribution of the three velocity components as a function of the position angle for MPs in 47 Tuc. Left and right panels refer to the FP and SP sub-population, respectively. The solid lines show the best-fitting trend in all panels.  \n<!-- image -->  \nFig. 8. 2D density maps of the stars selected in the GC 47 Tuc for the kinematic analysis. FPs are shown on the left panel, while SPs in the right panel. Overplotted to the density distributions are the best-fit ellipses derived as described in Section 4.  \n<!-- image -->  \n-  \n-", '5. Results: rotation': 'To obtain quantitative and homogeneous estimates of the possible kinematic di ff erences among MPs, to follow their evolution and eventually to compare the results obtained for all GCs in the sample with theoretical models and dynamical simulations, we introduced a few simple parameters described in detail in the following. These parameters are meant to incorporate in a meaningful way all the main relevant physical quantities at play in a single value. The general approach of our analysis is not to focus  \non the detection of specific and particularly significant kinematic di ff erences of specific targets, but rather to compare the general kinematic behaviors described by MPs in all targets in the sample in the most e ff ective way.  \nA description about the approach adopted to quantify the possible impact of the intrinsically limited statistical kinematic samples and of their incompleteness on the final results is discussed in Appendix B. Here we briefly stress that the main effects are not on the derived best-fit values, but rather on their uncertainties. While the main focus of the following sections is  \non the MP kinematics, we analyzed also the entire sample of stars with kinematic information (hereafter labeled as TOT) for comparison with previous works and we present the main results in Appendix C.', '5.1. Observations': "To measure the rotation di ff erences between the SP and FP subpopulations for all clusters in the sample for both the LOS and TAN velocity components, we introduced a parameter, hereafter referred as α , defined as the area subtended by the ratio between the best-fit rotation velocity profile and the best-fit velocity dispersion profile for each sub-population in a cluster (Section 3.1) after rescaling the cluster-centric distance to the value of the peak ( R α ) of such a distribution:  \nα X = Z 1 0 V rot( R l) /σ ( R l) dR l (9)  \nwhere V rot( R l) and σ ( R l) can either refer to the LOS or the TAN velocity components, R l is the cluster-centric distance normalized to the R α and the index X refers to the di ff erent sub-populations (i.e., FP or SP).  \nThis parameter has the advantage of providing a robust measure of the relative strength of the rotation signal over the disordered motion at any radial range without making any assumption about the underlying star or mass distribution. By construction α depends on the considered cluster-centric distance and therefore a meaningful cluster-to-cluster comparison requires that the parameter is measured over equivalent radial portions in every system. As shown in a number of numerical studies (see, e.g., Hénault-Brunet et al. 2015; Tiongco et al. 2019), dynamical evolution is expected to smooth out primordial kinematic and structural di ff erences in the innermost regions first and then in the cluster's outskirts. Therefore, capturing rotation di ff erences between MPs require a compromise between probing a fairly wide radial coverage, thus to trace regions where kinematic di ff erences should be present for a longer time, and sampling distances from the cluster center where rotation is more prominent. With this in mind, we decided to measure α within R α . We verified also that the adoption of di ff erent radial selections does not have a significant impact on the overall relative distribution of α values. Errors on α were obtained by propagating the posterior probability distributions obtained from the MCMC analysis for the best-fit rotation and velocity dispersion profiles derivation (see Section 3.1). Di ff erences between SP and FP kinematic patterns are constrained simply by ( α SP -α FP ). With such a definition, a more rapidly rotating SP yields positive values of ( α SP -α FP ).  \nFollowing similar lines, we defined a parameter to describe the 3D rotation:  \nω X 3 D = ( A 3D /σ 3D 0 ) / ( R m / R h) (10)  \nwhere σ 3D 0 represents the 3D central velocity dispersion and it is defined as the quadratic average of the σ 0 values obtained for the three velocity components (i.e., LOS , TAN and RAD - Section 3.1); R m is the average cluster-centric distance of stars for which we have tri-dimensional velocity measures, and R h is the system half-light radius (from Harris 1996). Here R h is adopted as a meaningful radial normalization factor to secure a direct comparison among di ff erent GCs attaining significantly di ff erent projected radial extension. In the assumption of a pure solid-body rotation, ω 3D would represent the best-fit angular rotation. As  \nfor the 1D analysis, the introduction of ω 3D is primarily meant to provide a direct and reliable characterization of the 3D rotation based only on quantities that are directly derived from the observations. Di ff erences in the 3D rotation of SP and FP are given by ( ω SP 3D -ω FP 3D ), which yields positive values for a more rapidly rotating SP.  \nSeveral works have shown that the rotation strength observed in GCs is primarily shaped by their dynamical age, with dynamically young systems typically showing the larger degree of rotation (e.g. Fabricius et al. 2014; Bianchini et al. 2018; Kamann et al. 2018; Sollima et al. 2019). We used N h = t age / t rh as a proxy of the clusters' dynamical age. We adopted the the half-mass relaxation time ( t rh) values reported by Harris (1996) and ages derived by Dotter et al. (2010) for all clusters but NGC 1904, for which we used the age inferred by Dalessandro et al. (2013).  \nIn Figure C.1 in the Appendix we show the distribution of the α values (along both the LOS and TAN velocity components) and of ω 3D as a function of N h for the TOT population. In general, we find a very good agreement with previous analysis (e.g., Bianchini et al. 2018; Kamann et al. 2018; Baumgardt et al. 2019; Sollima et al. 2019) in terms of correlation between cluster rotation strength and dynamical age, thus further strengthening the idea that the present-day cluster rotation is the relic of that imprinted at the epoch of cluster formation, and that it has been then progressively dissipated via two-body relaxation (Einsel & Spurzem 1999; Hong et al. 2013; Tiongco et al. 2017; Livernois et al. 2022; Kamlah et al. 2022). Interestingly, such an agreement provides also an independent assessment about the reliability of the adopted kinematic parameters.  \nIn Figures 9, 10 and 11 we show the distributions of α (for both the LOS and PM components) and of ω 3D as a function of N h for both SP and FP stars (bottom and middle panels). Interestingly, a number of common patterns can be highlighted in the three Figures. First, we note that in a large fraction of clusters (up to ∼ 50%) both FP and SP show evidence of non-negligible rotation. Secondly, both sub-populations show evidence of anticorrelation with N h in all three analyzed velocity components, with dynamically young clusters being characterized by a larger rotation strength. This behavior turns out to be more prominent when the PM and 3D analysis are considered. This is somehow expected as PMs are less a ff ected than the LOS velocities by the smoothing introduced by the superposition of stars located at di ff erent cluster-centric distances and attaining di ff erent rotation velocities. In addition, the 3D analysis accounts for any rotation axis inclination and projection e ff ects. Finally, we observe that for all velocity components, in dynamically young clusters the SP is characterized by larger α and ω 3D values (i.e., more rapid rotation) than that observed for the FP at similar N h, and it shows a more rapid decline than the FP as a function of N h. In fact, a Spearman correlation test gives a probability P spear larger than ∼ 99% of correlations between ( α SP )PM or ω SP 3D and N h, while probabilities are smaller when either the LOS component or the FP is considered. We note here that by using the approach described in (Curran 2015) we have verified that the results of the Spearman rank correlation tests performed in our analysis and reported in the following are robust against possible outliers and errors associated to the adopted kinematic parameters.  \nTo better highlight such a di ff erential behavior, the upper panels of Figures 9, 10 and 11 show the di ff erence between the rotation strength of the SP and FP as given by ( α SP -α FP ) and ( ω SP 3D -ω FP 3D ). Admittedly, nor ( α SP -α FP ) or ( ω SP 3D -ω FP 3D ) show striking variations in the dynamical age range sampled by the target clusters (2 < N h < 25). In fact, Spearman rank correlation tests give probabilities of correlation of P spear ∼ 90% ( ∼ 95% in  \nFig. 9. Bottom and middle panels show the distribution of the ( α )LOS parameter for the FP and SP as a function of the dynamical age Nh for all clusters in the sample (grey circles). The upper panel shows instead the distribution of the rotation di ff erences ( α SP -α SP )LOS. The dashed-lines represent the linear best-fit to the GC distribution. In all panels, the starry symbols refer to the results obtained for the stacked analysis on the dynamically young (green) and old (orange) samples. The size of the star matches the amplitude of the errorbars.  \n<!-- image -->  \nthe 3D case). Interestingly however, a negative trend between the rotation strength di ff erences and N h is consistently observed in all velocity components. In fact, both ( α SP -α FP ) and ( ω SP 3D -ω FP 3D ) show positive values for dynamically young GCs and then they progressively approach 0 for dynamical older clusters, meaning that FP and SP rotate at the same velocity. The nice agreement between the results obtained in the three analysis definitely provides support to the fact that there is a real correlation between the SP and FP rotation strength di ff erences and N h, and that in general SP shows a more rapid rotation than the FP at dynamically young ages. These results represent the first observational evidence of the link between MP rotation patterns and clusters' long-term dynamical evolution.  \nWe do not find any significant di ff erence between the MP rotation axis orientation both for the LOS and 3D analysis. In fact, the mean di ff erence between the best-fit PA 0 values of the SP and FP is -2 · ± 19 · , and those between θ 0 and i are 4 · ± 24 · and 2 · ± 12 · , respectively.  \nAs an additional way to analyze the data and search for possible trends, we divided the clusters in two sub-groups according to their dynamical ages. In particular, we defined a group of clusters with N h < 8 and a complementary one ( N h > 8) including all the remaining GCs. In this way the two sub-groups turn out to be populated by the same number of systems. Within each sub-group and sub-population, we then stacked all the available kinematic information after normalizing the cluster-centric dis-  \nFig. 10. Same as in Figure 9, but now for ( α ) PM (Section 5.1).  \n<!-- image -->  \ntances to the cluster's R h (from Harris 1996), the velocities to the central velocity dispersion in a given velocity component (as obtained by the analysis described in Section 3) and rotating all clusters to have the same PA 0 (Section 3.1). In this way, we could jointly compare the behavior of MPs for multiple clusters at once, thus increasing the number of stars that can be used to study the kinematics of each sub-population and narrowing down the uncertainties on the derived kinematic parameters The kinematics of MPs in the two stacked samples was then analyzed by following the same approach described in Section 3 and previously adopted for single GCs. Figure 12 shows the results of the kinematic analysis for the two stacked samples for both the LOS and TAN components (lower and upper panels, respectively). In both cases, a significant di ff erence between the FP and SP rotation profiles is observed for the dynamically young ( N h < 8)  \nstacked sample, while they almost disappear for the dynamically old GCs. We then derived the same kinematic parameters described by Equation 9 for a direct comparison with single GCs. Results are shown in Figures 9 and 10 by the two starry symbols. Both in the LOS and TAN components and for each subpopulations, results are fully consistent with the general trend described by single clusters. Interestingly, the reduced uncertainties strengthen the significance of the observed di ff erences discussed above. In particular, the stacked analysis shows that the observed trends between α FP , α SP and N h along both the LOS and TAN components are significant at a large confidence level ( P stacked > 5 σ ). Also, while the ( α SP -α FP ) di ff erence between dynamically young and old clusters is only marginally significant along the LOS component, it turns out to be significant at a ∼ 6 σ level for the tangential component.  \nFig. 11. As in Figures 9 and 10, but now for the ω 3D parameter (Section 5.1).  \n<!-- image -->  \nUnfortunately, we could not apply the 3D rotation analysis (Section 3.3) to the stacked samples as it is not possible to report all clusters to the same values of θ 0 and i .", '5.2. Ellipticity': 'In general a rotating system is also expected to be flattened in the direction perpendicular to the rotation axis (Chandrasekhar 1969). Under the assumption that GCs can be described by the same dynamical model, such as an isotropic oblate rotator (e.g., Varri & Bertin 2012), stronger rotation would be expected in more flattened systems. However, various e ff ects can dilute a possible correlation, the most important ones being anisotropies, inclination e ff ects or tidal forces from the Milky Way (see van den Bergh 2008, for an estimate of the impact of the latter). Nev-  \nertheless, Fabricius et al. (2014) and Kamann et al. (2018) were able to reveal a correlation between cluster rotation and ellipticity in a sample of Galactic GCs (see also Lanzoni et al. 2018a; Dalessandro et al. 2021a; Leanza et al. 2023 for similar analysis on specific clusters).  \nFollowing on those results, we searched for any link between MP rotation and ellipticity for all clusters in our sample. The results of such a comparison for population TOT are reported in Section B. Interestingly, we find a clear and significant correlation between these two quantities (Spearman rank correlation test gives probabilities ∼ 99 . 9%), in nice agreement with previous results by Fabricius et al. (2014) and Kamann et al. (2018).  \nHere we explore in some detail whether also the MP rotation patterns are possibly linked to their morphological 2D spatial distributions. We compared the ellipticity estimates obtained as  \n<!-- image -->  \nFig. 12. Best-fit results of the kinematic analysis of the dynamically young ( N h < 8) and old ( N h > 8) stacked samples. FP is in red and SP in blue. The lower panels refer to the LOS rotation, while the upper row shows the results for the TAN velocity component.  \n<!-- image -->  \ndescribed in Section 4 with the ( α )LOS MP values in Figure 13. Both populations show a positive correlation between ( α )LOS and ϵ , with Spearman rank correlation probabilities P spear ∼ 80% and larger than 99 . 9% for the FP and SP, respectively. In detail, the FP shows a pretty flat distribution of ( α FP )LOS up to ellipticity values ϵ ∼ 0 . 2. Then, ( α FP )LOS starts to increase almost linearly with ϵ . As discussed in Section 5.1, the SP tends to show larger values of rotation than the FP. Likely driven by such a stronger rotation, in the upper panel of Figure 13, we observe that ( α SP )LOS follows a nicely linear correlation with ϵ for the entire range of ellipticity values sampled by the target GCs.  \nFollowing the analysis by Fabricius et al. (2014) and Kamann et al. (2018) we computed the di ff erences between the PA values obtained in Section 4 for the 2D stellar spatial distribu-  \ntion and the best-fit rotation axis position angles ( PA 0). Interestingly, while the distribution of the di ff erence is pretty scattered we find that the average value for the systems in our sample is ∼ 85 · , thus implying that the stellar density distribution is on average flattened in the direction perpendicular to the rotation axis. This behavior is in general agreement with what is expected for a rotating system and it is qualitatively consistent with what predicted, for example, by the models introduced by Varri & Bertin (2012) and previously found in other observational studies (e.g., Bianchini et al. 2013; Bellini et al. 2017; Lanzoni et al. 2018a; Dalessandro et al. 2021a; Leanza et al. 2022).  \nFig. 13. Distributions of the ( α )LOS parameter for the FP and SP (lower and upper panel respectively) as a function of the best-fit ellipticity values (Section 4) for the two sub-populations.  \n<!-- image -->', '5.3. Dynamical simulations': "To conclude the discussion about the rotational properties of MPs in our target GCs, we briefly present the results of a set of N -body simulations aimed at exploring the evolution of rotating MP clusters. A full discussion and detailed description of the results of these simulations will be presented in White et al. (in prep.), here we just report the evolutionary path followed by the α parameter introduced in this paper to trace the strength of rotation of FP and SP stars and their di ff erence. In our simulations we focus only on the long-term dynamics driven by the e ff ects of two-body relaxation for star clusters evolving in the external tidal field of their host galaxy. Each system starts with 10 5 stars with masses between 0.1 and 1 M ⊙ distributed according to a Kroupa (2001) stellar initial mass function. Our systems start with an equal number of FP and SP stars; following the general properties emerging from a few studies of the formation of SP stars in rotating clusters (see e.g., Bekki 2010, 2011; Lacchin et al. 2022) the SP is initially more centrally concentrated and more rapidly rotating than the FP. The two populations rotate around a common axis. To explore the interplay between internal dynamics and the e ff ects due to the external tidal field we have explored models with di ff erent angles δ between the internal rotation axis and the rotation axis of the cluster orbital motion around the center of the host galaxy. In particular, we have explored systems with values of δ equal to 0 · , 45 · , 90 · , and 180 · . The simulations were run with the NBODY6++GPU code (Wang et al. 2015). In Fig. 14 we show the time 5 evolution of α FP , α SP and ( α SP -α FP ) for these models using rotational velocity profiles calculated for both the LOS (left panel) and the TAN component (right panel). For these plots we adopt an ideal line of sight perpendicular to the cluster angular momentum or parallel to it. We emphasize  \nthat these simulations are not aimed at a detailed comparison with observations but they rather serve as a guide to illustrate the extent of the e ff ects of dynamical processes on the initial di ff erences between the rotational kinematics of the FP and SP populations. We also reiterate that our simulations are focussed on the e ff ects of two-body relaxation and do not include early dynamical phases such as those during which a star cluster responds to the mass loss due to stellar evolution, which can have an e ff ect on the sub-populations' dynamical di ff erences (see e.g. Vesperini et al. 2021; Sollima 2021).  \nBy construction, the α values derived for the SP are significantly larger by about a factor of 2-3 than those of the FP. The results of our simulations show that the e ff ects of two-body relaxation lead to a rapid and significant reduction of the initial di ff erence between the FP and the SP rotation in the first 2-3 half-mass relaxation times reaching values of ( α SP -α FP ) similar to those found in our observational analysis. Then at later dynamical ages, ( α SP -α FP ) keeps decreasing at a lower pace and it progressively approaches values close to 0 around 10 relaxation times, when FP and SP stars rotate at the same velocity. We note that, as already discussed in Section 5.1 and in agreement with what was found in the observations, the rotation strength for both the FP and SP along with their di ff erence is stronger when a PMlike projection is considered (Fig. 14). The behavior described by the simulations is in nice general agreement with the observed trends (Figures 9, 10 and 11). Such an agreement strongly suggests that both the rotational di ff erences and the mild trend between the rotation strength and the dynamical age revealed by our observational analysis are consistent with those expected for the long-term dynamical evolution of GCs born with a SP initially rotating more rapidly than the FP.", '6.1. Observations': "Weinvestigated the di ff erences in the mean anisotropy properties of MPs by defining the following parameter:  \nαβ = Z 1 0 ( β ( Rl ) X ) dRl (11)  \nhere Rl represents the distance from the cluster center normalized to the value of the clusters' Jacobi radius (from Baumgardt et al. 2019) to allow a meaningful comparison among different systems. From a geometrical point of view, Equation 11 represents the area subtended by the best-fit radially normalized anisotropy profile (as defined in Section 3.2). By definition, if a sub-population is radially anisotropic, αβ will have positive values, while it will attain negative values in case of tangential anisotropy. Di ff erences between the anisotropy properties of MPs are given by ( α S P β -α FP β ). As before, errors on αβ were obtained by propagating the posterior probability distributions obtained from the MCMC analysis for the best-fit anisotropy profiles.  \nSimilarly to the rotation analysis, Figure 15 shows the variation of αβ as a function of the dynamical age for the FP and SP in the bottom and middle panel, while the distribution of ( α S P β -α FP β ) is shown in the top panel. α FP β shows a flat distribution around ∼ 0 at any N h, thus suggesting that the FP is on average isotropic for all GCs in the sample, independently on their dynamical ages. On the contrary, the α S P β shows a non negligible variation as a function of time. In fact, for dynamically young  \n<!-- image -->  \nFig. 14. Time evolution of the rotational parameters α FP , α SP , and their di ff erences (see Section 5.1) for the simulations described in Section 5.3. The blue line corresponds to the model with δ = 0 · , the pink, cyan and purple correspond to δ = 45 · , 90 · and 180 · , respectively.  \n<!-- image -->  \nclusters, α S P β is positive (i.e., SP is radially anisotropic) and it becomes isotropic only for the dynamically older systems. To better highlight such a trend, we can use the two GC sub-groups defined in Section 5.1. At a face value, GCs with N h < 8 (dynamically young) have a mean α S P β value of 0 . 40 ± 0 . 20, while for dynamically older clusters the mean value is -0 . 09 ± 0 . 23 (dashed lines in Figure 15). Such a di ff erence is very nicely described by the comparison between the best-fit results obtained for MPs in the two stacked samples as shown in Figure 16. Indeed, when the anisotropy analysis is performed on the two stacked sub-samples, we find that α S P β = 0 . 23 ± 0 . 03 at dynamically young ages, while for dynamically older clusters we obtain α S P β = 0 . 04 ± 0 . 03, thus implying that the di ff erence in the average SP orbital anisotropy is significant at a ∼ 5 σ level. For comparison, we computed the same quantities for the FP on the stacked clusters. As expected, in this case mean values are virtually indistinguishable within the errors. The di ff erent behavior between the average orbital properties of FP and SP stars is further highlighted by the distribution of ( α S P β -α FP β ). In dynamically young GCs SP is systematically more radial anisotropic than the FP (mean di ff erence being 0 . 44 ± 0 . 31), then they both become compatible with being isotropic at later dynamical stages ( -0 . 08 ± 0 . 37). The analysis on the two stacked sub-samples provides values a ( α S P β -α FP β ) = 0 . 28 ± 0 . 07 and 0 . 09 ± 0 . 07 for the dynamically young and old groups, respectively, thus implying that there is a real orbital velocity distribution di ff erence between FP and SP stars with a significance of ∼ 3 . 5 σ .", '6.2. Dynamical models': "The di ff erence between the anisotropy of the FP and SP populations is consistent with the predictions of a number of theoretical investigations of the dynamics of MP clusters. These studies find that in clusters where the SP forms more centrally concen-  \ntrated than the FP, the SP is generally characterized by a more radially anisotropic velocity distribution than the FP (see e.g., Bellini et al. 2015; Hénault-Brunet et al. 2015; Tiongco et al. 2019; Vesperini et al. 2021; Sollima 2021). In order to illustrate the expected strength and evolution of the di ff erences between the anisotropy of the FP and SP populations, as measured by the parameter ( α S P β -α FP β ) introduced in this paper, we show in Figure 17 the time evolution of α FP β , α S P β and ( α S P β -α FP β ) for three Monte Carlo simulations following the dynamical evolution of MP clusters with di ff erent initial conditions. The Monte Carlo simulations were run with the MOCCA code (Hypki & Giersz 2013; Giersz et al. 2013; see also Vesperini et al. 2021; Hypki et al. 2022 for the specific case of GCs with MPs). A full presentation of the kinematics and phase space properties of these systems will be presented in a separate paper (Aros et al. in prep.). The three models for which we present the anisotropy evolution in Fig. 17 all start with a ratio of the SP to total mass equal to 0 . 2. The total number of stars is equal to N = 10 6 (hereafter we use the id. N1M for these model) or 0 . 5 × 10 6 (with id. N05M ) and stars have been extracted from a Kroupa (2001) initial stellar mass function between 0.1 and 100 M ⊙ . The SP is initially more concentrated than the FP: for two models the initial ratio for the FP to SP half-mass radii is equal to 20 ( N1Mr20 and N05Mr20 ) and is equal to 10 for one model ( N1Mr10 ). Both populations are initially characterized by a non-rotating, isotropic velocity distribution. All models include the e ff ects of mass loss due to stellar evolution, two-body relaxation, binary interactions and a tidal truncation (see Vesperini et al. 2021 for further details on the evolutionary processes included in these simulations). As shown in Figure 17, in all models the SP is characterized by a more radially anisotropic velocity distribution than the FP. Even if both populations are initially isotropic, the development of a stronger radial anisotropy for the SP is the consequence of its initial more centrally concentrated spatial distribution (see e.g., Tiongco et al. 2019; Vesperini et al. 2021). The di ff erences in anisotropy  \nFig. 15. The bottom and middle panels show the distribution of the anisotropy parameter αβ for the FP and SP sub-populations of GCs in the sample (grey circles) as a function of N h. The top panel, shows the di ff erential trend α SP β -α FP β . The dashed lines represent the mean of the observed αβ values derived for the dynamically young ( N h < 8) and old ( N h > 8) GCs. The starry symbols are the result of the analysis on the two stacked samples and their sizes correspond to the errorbars.  \n<!-- image -->  \nbetween the two populations gradually decrease as the systems evolves. As shown by the results of our simulations, the extent of the di ff erences between the FP and SP anisotropy depends on the initial relative concentration of the two populations being smaller for the model where the initial ratio of the FP to SP half-mass radii equal to 10 ( N1MR10 ) than that for models where this ratio is initially equal to 20 ( N1MR20 ). While also in this case we point out that these simulations have not been specifically designed to provide a detailed fit to observations but rather to provide guidance in the general interpretation of the observational results, it is interesting to note that the extent of the di ff erence between the anisotropy of the FP and the SP popula-  \ntions as measured by the ( α S P β -α FP β ) and its time evolution is generally consistent with the findings of our observational analysis (Figure 15). The general agreement between observations and simulations lends further support to the interpretation of the observed di ff erence between the FP and SP anisotropy and its variation with the clusters' dynamical age as the manifestation of the di ff erent dynamical paths followed by subsystems with di ff erent initial spatial distributions.  \nFig. 16. Best-fit results of the anisotropy analysis of MPs in the two stacked samples. The upper panels shows the best-fit anisotropy profiles of the FP (red) and SP (blue) sub-populations as obtained for the stacked sample of dynamically young GCs, while the lower panel refers to dynamically old systems.  \n<!-- image -->", '7. Comparison with the literature': "In Appendix C we report on a quantitative comparison between the results obtained in the present work and those available in literature for the TOT population in each cluster. Here we detail on the comparison with previous works focusing on the kinematics of MPs. We stress that the full 3D kinematic analysis presented in this paper is the first ever obtained for MPs. Hence, in the following our comparison will be limited to studies considering a single velocity component.  \nOur sample has 6 GCs in common with the recent analysis by Martens et al. (2023) based on MUSE LOS velocities. The authors were able to find MP best-fit solutions for both the velocity dispersion and rotation profiles for three of them (namely 47 Tuc,  \nNGC5904 and NGC 6093), while for NGC 3201 and NGC 6254 they report only conservative upper limits for the MP rotation amplitudes, and for NGC 1904 they provide information only for the TOT population. Figure 18 shows the distributions of the di ff erences in terms of σ 0 (which is defined as σ max in Martens et al. 2023), A rot ( v max in Martens et al. 2023) and R peak for the FP and SP sub-populations. A good agreement is observed both for the MP central velocity dispersion values (left panel of Figure 18) and the rotation amplitude (middle panel), with the most discrepant value being that corresponding to the rotation of the SP in 47 Tuc. Within the errors there is also a reasonable match between the values of R peak. We note however that, while the sample is certainly small, the estimates of R peak by Martens et al.  \nTable 2. MP kinematic parameters describing the radial anisotropy and rotation along the LOS, PM and 3D components.  \n(2023) tend to be slightly smaller than those derived in this work. This might be somehow linked to the smaller radial coverage of the Martens et al. (2023) analysis. In fact, all values derived by Martens et al. (2023) for the clusters in common are located well outside the MUSE field of view and therefore they might be only partially constrained by their analysis.  \nThe sample analyzed in this work counts also 4 GCs in common with Cordoni et al. (2020), namely 47 Tuc, NGC 288, NGC 5904 and NGC 6254. Our results are in agreement with  \ntheirs in that 47 Tuc and NGC 5904 are the systems showing the larger rotation among the clusters in common. However, at odds with their results, we also find that in both GCs the SP shows a larger rotation (as inferred both by A rot and α values) than the FP. Also, within the uncertainties, we do not find evidence of any significant mis-alignment of the FP and SP rotation curves for these systems, neither in terms of position angles and inclination, as constrained by both the LOS and the full 3D analysis.  \nFig. 17. Time evolution of the anisotropy parameters defined in Section 6.1 for the FP (bottom panel), the SP (middle panel) and their di ff erence (top panel) in the three Monte Carlo simulations described in Section 6.2 (N1Mr10 blue line, N1M purple line; N05M pink line). Time is normalized by the half-mass relaxation time.  \n<!-- image -->  \nFinally, our analysis is in nice agreement with the results presented by Cordero et al. (2017) for the MP rotation patterns of M13.  \nAs for the anisotropy, the results obtained in this work are in qualitative agreement with those found for the GCs M 13, NGC 2808 and 47 Tuc by Richer et al. (2013); Bellini et al. (2015); Milone et al. (2018). More in general, they are in nice qualitative agreement with those recently obtained in the innermost regions (R < 100 '' ) for a sample of Galactic GCs by Libralato et al. (2023) by using HST PMs to search for average di ff erences between the anisotropy profiles of MPs.", '8. Summary and conclusions': "Wehave presented the first self-consistent 3D kinematic analysis of MPs for a sample of Galactic GCs. The study targets 16 systems spanning a broad range of dynamical ages (2 < N h < 25) and it is based on a large and mostly homogeneous observational dataset securing several hundreds of accurate LOS velocities and PMs for each cluster and sampling virtually their entire extension.  \nOur study is mainly focused on the analysis of the MP rotation along the three velocity components and the anisotropy of their velocity distributions. The adopted approach is aimed at providing new insights into the long term evolution of the kinematic properties of MPs (and their di ff erences) for the entire sample of GCs and for the entire dynamical age covered by our analysis instead of focusing on the kinematic di ff erences in specific clusters. To this aim, starting from the observed velocity distributions we defined a few key quantities, to quantitatively and homogeneously compare the results obtained for all the observed GCs.  \nOur analysis provides the first observational determination of the dynamical path followed by MP kinematic properties during  \ntheir long-term evolution. The main observational results we find can be schematically summarized as follows.  \n- -We observe evidence of di ff erential rotation between MPs with the SP preferentially rotating more rapidly than the FP. This result is consistent (although with di ff erent amplitudes) along both the LOS and TAN velocity components, as well as in the full 3D analysis. In all GCs in our sample we find that the rotation position and inclination angles are consistent within the uncertainties between FP and SP.\n- -The strength of the rotation signal of both FP and SP subpopulations nicely correlate with the ellipticity values derived for the two sub-populations. In addition, we find that the rotation axis position angle is typically perpendicular to the ellipses major axis.\n- -The di ff erence in the rotation strength between MPs is mildly anti-correlated with the cluster dynamical age. In particular, di ff erences are larger for dynamically young clusters and they become progressively indistinguishable as dynamical evolution proceeds.\n- -Stars belonging to di ff erent sub-populations show di ff erent average velocity distribution properties. FPs are characterized by isotropic velocity distributions at any dynamical age probed by our sample. On the contrary, the SP of dynamically young GCs have a radially anisotropic velocity distribution which then becomes isotropic in more advanced evolutionary stages.  \nThe combination of these results with the analysis of the MP radial distributions of a representative sample of GCs carried out in our previous paper (Dalessandro et al. 2019), provides a full picture of the present-day kinematic and structural properties of MPs in GCs and of their evolution. The comparison with dynamical models following the long-term evolution of MPs in GCs, suggests that these properties, and their evolution with dynamical age are in general good agreement with those expected in clusters forming with a SP subsystem initially more centrally concentrated and more rapidly rotating than the FP (see e.g., D'Ercole et al. 2008; Bekki 2010; Calura et al. 2019; Lacchin et al. 2022). In turn, this would possibly suggest (see e.g., Hénault-Brunet et al. 2015 and discussion in Martens et al. 2023) that GCs experienced multiple events of star formation during their early phases of evolution, with the rotation properties being the more stringent discriminating factors. In fact, according to multi-epoch formation models, the SP is expected to form a low-mass, more centrally concentrated and more rapidly rotating stellar sub-system than the FP. In such a configuration, dissipative accretion processes of material ejected by FP stars and angular momentum conservation in sub-systems with different initial concentration can produce a larger initial rotation of the SP sub-system (Bekki 2011; Hénault-Brunet et al. 2015; Tiongco et al. 2017). Interestingly, it has been shown (e.g., Bekki 2011) that even if only a very small fraction of the kinetic energy of the FP is in the form of bulk rotation energy, SP stars can acquire a much stronger rotation than what remains in the FP.  \nIt is important to note that, as shown in this paper and in a number of previous studies, early and long-term evolution can significantly reduce the initial di ff erences between the FP and SP rotational velocities making them virtually indistinguishable in dynamical old systems. The combination of these physical e ff ects with the observational uncertainties arising from the limited available stellar samples, partial cluster coverage and possible biases introduced by the di ff erent rotation inclination angles, can make extremely di ffi cult to capture present-day kinematic  \n<!-- image -->  \n<!-- image -->  \nFig. 18. Comparison between the best-fit σ 0, A rot and R peak values obtained for MPs in the present-work and by Martens et al. (2023). Blue dots refer to SP and red ones to FP results.  \n<!-- image -->  \ndi ff erences between MPs, in particular in the LOS velocity component. As a consequence, it is important to use some caution in drawing conclusions about the physical mechanisms at the basis of GC formation and early evolution based on the presentday kinematic and structural properties of individual systems or small samples.  \nAn homogeneous combination between data obtained with multi-object spectrographs and Gaia (as presented in this paper) which are particularly sensitive to the intermediate and large cluster-centric distances in GCs, with Integral Field Unit LOS RVs and HST absolute PMs sampling more e ffi ciently the clusters' innermost regions (see for example Martens et al. 2023; Libralato et al. 2023), along with an increase of the cluster sample size, represents a natural and promising next step that could potentially enable the exploration of additional physical ingredients at play.  \nThe results of this study clearly demonstrates that significant advances in our understanding of cluster formation and early evolution is only possible through a multi-faceted and multidiagnostic approach and by combining state-of-art observations and simulations.  \nAcknowledgements. The authors thank the anonymous referee for the careful reading of the paper and the useful comments that improved the presentation of this work. E.D. acknowledges financial support from the Fulbright Visiting Scholar program 2023. E.D. and A.D.C. are also grateful for the warm hospitality of the Indiana University where part of this work was performed. E.D. and M.C. acknowledge financial support from the project Light-on-Dark granted by MIUR through PRIN2017-2017K7REXT. E.V. acknowledges support from NSF grant AST-2009193. E.V. acknowledges also support from the John and A-Lan Reynolds Faculty Research Fund. The research activities described in this paper have been co-funded by the European Union - NextGenerationEU within PRIN 2022 project n.20229YBSAN - Globular clusters in cosmological simulations and in lensed fields: from their birth to the present epoch.", 'References': "Aarseth, S. J. 2003, Gravitational N-Body Simulations, by Sverre J. Aarseth, pp. 430. ISBN 0521432723. Cambridge, UK: Cambridge University Press, November 2003., 430  \nBastian, N., Lamers, H. J. G. L. M., de Mink, S. E., et al. 2013, MNRAS, 436, 2398  \n- Bastian, N., & Lardo, C. 2015, MNRAS, 453, 357\n- Bastian, N. & Lardo, C. 2018, ARA&A, 56, 83. doi:10.1146 / annurev-astro081817-051839  \nBaumgardt, H., Hilker, M., Sollima, A., et al. 2019, MNRAS, 482, 5138 Bekki, K. 2010, ApJ, 724, L99  \n- Bekki, K. 2011, MNRAS, 412, 2241. doi:10.1111 / j.1365-2966.2010.18047.x\n- Bellazzini, M., Bragaglia, A., Carretta, E., et al. 2012, A&A, 538, A18  \nBellini, A., Vesperini, E., Piotto, G., et al. 2015, ApJ, 810, L13. doi:10.1088 / 2041-8205 / 810 / 1 / L13  \n- Bellini, A., Milone, A. P., Anderson, J., et al. 2017, ApJ, 844, 164. doi:10.3847 / 1538-4357 / aa7b7e\n- Bianchini, P., Varri, A. L., Bertin, G., et al. 2013, ApJ, 772, 67. doi:10.1088 / 0004-637X / 772 / 1 / 67\n- Bianchini, P., van der Marel, R. P., del Pino, A., et al. 2018, MNRAS, 481, 2125. doi:10.1093 / mnras / sty2365\n- Bragaglia, A., Carretta, E., D'Orazi, V., et al. 2017, A&A, 607, A44\n- Cadelano, M., Saracino, S., Dalessandro, E., et al. 2020, ApJ, 895, 54. doi:10.3847 / 1538-4357 / ab88b3\n- Cadelano, M., Dalessandro, E., Salaris, M., et al. 2022, ApJ, 924, L2. doi:10.3847 / 2041-8213 / ac424a\n- Cadelano, M., Pallanca, C., Dalessandro, E., et al. 2023, A&A, 679, L13. doi:10.1051 / 0004-6361 / 202347961\n- Cadelano, M., Dalessandro, E., & Vesperini, E. 2024, arXiv:2402.09514. doi:10.48550 / arXiv.2402.09514\n- Carretta, E., Bragaglia, A., Gratton, R., et al. 2009, A&A, 505, 139. doi:10.1051 / 0004-6361 / 200912097\n- Chung, C., Yoon, S.-J., & Lee, Y.-W. 2011, ApJ, 740, L45\n- Cordero, M. J., Hénault-Brunet, V., Pilachowski, C. A., et al. 2017, MNRAS, 465, 3515\n- Cordoni, G., Milone, A. P., Mastrobuono-Battisti, A., et al. 2020, ApJ, 889, 18. doi:10.3847 / 1538-4357 / ab5aee\n- Della Croce, A., Dalessandro, E., Livernois, A., et al. 2023, A&A, 674, A93. doi:10.1051 / 0004-6361 / 202346095\n- Dotter, A., Sarajedini, A., Anderson, J., et al. 2010, ApJ, 708, 698. doi:10.1088 / 0004-637X / 708 / 1 / 698\n- Calura, F., D'Ercole, A., Vesperini, E., Vanzella, E., & Sollima, A. 2019, arXiv:1906.09137\n- Chandrasekhar, S. 1969, ApJ, 157, 1419. doi:10.1086 / 150161\n- Curran, P. A. 2015, Astrophysics Source Code Library. ascl:1504.008\n- D'Antona, F., Vesperini, E., D'Ercole, A., et al. 2016, MNRAS, 458, 2122\n- Dalessandro, E., Salaris, M., Ferraro, F. R., et al. 2011, MNRAS, 410, 694. doi:10.1111 / j.1365-2966.2010.17479.x\n- Dalessandro, E., Schiavon, R. P., Rood, R. T., et al. 2012, AJ, 144, 126\n- Dalessandro, E., Ferraro, F. R., Massari, D., et al. 2013, ApJ, 778, 135\n- Dalessandro, E., Massari, D., Bellazzini, M., et al. 2014, ApJ, 791, L4. doi:10.1088 / 2041-8205 / 791 / 1 / L4\n- Dalessandro, E., Lapenna, E., Mucciarelli, A., et al. 2016, ApJ, 829, 77\n- Dalessandro, E., Cadelano, M., Vesperini, E., et al. 2018, ApJ, 859, 15\n- Dalessandro, E., Mucciarelli, A., Bellazzini, M., et al. 2018a, ApJ, 864, 33\n- Dalessandro, E., Lardo, C., Cadelano, M., et al. 2018, A&A, 618, A131. doi:10.1051 / 0004-6361 / 201833650\n- Dalessandro, E., Cadelano, M., Vesperini, E., et al. 2019, ApJ, 884, L24\n- Dalessandro, E., Raso, S., Kamann, S., et al. 2021a, MNRAS, 506, 813. doi:10.1093 / mnras / stab1257\n- Dalessandro, E., Varri, A. L., Tiongco, M., et al. 2021b, ApJ, 909, 90. doi:10.3847 / 1538-4357 / abda43\n- Decressin, T., Meynet, G., Charbonnel, C., Prantzos, N., & Ekström, S. 2007, A&A, 464, 1029\n- de Mink, S. E., Pols, O. R., Langer, N., et al. 2009, A&A, 507, L1. doi:10.1051 / 0004-6361 / 200913205\n- D'Ercole, A., Vesperini, E., D'Antona, F., McMillan, S. L. W., & Recchi, S. 2008, MNRAS, 391, 825  \nEinsel, C. & Spurzem, R. 1999, MNRAS, 302, 81. doi:10.1046 / j.13658711.1999.02083.x  \n- Fabricius, M. H., Noyola, E., Rukdee, S., et al. 2014, ApJ, 787, L26. doi:10.1088 / 2041-8205 / 787 / 2 / L26  \nFerraro, F. R., Mucciarelli, A., Lanzoni, B., et al. 2018b, ApJ, 860, 50  \nForeman-Mackey, D., Hogg, D. W., Lang, D., Goodman, J. 2013, PASP, 125,  \n306  \nGaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1. doi:10.1051 / 0004-6361 / 202243940  \nGieles, M., Charbonnel, C., Krause, M. G. H., et al. 2018, MNRAS,  \n- Giersz, M., Heggie, D. C., Hurley, J. R., et al. 2013, MNRAS, 431, 2184. doi:10.1093 / mnras / stt307  \nGoldsbury, R., Richer, H. B., Anderson, J., et al. 2010, AJ, 140, 1830. doi:10.1088 / 0004-6256 / 140 / 6 / 1830  \n- Gratton, R., Bragaglia, A., Carretta, E., et al. 2019, A&A Rev., 27, 8. doi:10.1007 / s00159-019-0119-3  \nHarris, W. E. 1996, AJ, 112, 1487  \nHénault-Brunet, V., Gieles, M., Agertz, O., et al. 2015, MNRAS, 450, 1164. doi:10.1093 / mnras / stv675  \n- Hypki, A. & Giersz, M. 2013, MNRAS, 429, 1221. doi:10.1093 / mnras / sts415 Hypki, A., Giersz, M., Hong, J., et al. 2022, MNRAS, 517, 4768. doi:10.1093 / mnras / stac2815\n- Hong, J., Kim, E., Lee, H. M., et al. 2013, MNRAS, 430, 2960. doi:10.1093 / mnras / stt099  \nKamann, S., Dalessandro, E., Bastian, N., et al. 2020, MNRAS, 492, 966  \nKamann, S., Husser, T.-O., Dreizler, S., et al. 2018, MNRAS, 473, 5591.  \ndoi:10.1093  \n/  \nmnras  \n/  \nstx2719  \nKamlah, A. W. H., Spurzem, R., Berczik, P., et al. 2022, MNRAS, 516, 3266. doi:10.1093 / mnras / stac2281  \nKüpper, A. H. W., Maschberger, T., Kroupa, P., et al. 2011, MNRAS, 417, 2300. doi:10.1111 / j.1365-2966.2011.19412.x  \nKroupa, P. 2001, MNRAS, 322, 231. doi:10.1046  \n/  \nj.1365-8711.2001.04022.x  \nLacchin, E., Calura, F., Vesperini, E., et al. 2022, MNRAS, 517, 1171.  \ndoi:10.1093  \n/  \nmnras  \n/  \nstac2328  \nLardo, C., Salaris, M., Bastian, N., et al. 2018, arXiv:1805.09599  \nLanzoni, B., Ferraro, F. R., Mucciarelli, A., et al. 2018a, ApJ, 861, 16  \nLanzoni, B., Ferraro, F. R., Mucciarelli, A., et al. 2018b, ApJ, 865, 11  \nLarsen, S. S., Strader, J., & Brodie, J. P. 2012, A&A, 544, L14  \nLarsen, S. S., Brodie, J. P., Wasserman, A., & Strader, J. 2018, arXiv:1801.03140 Leanza, S., Pallanca, C., Ferraro, F. R., et al. 2022, ApJ, 929, 186. doi:10.3847 / 1538-4357 / ac5d4e  \nLeanza, S., Pallanca, C., Ferraro, F. R., et al. 2023, ApJ, 944, 162. doi:10.3847 / 1538-4357 / acaaa6  \nLeitinger, E., Baumgardt, H., Cabrera-Ziri, I., et al. 2023, MNRAS, 520, 1456.  \ndoi:10.1093  \n/  \nmnras  \n/  \nstad093  \n- Libralato, M., Bellini, A., Piotto, G., et al. 2019, ApJ, 873, 109. doi:10.3847 / 1538-4357 / ab0551  \nLibralato, M., Vesperini, E., Bellini, A., et al. 2023, ApJ, 944, 58. doi:10.3847 / 1538-4357 / acaec6  \nLivernois, A. R., Vesperini, E., Varri, A. L., et al. 2022, MNRAS, 512, 2584. doi:10.1093 / mnras / stac651  \nLynden-Bell, D. 1967, MNRAS, 136, 101  \nMartens, S., Kamann, S., Dreizler, S., et al. 2023, A&A, 671, A106. doi:10.1051 / 0004-6361 / 202244787  \nMartocchia, S., Cabrera-Ziri, I., Lardo, C., et al. 2018, MNRAS, 473, 2688  \nMartocchia, S., Niederhofer, F., Dalessandro, E., et al. 2018, MNRAS, 477,  \n4696. doi:10.1093  \n/  \nmnras  \n/  \nsty916  \nMerritt, D. 1985, AJ, 90, 1027. doi:10.1086  \n/  \n113810  \nMilone, A. P., Piotto, G., Renzini, A., et al. 2017, MNRAS, 464, 3636  \nMilone, A. P., Marino, A. F., Mastrobuono-Battisti, A., et al. 2018, MNRAS,  \n479, 5005. doi:10.1093  \n/  \nmnras  \n/  \nsty1873  \nMiocchi, P., Lanzoni, B., Ferraro, F. R., et al. 2013, ApJ, 774, 151  \nMonelli, M., Milone, A. P., Stetson, P. B., et al. 2013, MNRAS, 431, 2126. doi:10.1093 / mnras / stt273  \nMuratov, A. L., & Gnedin, O. Y. 2010, ApJ, 718, 1266  \nMucciarelli, A., Origlia, L., Ferraro, F. R., & Pancino, E. 2009, ApJ, 695, L134 Mucciarelli, A., Dalessandro, E., Massari, D., et al. 2016, ApJ, 824, 73. doi:10.3847 / 0004-637X / 824 / 2 / 73  \nNardiello, D., Piotto, G., Milone, A. P., et al. 2015, MNRAS, 451, 312. doi:10.1093 / mnras / stv971  \nNardiello, D., Libralato, M., Piotto, G. et al. 2018, MNRAS, 481, 3382  \nNiederhofer, F., Bastian, N., Kozhurina-Platais, V., et al. 2017, MNRAS, 465, 4159  \nNitadori, K. & Aarseth, S. J. 2012, MNRAS, 424, 545. doi:10.1111 / j.13652966.2012.21227.x  \nOnorato, S., Cadelano, M., Dalessandro, E., et al. 2023, A&A, 677, A8. doi:10.1051 / 0004-6361 / 202346792  \nOsipkov, L. P. 1979, AZh, 56, 378  \nPedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825. doi:10.48550 / arXiv.1201.0490  \nArticle number, page 22 of 26  \n- Petralia, I., Minniti, D., Fernández-Trincado, J. G., et al. 2024,\n- arXiv:2404.10902. doi:10.48550 / arXiv.2404.10902 Piotto, G., Bedin, L. R., Anderson, J., et al. 2007, ApJ, 661, L53. doi:10.1086 / 518503 Piotto, G., Milone, A. P., Bedin, L. R., et al. 2015, AJ, 149, 91 Plummer, H. C. 1911, MNRAS, 71, 460. doi:10.1093 / mnras / 71.5.460 Pryor, C., & Meylan, G. 1993, Structure and Dynamics of Globular Clusters, 357 Richer, H. B., Heyl, J., Anderson, J., et al. 2013, ApJ, 771, L15. doi:10.1088 / 2041-8205 / 771 / 1 / L15 Saracino, S., Martocchia, S., Bastian, N., et al. 2020, MNRAS, 493, 6060. doi:10.1093 / mnras / staa644 Sbordone, L., Salaris, M., Weiss, A., et al. 2011, A&A, 534, A9 Szigeti, L., Mészáros, S., Szabó, G. M., et al. 2021, MNRAS, 504, 1144. doi:10.1093 / mnras / stab1007 Schiavon, R. P., Caldwell, N., Conroy, C., et al. 2013, ApJL, 776, L7 Sills, A., Dalessandro, E., Cadelano, M., et al. 2019, MNRAS, 490, L67. doi:10.1093 / mnrasl / slz149 Sollima, A., Baumgardt, H., & Hilker, M. 2019, MNRAS, 485, 1460 Sollima, A. 2021, MNRAS, 502, 1974. doi:10.1093 / mnras / stab154 Stetson, P. B., Pancino, E., Zocchi, A., et al. 2019, MNRAS, 485, 3042. doi:10.1093 / mnras / stz585 Tiongco, M. A., Vesperini, E., & Varri, A. L. 2016, MNRAS, 455, 3693. doi:10.1093 / mnras / stv2574 Tiongco, M. A., Vesperini, E., & Varri, A. L. 2017, MNRAS, 469, 683. doi:10.1093 / mnras / stx853 Tiongco, M. A., Vesperini, E., & Varri, A. L. 2019, MNRAS, 1619 van den Bergh, S. 2008, AJ, 135, 1731. doi:10.1088 / 0004-6256 / 135 / 5 / 1731 van Leeuwen, F. 2009, A&A, 497, 209. doi:10.1051 / 0004-6361 / 200811382 Varri, A. L. & Bertin, G. 2012, A&A, 540, A94. doi:10.1051 / 00046361 / 201118300 Vasiliev, E. & Baumgardt, H. 2021, MNRAS, 505, 5978. doi:10.1093 / mnras / stab1475 Vesperini, E., McMillan, S. L. W., D'Antona, F., et al. 2013, MNRAS, 429, 1913. doi:10.1093 / mnras / sts434 Vesperini, E., Hong, J., Giersz, M., et al. 2021, MNRAS, 502, 4290. doi:10.1093 / mnras / stab223 Wang, L., Spurzem, R., Aarseth, S., et al. 2015, MNRAS, 450, 4070.\n- doi:10.1093 / mnras / stv817\n- Zemp, M., Gnedin, O. Y., Gnedin, N. Y., et al. 2011, ApJS, 197, 30. doi:10.1088 / 0067-0049 / 197 / 2 / 30", 'Appendix A: Additional table': 'Table A.1. Best-fit kinematic parameters for the MPs in the target GCs.', 'Appendix B: Incompleteness effects': "We constrained the possible impact of the kinematic samples' size and of their (radially dependent) incompleteness (mainly caused by the intrinsically limited allocation e ffi ciency of multiobject spectrographs) on the results obtained in this work by using the dynamical simulations described in Sections 5 and 6.  \nIn detail, we estimated the completeness ( C ) of the observed kinematic samples and its radial variation as the ratio between the number of RGB stars detected in the photometric catalogs and those with LOS RVs and / or PMs measures within concentric radial annuli at di ff erent cluster-centric distances. While we acknowledge that these estimates represent a lower-limit to the real incompleteness, as photometric catalogs are not fully complete, we note that our targets are RGB stars, which are among the brightest stars in GC CMDs and therefore they are only moderately a ff ected by incompleteness even when ground-based catalogs are adopted.  \nWe then extracted randomly from the simulations subsamples of stars with similar sizes as the observed ones for a large number of times. We applied the derived completeness curves to these sub-samples to make them as similar as possible to the observed catalogs. Finally we run the same kinematic analysis described in Section 3. We find that, while the limited sample sizes and incompleteness have an impact on the overall noise of the observed kinematic profiles and as a consequence on the uncertainties associated to the derived parameters, they do not have a significant impact on the final results.", 'Appendix C: Global kinematics and comparison with the literature': 'Figure C.1 shows the distribution of the α values derived for the entire population along both the LOS and TAN velocity components, and of the ω TOT 3 D , as a function of N h. As expected, both α TOT (for both LOS and TAN ) and the ω TOT 3 D distributions show pretty clear anti-correlations with N h. In fact, while dynamically young GCs attain larger values of rotation parameters, with the only significant exception being M 3, the rotation strength progressively decreases as dynamical evolution proceeds. By performing a Spearman rank correlation test we find that such anti-correlations have a significance of ∼ 98% for the LOS and > 99 . 9% for both the TAN and the 3D analysis, with the 3D case being the most significant. In general, these results are in very good agreement with previous analysis (e.g., Kamann et al. 2018; Sollima et al. 2019) and they further strengthen the conclusion that the present-day cluster rotation is the relic of that imprinted at the epoch of cluster formation, which has been then progressively dissipated via two-body relaxation.  \nFigure C.2 shows the distribution of α TOT LOS with the cluster ellipticity obtained as described in Section 3. As expected (see discussion and references in Section 5.2), a nice correlation between rotation and ellipticity is observed also for the total population in very good agreement with previous findings by Fabricius et al. (2014); Kamann et al. (2018).  \nWhile the focus of this work is on the MP kinematics, nevertheless it is useful to compare the results obtained for the TOT population with those largely available in the literature to have an indication about the general performance of the adopted approach and data-sets.  \nDetailed one-to-one comparisons with recent results obtained in the literature (Bellazzini et al. 2012; Ferraro et al. 2018; Lanzoni et al. 2018a,b; Baumgardt et al. 2019; Sollima et al. 2019) for the TOT population are shown in Figures C.3 and C.4.  \nFig. C.1. Distribution of the three rotation parameters defined in this work for the LOS , PM and 3D velocity components, as a function of the dynamical age ( N h) for the total population (TOT) of GCs in our sample.  \n<!-- image -->  \nFig. C.2. Distribution of the ( α TOT )LOS parameter for the total population as a function of the best-fit ellipticity values derived as described in Section 4.  \n<!-- image -->  \nIn general, a quite good agreement is observed with all the compilations considered here.  \nOur sample has 6 GCs in common with Bellazzini et al. (2012). A nice match is observed both in terms of σ 0 and A rot (top row of Figure C.3) with the only exception being NGC 6171 for which Bellazzini et al. (2012) finds a rotation amplitude 4-5 times larger than the one obtained in this work. Given the estimate by Bellazzini et al. (2012), NGC 6171 would be a very fast rotator, with Arot /σ 0 ∼ 0 . 7. However, it is important to note, that the sample of LOS RVs used by Bellazzini et al. (2012) for this cluster includes only 31 stars in total, resulting the smallest sample of LOS RVs in their analysis. Here we sample the kinematic profile of NGC 6171 with 184 LOS RVs (see Table 1). We note also that NGC 6171 results to have a significantly smaller rotation amplitude (1 . 2 km / s) than what found by Bellazzini et al. (2012) in the analysis by Ferraro et al. (2018) and it is classified as non rotator by Sollima et al. (2019)  \nAs for the comparison with results by Ferraro et al. (2018), we stress that while the spectroscopic sample is largely similar, the adopted kinematic analysis (both the discrete and continuous  \n<!-- image -->  \nFig. C.3. Comparison between the best-fit σ 0 and A rot values obtained for the TOT sample for the clusters in common between the present work and Bellazzini et al. (2012) - B12, Ferraro et al. (2018) - F18 and Petralia et al. (2024) - P24.  \n<!-- image -->  \nones) is significantly di ff erent for the rotation study in particular (see Ferraro et al. 2018 for details). Hence, it is not surprising that while the derived central velocity dispersion values are in excellent agreement for the entire sample (middle row of Figure C.3), the distribution of di ff erences for A rot is more scattered, while still showing a reasonable match within the errors. In this case, the most significant discrepancy is observed for NGC 5927, for which Ferraro et al. (2018) derived A rot = 2 . 3 km / s, while we find Arot = 0 . 76 + 0 . 90 -0 . 54 km / s. For this cluster also Sollima et al. (2019) derived a low probability of rotation.  \nIn the bottom row of Figure C.3 we compare the results of this work with those recently obtained by Petralia et al. (2024) by using APOGEE spectra for a sample of Galactic GCs. For A rot we use the semi-amplitude of the A fit values reported in their work. A reasonable overall agreement is observed also in this case for the clusters in common, however 47 Tuc and NGC 5904 result to have larger central velocity dispersion values and peak of rotation than in our work.  \nFinally, as shown in Figure C.4 (left panel) a reasonably good match is also found with the σ 0 estimates by Baumgardt  \n<!-- image -->  \nFig. C.4. Left panel: comparison with the best-fit σ 0 values obtained by Baumgardt et al. (2019) - B19 - for the total population of GCs in common with the present work. Right-panel: one-to-one comparison of the 3D rotation amplitude values A 3D found by Sollima et al. (2019) S19 - and the present analysis.  \n<!-- image -->  \net al. (2019). Among the clusters in common with Sollima et al. (2019), the only significantly discrepant result is that of 47 Tuc, which results to have a ∼ 25% larger rotation in this work. However, we note that in this case, as for the entire sample, both the i and θ 0 values are in very good agreement. In this respect, it is also interesting to highlight the nice match in terms of both the observed rotation amplitude and angles of the 3D rotation of 47 Tuc obtained in this work and those inferred by means of a detailed comparison between HST PMs and theoretical models of rotating clusters by Bellini et al. (2017).'}