Internal Dynamics of Multiple Populations in 28 Galactic Globular Clusters: A Wide-Field study with Gaia and the Hubble Space Telescope
Giacomo Cordoni Luca Casagrande Antonino Milone Emanuele Dondoglio Alessandra Mastrobuono-Battisti Sohee Jang Anna Marino Edoardo Lagioia Maria Vittoria Legnardi Tuila Ziliotto Fabrizio Muratore Vernica Mehta Elena Lacchin Marco Tailo


firstpage–lastpage
A generalized adaptive central-upwind scheme for compressible flow simulations and preventing spurious vortices
    [
    
===============================================================================================================




§ ABSTRACT 
    We present the first comprehensive analysis of the internal dynamics of multiple stellar populations (MPs) in 28 Galactic Globular Clusters (GCs) across a wide field of view, extending from the innermost regions to the clusters' outskirts. Using astro-photometric catalogs from ground-based observations, Gaia, and the Hubble Space Telescope (HST), we identify first- (1P) and second-population (2P) stars, and study the internal dynamics of MPs using high-precision Gaia DR3 and HST proper motions. Our results reveal that while the 1P transitions from isotropy to slight tangential anisotropy toward the outer regions, 2P stars become increasingly radially anisotropic beyond the half-light radius. We also explore the connection between the dynamics of MPs and the clusters' structural and dynamical properties, finding statistically significant differences in the anisotropy profiles of dynamically young and non-relaxed clusters, particularly beyond the 1-2 half-light radii. In these regions, 1P stars transition from isotropic to slightly tangentially anisotropic motion, while 2P stars become more radially anisotropic.
    In contrast, dynamically older clusters, with mixed MPs, exhibit weaker relative differences. Furthermore, clusters with orbits closer to the Galactic center exhibit larger dynamical differences between 1P and 2P stars than those with larger peri-Galactic radii. These findings are consistent with a scenario where 2P stars form in a more centrally concentrated environment, where the interaction with the Milky Way tidal field plays a crucial role in the dynamical evolution of MPs, especially of 1P.





stars: Hertzsprung-Russell and colour-magnitude diagrams; stars: kinematics and dynamics; Galaxy: globular clusters: general; Galaxy: kinematics and dynamics








§ INTRODUCTION




More than two decades of Hubble Space Telescope (HST) observations have revealed that nearly all Galactic Globular clusters (GCs) host Multiple Populations (MPs) with at least two main groups of stars: first-population (1P) and second-population <cit.>. These groups exhibit specific chemical compositions, with 2P stars being for example enriched in Sodium, Nitrogen and Helium  and depleted in Carbon and Oxygen <cit.>. Despite several models have been proposed to explain the formation of MPs, none can simultaneous match the numerous observational constraints <cit.>. Some scenarios predict that 2P stars formed from the ejecta of short-lived massive 1P stars <cit.>. As the expelled gas accumulated toward the cluster center, 2P stars would dominate the composition of the central cluster regions. Alternatively, the chemical composition of 2G stars could be the result of the accretion of chemically enriched material emitted from 1P stars in binary configuration <cit.>.

While some clusters indeed exhibit a more centrally concentrated 2P, there are numerous examples where 1P and 2P stars are spatially mixed, and even clusters with a more centrally concentrated 1P are observed <cit.>. Recent studies suggest that the extent of this mixing between 1P and 2P stars may correlate with the dynamical age of the cluster, indicating that dynamically older clusters tend to show a more complete mixing of MPs. Nonetheless, despite the numerous observational constraints, a comprehensive picture of how MPs interact form and evolve over time remains elusive.

In recent years, a promising new path has emerged: the study of the internal kinematics of cluster stars. This approach is crucial for uncovering the physical processes behind the formation of multiple populations. Specifically, N-body simulations indicate that the dynamical evolution of 2P stars formed in different environment should differ significantly from that of more spatially extended 1P stars. Such differences may still be detectable in present-day GC kinematics <cit.>, provided the populations are not completely relaxed. 



Over the past decade, a number of studies have analyzed the spatial distribution and the internal kinematics of MPs in various GCs, using HST photometry and proper motions <cit.>, MUSE/APOGEE line-of-sight (LoS) velocities <cit.>, ground-based photometry coupled with Gaia proper motions, or a combination of these methods <cit.>.
Such works presented the first evidence of significant dynamical differences between 1P and 2P stars in some GCs, qualitatively aligning with N-body simulations' results. However, the limited number of clusters and the small field of view prevented a complete and detailed characterization of the phenomenon.

In this work we aim to overcome the limitation of previous studies by investigating the internal dynamics of MPs in a sample of 28 GCs, combining photometric information from HST <cit.>, ground-based <cit.> and Gaia <cit.> and accurate HST <cit.> and Gaia <cit.> proper motions. This work presents the first comprehensive dynamical investigation of MPs in a large sample of clusters, extending from the central arcminutes to nearly the tidal radius.  The paper is structured as follows: in Sec. <ref> we introduce the dataset and MPs selection process. Section <ref> presents the study of the internal dynamics in the analyzed clusters and the global dynamical profiles of MPs. Finally, the discussion and summary are presented in Sec. <ref> and  <ref>, respectively. 




§ DATA


To investigate the internal dynamics of MPs in Galactic GCs over a wide field of view, we combined HST photometry <cit.>, ground-based UBVI photometry <cit.>, Gaia XP spectro-photometry <cit.>[Gaia XP spectra could be effectively used in only 4 GCs, namely NGC 0104, NGC 3201, NGC 6121 and NGC 6752.], with HST <cit.> and Gaia DR3 proper motions <cit.>. We refer to the different works for a detailed description of the astro-photometric datasets. The total sample used in this work consists of 28 GCs for which we could reliably select MPs stars over a wide field of view. 


To identify 1P and 2P among RGB stars, we have made use of the Chromosome Maps (ChMs) presented in <cit.> for HST data, photometric diagrams made with the  index[defined as = (U-B) -(B-I)<cit.>], <cit.> for ground-based observations, and the ChMs introduced in <cit.> from Gaia spectro-photometry. ChMs effectively and accurately separate MPs in GCs, with 1P stars defining a clump of stars centered around (0, 0), and 2P stars spread over larger absolute values of Δ_BI and  Δ_C_UBI (Δ_F275W, F814W and  Δ_C_F275W, F336W, F438W in HST filters).  In the next section, we discuss in detail the procedure used to identify 1P and 2P stars in the ground-based ChMs of 28 GCs, whereas we refer to <cit.> and <cit.> for the identification of MPs in HST and Gaia XP photometry. 

The final selection of 1P and 2P stars for each cluster and the analyzed field of view are shown in App. <ref> as online supplementary material. [The sample of analyzed clusters also includes type-II GCs, i.e. clusters with internal variation of heavy-elements <cit.>. However, considering the small number of stars belonging to the anomalous populations, i.e. enriched stars, we only study stars with light-elements abundance variations.]



 §.§ Multiple Populations from ground-based photometry


In a nutshell, the I vs. B-I and I vs. C_UBI=U - 2B + I CMDs of Red Giant Branch stars (RGB) have been used to derive the verticalized Δ_BI and  Δ_C_UBI colors, that are respectively shown on the x and y axis of the ChMs. We refer to <cit.> for a detailed discussion on how these quantities have been derived. 

To select 1P and 2P stars from the ChMs we adapted the procedure introduced in <cit.> on HST ChMs to our ground-based ChMs, as done in Jang et al., in preparation. Figure <ref> illustrates the procedure for the cluster NGC 2808. Briefly, we first selected by eye the bulk of 1P stars among the stars with (Δ_BI, Δ_C_UBI) ∼ (0, 0), and fitted a straight line to 1P stars (black solid line in Fig. <ref>, leftmost panel). We then rotated the ChMs by the angle θ defined by slope of the 1P best-fit line. We refer to the new rotated reference frame as  and  (Fig. <ref>, middle panel). By construction, 1P stars define a horizontal distribution in  vs , and we can use the  distribution to disentangle 1P and 2P stars. We used the  python package <cit.> and employed Gaussian Mixture Models (GMM) to fit 2 Gaussians to the  distribution. We remind here that GMM is applied directly to the data points, so that the results do not depend on the bin choice.
Finally, 1P stars (red points in the rightmost panel) have been selected as stars with < ^sep, with ^sep being the intersection between the two best-fit Gaussians, indicated by the horizontal line in the middle panel of Fig. <ref>. 2P stars are marked with blue points. Additionally, stars with unusual values of  and  have been excluded from the selection.

The ground-based ChMs with selected 1P and 2P stars of the 28 analyzed Galactic GCs are shown in App. <ref> as online supplementary material.

Additionally, we employed the V vs.  pseudo-CMD to extend the regions analyzed in <cit.>, including stars at larger distance from the cluster centers. We refer to <cit.> for a detailed description of the selection of MPs.
 





 §.§ Multiple populations from HST photometry

To study MPs among RGB stars in the innermost cluster regions we exploited the ChMs and selection carried out in <cit.>. In a nutshell, HST ChMs have been derived by means of the verticalized m_F814W vs. m_F275W - m_F814W (for the abscissa) and m_F814W vs. C_F275W, F336W, F438W[C_F275W, F336W, F438W = m_F275W - 2 m_F336W + m_F438W] CMDs (on the ordinate). 1P and 2P stars have been selected with the procedures detailed in <cit.>. Each HST ChM has then been cross matched with the proper motion catalogs of <cit.>.



 §.§ Multiple populations from Gaia XP synthetic photometry

The recent Gaia DR3, <cit.> publicly released low-resolution spectra, namely XP spectra, for ∼200 million sources <cit.>, allowing to derive synthethic photometry in virtually any photometric systems <cit.>. In a recent work, <cit.> exploited Gaia XP spectra to compute synthetic photometry special filters designed to maximize the separation between 1P and 2P stars in GCs' RGB. Their analyzed field of view extends well beyond that covered by the ground-based photometric dataset of <cit.>, reaching their tidal radius. We refer to <cit.> for a detailed description of the procedure used to derive and separate the populations, and for the validation of such identification. In a nutshell, exploiting synthethic spectra with the typical composition of 1P and 2P stars, Mehta and collaborators defined a series of photometric filters to maximize the separation between 1P and 2P stars in 5 Galactic GCs. To validate their MPs identification, <cit.> compared the classification with available spectroscopic information, finding consistent results. In this work, we exploit their identification to extend our ground-based dataset and investigate the dynamics of 1P and 2P stars up to the clusters' tidal radii.




§ INTERNAL DYNAMICS OF MULTIPLE POPULATIONS


After carefully selecting MPs in the 28 analyzed clusters, we exploited HST <cit.> and Gaia DR3 <cit.> proper motions to investigate the internal dynamics of stars belonging to different populations. Specifically, we followed the procedure described in <cit.> and <cit.>[see also <cit.>.]. We first transformed the celestial coordinates and their proper motion components (α, δ, , )[=μ_αcosδ] into a Cartesian reference frame (x, y, , ), using the orthographic projection <cit.>, and
then the proper motions into their sky-projected radial and tangential components, defined as = x  + y and = y  - x, with  pointing outward (positive) or inward (negative), and  positive when counterclockwise.  and  here indicate the proper motion of each star relative to the motion of the cluster, determined by <cit.>. Furthermore, considering that the sample of analyzed stars is different from that analyzed in <cit.>, we repeated the analysis determining the mean cluster motion from our sample of stars, and accounting for systematic errors as in <cit.>[routines available at <https://github.com/GalacticDynamics-Oxford/GaiaTools>.], finding consistent results.

The uncertainties on the radial and tangential components of the proper motions have been determined from the uncertainties on  and , accounting for the correlation between the proper motions. Finally, the radial proper motions have been corrected for the perspective expansion/contraction due the bulk motion of the cluster along the line of sight, as in <cit.>. 

The mean motion and velocity dispersion of 1P and 2P stars have then been computed by dividing the field of view into different concentric annuli containing approximately the same number of stars. For each annulus we determined the mean motion (μ_R/T) and dispersion (σ_R/T) from Gaia DR3 radial and tangential components by minimizing the negative log-likelihood defined in <cit.> including the covariance term as in <cit.>. Uncertainties on the determination of the mean motion and dispersion have been determined using Markov Chain Monte Carlo algorithm <cit.>. 

The dynamical profiles have been re-computed using different number of bins and  bins with the same radial width to determine the effect of arbitrary bin choice, with consistent results. 
Concerning the HST proper motions, we used the same likelihood without the covariance term, <cit.> and maximizing only for the radial/tangential dispersion. This is because HST proper motions are relative to the cluster's bulk motion and do not provide individual stellar absolute mean motion.

In the following sections, we present and discuss the results for selected individual clusters as well as the global dynamical profiles derived by combining data from all clusters. This approach enables us to explore the general dynamical properties of 1P and 2P stars. Our analysis focuses on the dispersion profiles, while we defer a detailed investigation of the rotation of 1P and 2P stars to a forthcoming paper that is currently in preparation.



 §.§ Individual dynamical profiles
 

In this section, we present the study of the dynamical profiles of 1P and 2P stars in a few individual clusters with sufficiently large samples, namely NGC 0104, NGC 2808, NGC 5904, and NGC 6205. The dynamical profile of the remaining clusters are shown in App. <ref> as online supplementary material.


Figure <ref> shows, from top to bottom, the mean radial and tangential profiles, the radial and tangential dispersion, and the anisotropy profiles, defined as β=. The red and blue solid lines indicate the dynamical quantities, determined as described in Sec. <ref>, while the shaded rectangles indicate the relative uncertainties and the extension of each radial bin. The distance from the cluster center are in units of  indicated in the middle panel <cit.>. Black diamonds indicate the values of the 1-dimensional dispersion from <cit.>, while the regions with available HST data is indicated by the vertical green lines. The top and right axis (shown in gray) of each panel display distances in units of parsec and km/s, converted considering cluster distances determined in <cit.>.

A visual inspection of the individual profiles reveals that, as expected, the mean radial motions of cluster stars are consistent with 0, as we already subtracted the perspective contraction/expansion. Consistently with <cit.>, we find that all four clusters exhibit rotation in the plane of the sky. The mean tangential profiles of NGC 0104 and NGC 5904 exhibit clear negative and positive values, respectively indicating counterclockwise and clockwise rotation. Possible differences in the rotation profile are consistent with the inferred uncertainties. NGC 2808 and NGC 6205, show possible hints of stellar rotation, even though the signal too weak to carry out a comparison. 



Concerning the dispersion profiles, 1P stars exhibit a larger tangential and similar radial dispersion compared with 2P stars, especially outside ∼ 1-2, while 1P and 2P have similar radial dispersion profiles. As a result, 1P stars exhibit an isotropic motion in the central regions, possibly shifting toward tangentially anisotropic motion in the outer regions. On the other hand, the motion of 2P stars is consistent with being isotropic in the cluster center, becoming radially anisotropic beyond 1. A similar pattern is also observed in NGC 3201, consistently with <cit.>. The mean, dispersion and anisotropy profiles of the remaining clusters are displayed in App. <ref>. Overall, we find a good agreement between HST and Gaia dynamical profiles in overlapping or neighbouring regions.

The observed dynamical profiles of 1P and 2P stars are qualitatively consistent with the predictions from N-body simulations where 2P stars form in a more centrally concentrated environment. We refer to <cit.> for a detailed description of the assumptions adopted in the simulations and the results. Additional discussion is carried out in Sec. <ref>.




 §.§ Global dynamical profiles

Given the low number of 1P and 2P stars in many clusters, a detailed analysis of the individual dynamical profiles as a function of cluster radius is often unfeasible. Therefore, similarly to <cit.>, we opted to investigate the internal dynamics in all clusters simultaneously, combining together the dispersion profiles of the individual clusters[We find worth mentioning that <cit.> normalized positions and proper motions of cluster stars and then derived the dispersion profiles. However, due to the non-relative nature of Gaia proper motions, we cannot use the  same procedure on our dataset. Hence, we first derive the dynamical profiles in each cluster, and then normalize and combine the profiles.]. In order to properly compare different clusters, we normalized the radial coordinates to the cluster  from <cit.> and the dispersion profiles to the central dispersion determined in <cit.> from HST data. Four clusters do not have determination of the central velocity dispersion from HST data, namely NGC 1904, NGC 4147, NGC 6712 and NGC 7492. For these clusters, we adopted the central dispersion determined in <cit.>. We remind here that such approach can be used only to compare dispersion profiles, but not mean motions, as the latter are not on a relative scale. 

The normalized global dispersion and anisotropy profiles as a function of distance from cluster center for the 28 analyzed GCs are shown in Fig. <ref>. Global profiles for 1P and 2P stars are shown in the first and second columns, while the comparison between the average trend is shown in the third column. The average profiles were determined using the Locally Weighted Regression (LOESS) as implemented in <cit.>. To estimate the uncertainties in the average trend, we repeated the LOESS fitting on 1000 bootstrapped samples, drawing the values of σ_R/T in each realization from a Gaussian distribution centered on the observed value with a dispersion equal to the observed uncertainties. The uncertainties were determined as the 16^th and 84^th percentiles of the 1000 LOESS fits. This is a more robust approach than simply bootstrapping the sample, as it accounts for observed uncertainties as well. Finally, we estimated the significance of the observed differences between 1P and 2P stars accounting for uncertainties in both data and fitting. The details of the procedure are outlined in App. <ref>. In a nutshell, for each of the 1000 realizations of the LOESS fit, we computed the difference between 1P and 2P as a function of the radial coordinate. The 1, 2 and 3σ confidence regions are shown as gray shaded regions in the rightmost panels, where we display the difference between 1P and 2P a function of normalized radius. Additionally, we quantified the statistical significance of each point as the fraction of simulations which returned a difference smaller than the observed one, at the same radial location. The color of the line in the rightmost panels is indicative of of the statistical significance, in units of σ, as indicated in the right colorbar.

While there are no differences between the global radial dispersion profiles of 1P and 2P stars, 2P stars exhibit a lower tangential dispersion with respect to 1P stars, especially in the outer regions, e.g. R > 1-2. As a consequence, the 1P displays an isotropic/tangentially anisotropic motion in the outer regions, while the 2P is radially anisotropic. Due to the large uncertainties of the mean trend, the observed differences are significant beyond the 3σ level only between 2 and 4 . Nonetheless, the magnitude of the differences is consistent with the predictions from the theoretical simulations <cit.>. Additionally, the global profiles within  are consistent with the results of <cit.>.





§ DISCUSSION

In the last decade, several theoretical works based on N-body simulations have investigated the dynamical evolution of MPs in GCs and their potential for distinguishing between different formation scenarios. The velocity dispersion and anisotropy profiles of MPs can provide fundamental insights into the origins of these stellar populations <cit.>. However, due to dynamical evolution and spatial mixing, possible differences in the internal dynamics of MPs disappear as the clusters reach relaxation. As shown in Fig.<ref> for individual clusters and Fig.<ref> for all clusters, the tangential dispersion and anisotropy profiles of 1P and 2P stars exhibit some differences, with 1P being isotropic across the entire cluster field and 2P stars shifting from isotropy to radially anisotropy beyond .

To investigate how these dynamical differences depend on the clusters' properties evolutionary state, we divided our sample of 28 clusters into different groups based on various internal and external factors: dynamical ages <cit.>, tidal filling factor defined as ℛ=/ <cit.> with  being the Jacobi radius derived in <cit.>[Clusters with ℛ < 0.05 are considered in “isolated regime” or tidally underfilling clusters, while clusters with ℛ > 0.05 are in “tidal regime”, or tidally filling.], 
spatial mixing of MPs <cit.>, tidal interaction with the Milky Way <cit.>[Consistently with literature works, we divide clusters with R_peri≶3.5 kpc], clusters' origin (e.g. in situ or accreted) as identified in <cit.> and escape velocity <cit.>[See e.g. the analysis of <cit.> and <cit.> and the connection between the cluster's  and MPs. The value of the limiting escape velocity, i.e. 20 km/s, is determined examining the distribution of escape velocities from <cit.>.]. The groups are shown in the table provided in App. <ref>. We find worth mentioning that different groups have overlap, as some clusters' properties are not independent one from the other <cit.>. Additionally, we also divided clusters based on their internal bulk rotation <cit.>, but, to keep the figure simpler, we not show the results in Fig. <ref>. Nonetheless, the difference between 1P and 2P anisotropy profiles for the groups and the statistical significance are included in App. <ref>.


The dynamical profiles for different groups of clusters are displayed in Fig. <ref>. Each group is indicated in the top inset, together with the number of clusters in each group. The global profiles (solid lines) have been computed with the LOESS algorithm, while the uncertainties (shaded regions, 1σ) have been determined by bootstrapping with replacements a 1000 times and accounting for the observational uncertainties (see Sec. <ref> for a detailed description). The significance of the differences between the anisotropy profiles of 1P and 2P for each group is presented in App. <ref>, with the description of the procedure adopted to estimate it.

We find significant differences, i.e. above the 3σ level, in the average anisotropy profiles of many of the analyzed groups. Overall, we find that the 2P is more radially anisotropic beyond , while 1P is isotropic and become slightly tangentially anisotropic beyond 2-3. There are important differences in the 1P-2P relative dynamics among some of the analyzed groups. For example, while dynamically young clusters and those with centrally concentrated 2P stars show clearly different dispersion and anisotropy profiles, these differences are less pronounced in intermediate-age and mixed-population clusters. These findings are qualitatively consistent with the simulations by <cit.> and <cit.>, which suggest that due to different initial spatial configurations (such as centrally concentrated 2P stars) or initial isotropy, 1P stars tend to develop isotropic motion, while 2P stars display more radially anisotropic motion, especially in the outer regions. Over time, these dynamical differences gradually diminish as spatial and dynamical mixing of MPs occurs.

The comparison between clusters with small and large peri-Galactic radii (shown in the first two columns of the bottom row of Fig. <ref>) reveals some interesting patterns. In clusters with small peri-Galactic radii, 1P stars exhibit tangential anisotropy beyond , while in clusters with larger peri-Galactic radii, 1P stars are either isotropic or slightly radially anisotropic. In both cases, 2P stars become increasingly radially anisotropic as they move outward. These findings suggest that the interaction with the Milky Way’s tidal field plays a significant role in the dynamical evolution of MPs. Moreover, in clusters with small peri-Galactic radii, 1P stars appear to behave like tidally filling systems, where the outer regions develop slight tangential anisotropy, whereas 2P stars resemble tidally underfilling systems with radial anisotropy. In contrast, weaker interactions with the Milky Way’s tidal field result in a more isotropic 1P population. This scenario aligns qualitatively with the theoretical models presented in <cit.>.

To further explore this hypothesis, we also examined the internal dynamics of MPs in clusters classified as tidally underfilling (ℛ < 0.05) and tidally filling (ℛ > 0.05). While tidally underfilling clusters exhibit relative differences consistent with other groups, such as isotropic 1P and radially anisotropic 2P, MPs in tidally filling clusters display a different pattern. Specifically, 1P stars are tangentially anisotropic, while 2P stars exhibit isotropic motion across the entire field. Interpreting these results is challenging, as the tidal filling factor ℛ is closely linked to other properties, such as , R_peri, and cluster mass. Specifically, tidally underfilling clusters tend to be more massive, with larger  and smaller R_peri.

A qualitatively similar pattern is observed in clusters with low and high escape velocities, with the latter often exhibiting tangentially anisotropic 1P stars. This difference may be due to the preferential loss of 1P stars on radial orbits in clusters with higher  <cit.>. However, the relationship between  and the fraction of MPs remains complex and is highly dependent on the specific formation scenario adopted for MPs.

1P and 2P stars exhibit comparable relative differences in clusters, regardless of whether they are in situ or accreted, and whether they are rotating or non-rotating (see, e.g. App. <ref>). In these cases, 2P stars remain radially anisotropic, while 1P stars display isotropic motion across the entire cluster field.

In summary, the results presented in Fig. <ref> indicate that clusters in a less dynamically evolved state show significant dynamical differences among MPs. Additionally, the interaction with the Galaxy appears to play a crucial role in shaping the evolution of different populations. These findings align with the conclusions of <cit.> regarding the internal dynamics of MPs in the innermost regions.







§ SUMMARY AND CONCLUSION


In this study, we present the first comprehensive investigation of the internal dynamics of multiple populations across a wide field of view, from the innermost arcmin to the clusters' outskirts, in a large sample of 28 Galactic GCs. Using HST, ground-based and XP synthethic photometry, we identified first- and second-population stars among RGB stars from the innermost regions to the clusters' outskirts. To achieve this, we exploited the pseudo CMD dubbed Chromosome maps as well as the  index, where 1P and 2P stars form well-separated groups. The internal dynamics of MPs was investigated using HST and Gaia DR3 proper motions, allowing us to determine the mean, dispersion, and anisotropy profiles as function of distance from the cluster center. The analyzed cluster regions range from 0.2 to more than 10 . The analysis presented in this work confirms and extends the results of previous studies based on HST and Gaia proper motions <cit.>. 

Our analysis of individual clusters reveals distinct dynamical differences in the anisotropy profiles of 1P and 2P stars in several cases. On average, 2P stars are more radially anisotropic beyond , whereas 1P stars generally exhibit isotropic motion, with some showing tangential anisotropy beyond 2-3. These anisotropy differences between 1P and 2P stars are primarily driven by a lower tangential dispersion of 2P, with no significant differences observed in the radial component. These results agree with the analysis of <cit.> for the innermost cluster regions.

Studying the global trends in dispersion and anisotropy profiles, derived by combining the dynamical profiles of various clusters, we observe significant differences between 1P and 2P stars, especially outside  where the differences are significant beyond the 3σ level. 
In dynamically young clusters, 1P stars are isotropic in the inner regions but become slightly tangentially anisotropic toward the outskirts. In contrast, 2P stars are isotropic in the cluster centers and become radially anisotropic in the outer regions. These patterns are also evident in clusters where 2P stars are more centrally concentrated and in clusters with escape velocities exceeding 20 km/s. However, dynamical evolved clusters, with spatially mixed MPs and lower escape velocities do not exhibit these dynamical differences. Rotating/non-rotating clusters, and in situ/accreted clusters show similar relative differences between 1P and 2P stars.

We also explore the influence of the Milky Way’s tidal field on the dynamical properties of MPs by analyzing clusters with different peri-Galactic radii. Our findings reveal a strong connection between the dynamical behavior of 1P stars and the strength of the Milky Way’s tidal field. In clusters with orbits closer to the Galactic center, where the tidal field is stronger, 1P stars tend to exhibit tangential anisotropy beyond 1-2. Conversely, in clusters with weaker interactions with the Milky Way’s tidal field, 1P stars display isotropic motion. This suggests that the Milky Way’s tidal field plays a crucial role in the dynamical evolution of MPs. Further supporting this conclusion, our analysis shows that tidally underfilling and filling clusters exhibit distinct relative patterns in the dynamical profiles of 1P and 2P stars.

The observed differences in the internal dynamics between 1P and 2P stars qualitatively align with the predictions of N-body and theoretical simulations by <cit.>. These studies indicate that the distinct dynamical properties of 1P and 2P stars are indicative of 2P stars forming in a more centrally concentrated environment. 


Overall, the analysis presented in this paper offers key insights into the formation scenarios of multiple stellar populations and their relationship with both internal factors (such as  and dynamical age) and external influences (such as interaction with the Milky Way’s tidal field).




§ ACKNOWLEDGMENTS


SJ acknowledges support from the NRF of Korea (2022R1A2C3002992, 2022R1A6A1A03053472). EPL acknowledges support from the “Science
& Technology Champion Project” (202005AB160002)
and from the “Top Team Project” (202305AT350002),
all funded by the “Yunnan Revitalization Talent Support Program”. TZ acknowledges funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie
Grant Agreement No. 101034319 and from the European Union – Next Generation EU. This work has received funding from  “PRIN 2022 2022MMEB9W - Understanding the formation of globular clusters with their multiple stellar generations” (PI Anna F. Marino),  and from INAF Research GTO-Grant Normal RSN2-1.05.12.05.10 -  (ref. Anna F. Marino) of the “Bando INAF per il Finanziamento della Ricerca Fondamentale 2022”. EL acknowledges financial support from the European Research Council for the ERC Consolidator grant DEMOBLACK, under contract no. 770017






§ DATA AVAILABILITY


Relevant data underlying this work is available in the article. All other data will be shared upon reasonable request to the corresponding author.






mnras
 




















§ SELECTION OF MULTIPLE POPULATIONS

In this Appendix we show the selection of MPs for all the 28 analyzed clusters. As discussed in Sec. <ref>, we used a combination of HST and ground-based photometry <cit.>, and Gaia XP synthetic photometry <cit.>. In Fig. <ref>- <ref> we display the analyzed field of view in the central panel, together with the photometric diagrams used to separate 1P and 2P, always shown with red and blue colors respectively. Specifically, ground-based photometry is shown in the leftmost panels, with the ChMs from <cit.> shown in the top one. HST and Gaia XP ChM are instead displayed in the right panels, respectively in the top and bottom one. The coloured circles in the central panels indicate the analyzed regions. 

Concerning ground-based observations, we first used the selection in the ChM whenever available, and we adopted the  selection for stars without ChM information. MPs outside the orange circles are selected by means of Gaia XP synthethic photometry. We remind here that Gaia XP synthethic photometry is only available for four clusters, namely NGC 0104, NGC 3201, NGC 6121 and NGC 6752 <cit.>.











§ INDIVIDUAL PROFILES OF ALL CLUSTERS


We present in this section the individual dynamical profiles of the analyzed clusters. The profiles are displayed in a similar fashion as Fig. <ref>. Mean trends are shown in Fig. <ref>-<ref>, while dispersion and anisotropy profiles are presented in Fig. <ref>-<ref>. Values computed from HST/Gaia data are displayed as crosses and circles respectively, while their uncertainties are indicated by the shaded rectangles. Mean motions and dispersion profiles have been converted to km/s adopting the clusters' distances derived in <cit.>. Red and blue colors indicate 1P and 2P stars, while black diamonds represent the values of the 1D-dispersion computed in <cit.>. Finally, the locations of the , , radii are indicated by the dash-dotted, solid and dashed lines respectively.  is indicated only if contained within the analyzed field of view.













§ STATISTICAL SIGNIFICANCE OF THE OBSERVED GLOBAL PROFILES.
 
To assess the statistical significance of the differences in the global dynamical profiles of 1P and 2P shown in Fig. <ref> and  <ref>, we adopted the following procedure. We refer for simplicity to the case with all clusters. First we created a 1000 realizations of bootstrapped samples of y=, , β, where each value has bee drawn from a Gaussian with dispersion equal to the observed uncertainty. For each realization we repeated the LOESS fit for 1P and 2P, and computed the difference the between the two fits. Finally, we estimate the 1, 2, 3σ confidence intervals from the distribution of simulated difference at each R/. The three confidence intervals are shown as gray shaded regions in the top panels of Fig. <ref>. In addition to this analysis, we also directly quantified the fluctuations introduced solely by the uncertainties. To do this, we repeated the same analysis on two identical distributions, testing the null hypothesis that 1P and 2P stars share the same dynamical profiles, and determined the fraction (f) of simulations with differences larger than the observed ones. Such fraction indicates the probability that uncertainties alone can reproduce a difference between 1P and 2P as large as the observed ones, and thus we computed the significance of the difference as p=1-f. The value of p is indicated by the color of the lines in the top panels of Fig. <ref>, as indicated by the colorbar.
We stress here that, instead of determining one single value of significance for each 1P/2P dynamical profiles, we compute the statistical significance for each radial coordinate.










