Dynamical Evolution of Multiple-Population Globular Clusters
MMNRAS , 1–16 (2020) Preprint 4 February 2021 Compiled using MNRAS L A TEX style file v3.0
Dynamical Evolution of Multiple-Population Globular Clusters
Enrico Vesperini, (cid:63) Jongsuk Hong, , Mirek Giersz, Arkadiusz Hypki Department of Astronomy, Indiana University, Bloomington, Swain West, 727 E. 3rd Street, IN, 47405, USA Department of Astronomy, Yonsei University 50 Yonsei-Ro, Seodaemun-Gu, Seoul 03722, Republic of Korea Korea Astronomy and Space Science Institute, Daejeon 34055, Republic of Korea Nicolaus Copernicus Astronomical Centre, Polish Academy of Sciences, ul. Bartycka 18, 00-716 Warsaw, Poland
ABSTRACT
We have carried out a set of Monte Carlo simulations to study a number of fundamentalaspects of the dynamical evolution of multiple stellar populations in globular clusterswith different initial masses, fractions of second generation (2G) stars, and structuralproperties. Our simulations explore and elucidate: 1) the role of early and long-termdynamical processes and stellar escape in the evolution of the fraction of 2G starsand the link between the evolution of the fraction of 2G stars and various dynamicalparameters; 2) the link between the fraction of 2G stars inside the cluster and in thepopulation of escaping stars during a cluster’s dynamical evolution; 3) the dynamicsof the spatial mixing of the first-generation (1G) and 2G stars and the details of thestructural properties of the two populations as they evolve toward mixing; 4) theimplications of the initial differences between the spatial distribution of 1G and 2Gstars for the evolution of the anisotropy in the velocity distribution and the expectedradial profile of the 1G and 2G anisotropy for clusters at different stages of theirdynamical history; 5) the variation of the degree of energy equipartition of the 1Gand the 2G populations as a function of the distance from the cluster’s centre and thecluster’s evolutionary phase.
Key words: globular clusters:general – stars:kinematics and dynamics –stars:chemically peculiar
Globular clusters host multiple stellar populations charac-terized by a variety of patterns in the abundances of a num-ber of light elements (e.g. C, N, Na, O, Mg, Al); a few clus-ters also show further complexities in the chemical propertiesof their populations including variations in the abundancesof Fe and s-process elements (see e.g. Gratton et al. 2019 fora recent review and references therein). The multiple stellarpopulations found in many clusters differ not only in theirchemical abundances but also in their structural and kine-matic properties (see e.g. Sollima et al. 2007, Bellini et al.2009, Lardo et al. 2011, Simioni et al. 2016, Dalessandro etal. 2019, Richer et al. 2013, Bellini et al. 2015, 2018 Corderoet al. 2017, Lee 2017, Cordoni et al. 2020a, 2020b); these dif-ferences clearly reveal the existence of a tight link betweenthe origin of the populations’ chemical properties and thecluster’s formation and dynamical history and add anotherfundamental node in the complex network of relationships (cid:63)
E-mail: [email protected] between dynamics and the properties of the stellar contentof globular clusters.The numerous observational and theoretical studies ofmultiple stellar populations are revealing an increasinglycomplex picture of the chemical and dynamical properties ofglobular clusters and raising many new questions. The mainissues in the study of multiple populations in globular clus-ters concern the origin of the processed gas from which starswith anomalous chemical properties form, the formation his-tory of multiple populations, the differences in their initialstructural and kinematic properties, how these differencesaffected the cluster dynamical evolution, and the extent towhich the current properties retain memory of these initialdifferences.As for the gas out of which second generation starsformed (hereafter 2G; notice that according to some of themodels listed below different populations would form at thesame time and in that case it is more appropriate to usethe terms first- and second-population to denote stars withdifferent chemical properties) a variety of different stellarsources have been proposed including AGB stars, rapidly-rotating massive stars and massive binaries, supermassive © a r X i v : . [ a s t r o - ph . GA ] F e b E. Vesperini et al. stars, black hole accretion disks, stellar mergers (see e.g.Ventura et al. 2001, Decressin et al. 2007, de Mink et al.2009, Bastian et al. 2013, Krause et al. 2013, D’Ercole et al.2008, 2010, 2012, Bekki 2010, Bekki et al. 2017, Denissenkov& Hartwick 2014, D’Antona et al. 2016, Elmegreen 2017,Kim & Lee 2018, Breen 2018, Gieles et al. 2018, Calura etal. 2019, Howard et al. 2019, Wang et al. 2020, McKenzie &Bekki 2021). No consensus has been reached on this funda-mental aspect of the study of multiple populations and manyof these models still require significant further developmentconcerning the nucleosynthesis of the proposed polluters andthe comparison with all the observational constraints (seeGratton et al. 2019).As pointed out above, different stellar populations havebeen found in many cases to be characterized by differentdynamical properties and the study of the dynamical impli-cations of these differences is essential to build a completepicture of globular clusters’ formation and evolution.In all models dynamics plays a key role in determiningthe initial differences in the structural and kinematic proper-ties of multiple populations and their subsequent evolution(see e.g. D’Ercole et al. 2008, Bekki 2010, 2011, Vesperiniet al. 2013, Calura et al. 2019), and, in some of the mod-els presented in the literature, the cluster’s dynamics playsa central role also in the formation of the stellar pollutersproviding gas for the formation of 2G stars (see e.g. Gieleset al. 2018, Wang et al. 2020).The first studies of multiple-population cluster dynam-ics have addressed a number of aspects related to the hy-drodynamics of the gas out of which 2G stars formed (seee.g. D’Ercole et al. 2008, 2016, Bekki 2010, 2011, Calura etal. 2019). These studies have shown that 2G stars tend toform spatially concentrated in the cluster’s inner regions;other models, although not based on detailed hydrodynami-cal simulations, also agree and provide qualitative considera-tions supporting this expected difference between the initialspatial distributions of 2G and 1G stars. A few previousstudies (see e.g. Vesperini et al. 2013, 2018, Miholics et al.2015, Henault-Brunet et al. 2015, Fare et al. 2018, Tiongcoet al. 2019) have shown that initial differences in the spatialdistribution of 1G and 2G stars are gradually erased duringa cluster’s evolution. Dynamically younger clusters may re-tain some memory of the initial differences until the presentday while in clusters that have reached the more advancedstages of their evolution the two populations are expected tobe completely mixed. Examples of both of these cases havebeen found in observational studies (see e.g. Sollima et al.2007, Bellini et al. 2009, Lardo et al. 2011, Beccari et al.2013, Cordero et al. 2014, Simioni et al. 2016, Boberg et al.2016, Lee 2017, Gerber et al. 2020b for studies showing clus-ters in which the 2G is more concentrated than the 1G andDalessandro et al. 2014, Nardiello et al. 2015, Cordero etal. 2015, Gerber et al. 2018, 2020a for clusters in which thetwo populations are instead completely mixed) and a recentanalysis by Dalessandro et al. (2019) has provided an ob-servational picture of the evolutionary path toward spatialmixing.Dynamical differences between the 2G and the 1G popu-lations are not limited to their structural properties: a num-ber of studies have shown that, either as a result of theformation process or the subsequent dynamical evolution,the kinematic properties of the 2G and 1G stars may differ: the 2G population may be characterized by a more rapid ro-tation and a more radially anisotropic velocity distributionthan 1G stars (see e.g. Bekki 2010, Mastrobuono-Battisti& Perets 2013, 2016, Bellini et al. 2015, Henault-Brunet etal. 2015, Tiongco et al. 2019). The first observational stud-ies of the kinematics of multiple populations have revealedsuch differences in a few clusters (see e.g. Richer et al. 2013,Bellini et al. 2015, 2018, Cordero et al. 2017, Milone et al.2018, Libralato et al. 2019, Cordoni et al. 2020a, 2020b; seealso e.g. Pancino et al. 2007, Cordoni et al. 2020b, for clus-ters in which the multiple populations share similar kine-matic properties). Significant additional efforts beyond theseinitial pioneering studies will be necessary to draw a morecomplete observational picture of the dynamical propertiesof multiple-population clusters. In particular, it will be nec-essary to extend the observational studies to cover a widerradial range of distances from clusters’ centres to include theouter regions where structural and kinematic differences, ifnot completely erased during a cluster’s long-term dynami-cal evolution, should be stronger.The differences between the 1G and the 2G structuralproperties have also been shown to have important implica-tions for the evolution and survival of binary stars. For themore centrally concentrated 2G population, binaries may bedisrupted and evolve more rapidly than 1G binaries leadingto different 1G and 2G binary fractions and differences be-tween the orbital properties of the surviving binaries (seee.g. Vesperini et al. 2011, Hong et al. 2015, 2016, 2019, Lu-catello et al. 2015, Dalessandro et al. 2018, Gratton et al.2019, Milone et al. 2020, Kamann et al. 2020 for some theo-retical and observational studies on the properties of binariesin multiple populations). The possible role of the differencesbetween the formation history of 1G and 2G populations onthe formation of binaries and their initial properties is stillunexplored and is certainly another important aspect in thestudy of binaries in multiple-population clusters.In this paper, we present the results of a set Monte Carlosimulations following the evolution of multiple-populationglobular clusters. We explore the role played by early andlong-term evolutionary processes in the dynamical evolutionof multiple populations, in driving the escape from the clus-ter of 1G and 2G stars, and determining the evolution ofthe fraction of 2G stars in the cluster. We follow the evolu-tion of the main structural and kinematic properties duringthe various stages of a cluster’s evolution and discuss whatfingerprints of the formation process may still be observabletoday along with the dynamical signatures and implicationsof the initial differences between the structural properties ofthe 1G and the 2G populations.This is the structure of the paper: in Section 2 we de-scribe the method and the initial conditions used for oursimulations. In section 3 we present our results, and we sum-marize our conclusions in section 4.
This study is based on a set of Monte Carlo simulationsrun with the MOCCA code (Hypki & Giersz 2013, Gierszet al. 2013). The code includes the effects stellar and binaryevolution modeled using the SSE and BSE codes (Hurleyet al. 2000,2002) and assuming supernovae kick velocities
MNRAS000
MNRAS000 , 1–16 (2020) ynamics of multiple populations folllowing a Mawellian distribution with a dispersion equalto 265 km/s (Hobbs et al. 2005), two-body relaxation andbinary star interactions, and a truncation radius mimickingthe effect of the tidal truncation due to the external tidalfield of the host galaxy for a cluster on a circular orbit (seeHypki & Giersz 2013, Giersz et al. 2013 for further detailson the MOCCA code).Table 1 summarizes the initial conditions of all our sim-ulations along with the id used throughout the paper to referto each simulation.The systems we have considered for our study start witha total number of stars equal to N = or N = × withno primordial binaries and with masses distributed accord-ing to a Kroupa (2001) initial mass function (IMF) between0.1 m (cid:12) and 100 m (cid:12) . In all the systems, the density profile of2G stars is modeled as a King model with W = (where W is the central dimensionless potential; King 1966) while the1G population is less concentrated and distributed accord-ing to a tidally limited King model with W = or W = . Inall cases the 2G population is confined in the inner regionsof the cluster; we have explored models with values of theinitial ratio of the 1G half-mass radius to the 2G half-massradius, r h , / r h , , equal to 20 and 10. As for the initial ra-tio of the total 2G mass to the total cluster mass, M / M tot ,we have studied systems with M / M tot =0.1, 0.25 and, fora few models, 0.4. The initial conditions adopted follow thegeneral trend emerging from the results of a number of hy-dro simulations of multiple population formation (see e.g.D’Ercole et al. 2008, Bekki 2010, 2011, Calura et al. 2019)showing that 2G stars form centrally concentrated in a moreextended 1G system. As for the initial 2G mass fraction, wehave explored a few different values since this quantity is stillpoorly constrained and may depend on a number of factorsincluding the amount of gas supplied by the 1G stellar pol-luters and the amount of pristine gas involved in the 2G starformation (see e.g. Calura et al. 2019).All the systems start with an isotropic velocity distri-bution but for the model sg025c20wf we have also exploredtwo initial conditions with an initially radially anisotropicvelocity distribution with anisotropy following an Osipkov-Merritt profile (see e.g. Binney & Tremaine 2008), β = − ( σ θ + σ φ ) / ( σ r ) = / ( + r a / r ) where σ θ , σ φ , and σ r arethe three components of the velocity dispersion in spheri-cal coordinates, and r a is an anisotropy scale radius beyondwhich the velocity distribution becomes increasingly rada-ially anisotropic. For the two anisotropic models we havestudied, the anisotropy radius, r a was set equal to r h / and r h / where r h is the cluster’s initial half-mass radius.Finally we have considered two different values for theinitial truncation radius corresponding to tidal fields withdifferent strength which we will refer to as strong field (sf)and weak field (wf). The two values of the truncation radiiadopted would correspond to the tidal radii of clusters withthe masses of our models and moving on circular orbits ina logarithmic potential at galactocentric distances of 4 kpcand 8 kpc for the strong and the weak field case respectively.We emphasize that the set of initial conditions ex-plored here are not meant to provide a complete coverage ofthe possible range of the various structural and kinematicproperties, and other parameters characterizing multiple-population clusters. The goal of this paper is to study thefundamental aspects of the early and long-term evolution of multiple-population clusters, and illustrate the role of vari-ous processes in the dynamics of multiple populations. Anextension of the set of simulations presented here is currentlyin progress. In the extended suite of simulations, besides con-sidering a broader range of initial structural properties, wewill include a population of primordial binaries. Althoughprimordial binaries are not expected to significantly affectthe general dynamical aspects studied in this paper, as dis-cussed in the Introduction the study of the differences in theevolution and survival of 1G and 2G binaries can provide akey additional dynamical insight into the formation and dy-namical of multiple-population clusters. Another importantaspect that will be further explored in a future study con-cerns the evolution of rotating multiple-population clusters.The Monte Carlo code used in this study can not follow theevolution of rotating stellar systems and the study of rotat-ing systems will be pursued by N-body simulations (see e.g.Mastrobuono-Battisti & Perets 2013, 2016, Henault-Brunetet al. 2015, Tiongco et al. 2019 for some early studies ofrotating multiple-population clusters). The planned exten-sions of the study presented here will allow us to furtherexplore the origin of the trends observed in the properties ofmultiple populations in the Galactic globular cluster systemand their possible connection with specific range of initialproperties set by the formation processes. Many formation models predict that the fraction of the totalcluster’s mass in 2G stars, M / M tot , was initially smallerthan its current value and increased to reach its presentvalue as a result of the preferential loss of 1G stars dur-ing the cluster’s dynamical evolution. Here we explore thepossible dynamical path and processes behind the evolutionof M / M tot .In Figs. 1 and 2 we show the time evolution of M / M tot for a few representative models starting with a range of dif-ferent initial values of M / M tot , different tidal radii, anddifferent initial cluster’s masses and structural parameters.In all of our models most of the increase in M / M tot occurs during the first few Gyrs of a cluster’s evolution. Dur-ing this early evolutionary phase the effects of two-body re-laxation are negligible and the preferential loss of 1G starsleading to the early increase of M / M tot is driven by thecluster’s expansion triggered by mass loss due to stellar evo-lution. This early expansion and the dynamical effects in-duced by this mass loss have been thouroghly investigatedin the context of single-population clusters in several earlypioneering studies (see e.g. Chernoff & Shapiro 1987, Cher-noff & Weinberg 1990, Fukushige & Heggie 1995) which haveshown that this process can cause a cluster to lose a signif-icant number of stars and even undergo an early completedissolution (see also Vesperini et al. 2009, Whitehead et al.2013, Contenta et al. 2015, Giersz et al. 2019 for more re-cent studies further exploring different aspects of this pathto cluster dissolution). In the context of the evolution ofmultiple-population clusters we have discussed the possi-ble role of this dynamical process for the evolution of thefraction of 2G stars in our early simulations presented in MNRAS , 1–16 (2020)
E. Vesperini et al.
Table 1.
Summary of initial conditionsid. N M / M tot r h , / r h , ( W , G , W , G ) r h / r tidal r tidal sg01c20sf × × × × × × . . . . Time (Myr) M G M T O T sg01c20sfsg01c20wfsg025c20sfsg025c20wfsg025c20w04sfsg025c20w04wf sg025c10sfsg025c10wfsg04c20sfsg04c20wfsg04c10sfsg04c10wf Figure 1.
Time evolution of the ratio of the total mass in 2Gstars to the total mass of the cluster for models starting with stars: sg01c20sf (solid black line), sg01c20wf (dashed black line),sg025c20sf (solid purple line), sg025c20wf (dashed purple line),sg025c20w04sf (dotted purple line), sg025c20w04wf (dot-dashedpurple line), sg025c10sf (solid green line), sg025c10wf (dashedgreen line), sg04c20sf (solid orange line), sg04c20wf (dashed or-ange line), sg04c10sf (solid cyan line), sg04c10wf (dashed cyanline). D’Ercole et al. (2008). The Monte Carlo simulations pre-sented in this paper allow us to carry out a significantly moredetailed characterization of this early dynamical phase anda more extensive study of its dependence on a number ofinitial parameters.As shown in Figs. 1 and 2, the early expansion triggered . . . . Time (Myr) M G M T O T Figure 2.
Time evolution of the ratio of the total mass in 2G starsto the total mass of the cluster for models starting with × stars: 4m-sg01c20sf (solid black line), 4m-sg01c20wf (dashed blackline), 4m-sg025c20sf (solid purple line), 4m-sg025c20wf (dashedpurple line), 4m-sg025c10sf (solid green line), 4m-sg025c10wf(dashed green line). by stellar evolution mass loss affects primarily the more ex-tended 1G population while the centrally concentrated spa-tial distribution of the 2G subsystem prevents any signifi-cant loss of 2G stars. After the first few Gyrs of evolution,clusters enter a dynamical phase driven by the effects of two-body relaxation; during this phase the loss of stars due totwo-body relaxation affects both 1G and 2G stars with aslight preferential loss of 1G stars causing M / M tot to con-tinue to slightly increase. In all the cases investigated the MNRAS000
Time evolution of the ratio of the total mass in 2G starsto the total mass of the cluster for models starting with × stars: 4m-sg01c20sf (solid black line), 4m-sg01c20wf (dashed blackline), 4m-sg025c20sf (solid purple line), 4m-sg025c20wf (dashedpurple line), 4m-sg025c10sf (solid green line), 4m-sg025c10wf(dashed green line). by stellar evolution mass loss affects primarily the more ex-tended 1G population while the centrally concentrated spa-tial distribution of the 2G subsystem prevents any signifi-cant loss of 2G stars. After the first few Gyrs of evolution,clusters enter a dynamical phase driven by the effects of two-body relaxation; during this phase the loss of stars due totwo-body relaxation affects both 1G and 2G stars with aslight preferential loss of 1G stars causing M / M tot to con-tinue to slightly increase. In all the cases investigated the MNRAS000 , 1–16 (2020) ynamics of multiple populations . . . . . . Time (Myr) MM Figure 3.
Time evolution of the total mass of 1G stars (redline), 2G stars (blue line), and of the total cluster’s mass (blackline) each normalized to their respective initial values for thesg025c20sf model. final values of M / M tot fall within the range of those foundin observational studies showing that for Galactic clusters M / M tot ∼ . − . (see Milone et al. 2017).The results of our simulations also clearly illustrate thedependence of the evolution of M / M tot on the structuralproperties of the 1G system. For models in which the 1Gstarts with a lower concentration density profile (modelssg025c20w04sf and sg025c20w04wf), the initial increase in M / M tot is more rapid and substantial as the early expan-sion of less concentrated 1G systems leads to a larger loss of1G stars (see the purple dotted and dot-dashed lines in Fig.1 for the evolution of the sg025c20w04sf and sg025c20w04wfmodels).In Fig. 3, we show the time evolution of the total massesof the 1G and 2G systems along with the total cluster’s mass(each normalized to their respective initial values) for one ofour models. The model shown in this figure is one of thosesignificantly affected by both early and long-term mass lossand serves to illustrate the early loss of stars suffered by the1G system and the transition from the early evolutionarystages dominated by the expansion triggered by mass lossdue to stellar evolution (and the loss of stars that ensuesfrom this expansion) to the long-term evolution driven bytwo-body relaxation.The link between the evolution of M / M tot and stel-lar escape is further illustrated in Fig. 4: this figure showsthe final values of M / M tot versus the ratio of the final tothe initial number of stars. Notice that, as expected, modelsevolving in a strong field are characterized by smaller val-ues of N f / N i ; this is a consequence of the stronger loss ofstars due to the effects of two-body relaxation during theclusters’ long-term evolution. We report in Table 2 the final l . . . . N f N init M G M T O T lll ll l l l Figure 4.
Final values of M / M tot versus the ratio of the final tothe initial number of stars: 4m-sg01c20sf (orange filled dot), 4m-sg01c20wf (orange square), 4m-sg025c20sf (brown filled dot), 4m-sg025c20wf (brown square), 4m-sg025c10sf (cyan filled dot), 4m-sg025c10sf (cyan square), sg01c20sf (blue filled dot), sg01c20wf(blue square), sg025c20sf (green filled dot), sg025c20wf (greensquare), sg025c20w04sf (filled yellow dot), sg025c20w04wf (yellowsquare), sg025c10sf (black filled dot), sg025c10wf (black square),sg04c20sf (red filled dot), sg04c20wf (red square), sg04c10sf (pur-ple filled dot), sg04c10wf (purple square). Dashed lines show thefull evolutionary tracks for a few representative models. values (at t = Gyr) of the total number of stars, N f , and M / M tot .As pointed out above, in our models, the loss of starsdue to the effects of two-body relaxation during the cluster’slong-term evolution causes only a slight additional increasein M / M tot compared to the more significant increase oc-curring during the cluster’s early evolution. The dominantrole played by early evolutionary processes in determiningthe current values of M / M tot implies that no significantdependence of M / M tot on the galactocentric distance isto be expected since, for the clusters with the tidally trun-cated structural properties adopted in our model, the extentof the early star loss depends very weakly on the strengthof the external tidal field. As shown in Figs. 1 and 2, thefinal values of M / M tot do not indeed depend significantlyon the strength of the external tidal field; for the modelsinvestigated in this paper, the differences between the valueof M / M tot at the end of the simulations for models in theweak field and those in the strong field range from ∼ to ∼ . .The small effect of stellar escape due two-body relax-ation on the evolution of M / M tot revealed by our simula-tions could possibly play the role of a second, and much lessimportant parameter in determining the current M / M tot of Galactic clusters. Possible observational evidence of thissecond-order dependence on the strength of the tidal fieldhas been suggested in the study by Zennaro et al. (2019) and MNRAS , 1–16 (2020)
E. Vesperini et al.
Table 2.
Final number of stars and 2G mass fractionid. N f M / M tot sg01c20sf 47,381 0.57sg01c20wf 111,517 0.57sg025c20sf 236,731 0.64sg025c20wf 352,947 0.58sg025c20w04sf 145,363 0.80sg025c20w04wf 239,498 0.77sg025c10sf 297,975 0.59sg025c10wf 388,698 0.55sg04c20sf 353,817 0.75sg04c20wf 477,971 0.70sg04c10sf 423,635 0.71sg04c10wf 524,268 0.684m-sg01c20sf 526,262 0.584m-sg01c20wf 651,973 0.534m-sg025c20sf 1,479,567 0.594m-sg025c20wf 1,727,911 0.544m-sg025c10sf 1,559,549 0.574m-sg025c10wf 1,699,293 0.54 Milone et al. (2020) who found that within the main correla-tion between fraction of 2G stars and cluster’s mass, clusterswith smaller pericentric distances tend to have slightly larger2G fractions.Finally we emphasize that while in this work the expan-sion of clusters triggered by mass loss due to stellar evolutionis the process responsible for the clusters’ early expansionand loss of 1G stars and evolution of M / M tot , there area number of additional dynamical ingredients not includedin our simulations which may further strengthen the earlyloss of stars. In particular, as shown in a number of stud-ies, the complex external environment and tidal field duringthe first few Gyrs of clusters’ evolution can lead to a signif-icant additional loss of stars and further contribute to theevolution of M / M tot (see e.g. Li & Gnedin 2019, Carlberg2020, Renaud 2018, 2020 and references therein). Primor-dial gas expulsion (see e.g. Boily & Kroupa 2003) could alsocontribute to a cluster’s early expansion and loss of stars. Fi-nally, primordial mass segregation can strengthen the earlyexpansion due to stellar evolution and further enhance theearly rate of star escape (see e.g. Vesperini et al. 2009, Haghiet al. 2014).In order to further explore and illustrate the role ofearly and long-term evolutionary processes in driving theevolution of M / M tot , Fig. 5 shows for a few representativemodels the evolution of M / M tot versus the power-law in-dex of the stellar mass function, α . − . , calculated from amaximum-likelihood fit of the stellar mass function for mainsequence stars with masses between 0.3 m (cid:12) and 0.8 m (cid:12) .As shown in several studies (see e.g. Vesperini & Heg-gie 1997, Baumgardt & Makino 2003, Trenti et al. 2010),stellar escape induced by the effects of two-body relaxationpreferentially affects low-mass stars and causes a gradualflattening of the stellar mass function (hereafter MF) as acluster evolves and loses stars. A correlation between thefraction of mass lost due to two-body relaxation and theslope of the MF has been found in several numerical stud-ies (see e.g. Vesperini & Heggie 1997, Trenti et al. 2010,Webb & Leigh 2015 and references therein; see also Sollimaet al. 2017, Ebrahimi et al. 2020 for observational studies). −2.0 −1.8 −1.6 −1.4 −1.2 −1.0 −0.8 . . . . a - M G M t o t lll l l l l l l l l l l l llll sg01c20sfsg01c20wfsg025c20sfsg025c20wfsg025c20w04sfsg025c20w04wf llll l lll l l l l l llllll ll l l l l l l lllllllllllllll −2.0 −1.8 −1.6 −1.4 −1.2 −1.0 −0.8 ( M G M t o t )( M G M t o t ) i n i t llll l l l l l l l l l l a - lllllllllll l l llllllllllll ll lllllllllllllll Figure 5.
Evolution of the ratio of the total mass of the 2G sys-tem to the total cluster’s mass, M / M tot , versus the power-lawindex of the stellar mass function for main sequence stars withmasses between 0.3 m (cid:12) and 0.8 m (cid:12) , α . − . , for the following mod-els: sg01c20sf (solid black line), sg01c20wf (dashed black line),sg025c20sf (solid purple line), sg025c20wf (dashed purple line),sg025c20w04sf (dotted purple line), sg025c20w04wf (dot-dashedpurple line). Points with different symbols overplotted on eachline are show the values of M / M tot and α . − . at times be-tween 0 Gyr and 13 Gyr with a time step of 1 Gyr. The figure inthe inset shows the same data but with the values of M / M tot normalized to their initial values. Assuming that the slope of the initial MF is known, thepresent-day MF slope can thus be used to estimate the frac-tion of stars lost due to two-body relaxation (Webb & Leigh2015, Baumgardt et al. 2019).Early loss of stars induced by a cluster’s expansion trig-gered by stellar evolution (or other early dynamical pro-cesses), on the other hand, does not significantly affect thecluster’s stellar mass function: even the loss of a large num-ber of stars during a cluster’s early evolution phases mightleave no fingerprint on the slope of the stellar MF in thelow-mass range (typically 0.3-0.8 m (cid:12) ) observed in old globu-lars clusters (see e.g. Vesperini et al. 2009). A very extremeprimordial mass segregation would be needed to affect theMF in this low-mass range with the early star loss discussedhere. This aspect will be further explored in a future inves-tigation.Fig. 5 illustrates these points and their connection withthe evolution of M / M tot : this figure shows that all themodels are characterized by an initial phase during which M / M tot increases while α . − . is approximately constant.During this phase, the escape of a significant fraction of 1Gstars leads to an increase in M / M tot but no significantevolution in α . − . . The evolutionary tracks on the α . − . - M / M tot plane clearly show that most of the loss of starsneeded to increase M / M tot from its initial value to valuesclose to those currently observed does not significantly affect MNRAS000
Evolution of the ratio of the total mass of the 2G sys-tem to the total cluster’s mass, M / M tot , versus the power-lawindex of the stellar mass function for main sequence stars withmasses between 0.3 m (cid:12) and 0.8 m (cid:12) , α . − . , for the following mod-els: sg01c20sf (solid black line), sg01c20wf (dashed black line),sg025c20sf (solid purple line), sg025c20wf (dashed purple line),sg025c20w04sf (dotted purple line), sg025c20w04wf (dot-dashedpurple line). Points with different symbols overplotted on eachline are show the values of M / M tot and α . − . at times be-tween 0 Gyr and 13 Gyr with a time step of 1 Gyr. The figure inthe inset shows the same data but with the values of M / M tot normalized to their initial values. Assuming that the slope of the initial MF is known, thepresent-day MF slope can thus be used to estimate the frac-tion of stars lost due to two-body relaxation (Webb & Leigh2015, Baumgardt et al. 2019).Early loss of stars induced by a cluster’s expansion trig-gered by stellar evolution (or other early dynamical pro-cesses), on the other hand, does not significantly affect thecluster’s stellar mass function: even the loss of a large num-ber of stars during a cluster’s early evolution phases mightleave no fingerprint on the slope of the stellar MF in thelow-mass range (typically 0.3-0.8 m (cid:12) ) observed in old globu-lars clusters (see e.g. Vesperini et al. 2009). A very extremeprimordial mass segregation would be needed to affect theMF in this low-mass range with the early star loss discussedhere. This aspect will be further explored in a future inves-tigation.Fig. 5 illustrates these points and their connection withthe evolution of M / M tot : this figure shows that all themodels are characterized by an initial phase during which M / M tot increases while α . − . is approximately constant.During this phase, the escape of a significant fraction of 1Gstars leads to an increase in M / M tot but no significantevolution in α . − . . The evolutionary tracks on the α . − . - M / M tot plane clearly show that most of the loss of starsneeded to increase M / M tot from its initial value to valuesclose to those currently observed does not significantly affect MNRAS000 , 1–16 (2020) ynamics of multiple populations the slope of the initial MF in the mass range observed in oldclusters; in general, the slope of the present-day MF can notbe used to infer the extent of any episode of early star lossor to rule out its role in determining the present value of M / M tot . The figure in the inset shows the same data butwith the values of M / M tot normalized to their initial val-ues to further illustrate the relative variation of M / M tot during the different evolutionary phases.Some of the models included in Fig. 5 lose a non-negligible number of stars due to the effects of two-body re-laxation during their long-term evolution and the MF grad-ually flattens as a consequence of the preferential loss of low-mass stars. Examples of these systems include the sg01c20sf,sg01c20wf, sg025c20sf, and the sg025c20w04sf models. Theevolution of these systems is characterized by the presenceof an almost horizontal portion of the track on the α . − . - M / M tot plane; during that phase the loss of stars driven bytwo-body relaxation preferentially removes low-mass starsand affects α . − . but, as shown in Fig. 5 and discussedabove, this loss of stars only slightly affects M / M tot .We conclude this section with a brief discussion of thefraction of 2G stars among the escaping stars that populatea cluster’s tidal tails. Several recent investigations have sig-nificantly extended the early observational studies (see e.g.Leon et al. 2000, Grillmair & Dionatos 2006) of globular clus-ters’ tidal tails (see e.g. Carballo-Bello et al. 2018, Grillmair2019, Ibata et al. 2019, 2020, Kaderali et al. 2019, Bonacaet al. 2020a, 2020b, Thomas et al. 2020, Lane et al. 2020,Sollima 2020, Starkman et al. 2020) and revealed a varietyof complex morphological and kinematic features that canshed light on the clusters’ dynamics and the properties ofthe Galactic halo (see e.g. Sanderson et al. 2017, Bonaca &Hogg 2018, Bonaca et al. 2019, Webb & Bovy 2019, Reinoet al. 2020). All the observational studies of multiple popu-lations have so far focussed their attention on the clusters’more populated inner regions; future photometric and spec-troscopic studies focussed on the properties of stars in thetidal tails could shed light on several key aspects of the dy-namics of multiple populations although the relatively smallnumber of stars spread in extended regions make these stud-ies particularly challenging (see e.g. Ji et al. 2020, Chun et al.2020, Hansen et al. 2020, for recent studies aimed at charac-terizing the chemical properties of extratidal stars and starsin tidal streams possibly originating from globular clusters).Tidal tails found in observational studies are typically popu-lated by stars escaped from a cluster in the last ≈ . − Gyr(depending on the extension of the tail revealed by the obser-vations and the cluster’s orbit). Here we focus our attentionon stars escaping from clusters during a more extended timeinterval in order to better illustrate the link with the evolv-ing fraction of 2G stars inside the cluster and the spatialmixing of the two populations. In Fig. 6 we show, for a fewmodels, ( N / N tot ) esc (defined as the fraction of 2G stars inthe population of stars escaping from a cluster in a 1 Gyrtime interval) measured at different times between 8 and 13Gyr versus the fraction of 2G stars inside the cluster at thesame times (top panel) and the fraction of 2G stars insidethe cluster’s half-mass radius (bottom panel). For example,a point corresponding at t = Gyr shows on the x-axis thefraction of 2G stars inside the cluster at t = Gyr and, onthe y-axis, the fraction of 2G stars in population of stars that escaped from the cluster between t = Gyr and t = Gyr.In all the models presented here, at t = Gyr clus-ters have already entered the phase dominated by relaxationmass loss. As discussed above, until the two populations arespatially mixed, mass loss due to relaxation still preferen-tially removes 1G stars although the preferential loss of 1Gstars is much weaker than that during the cluster’s earlyevolutionary phases dominated by the expansion triggeredby stellar evolution effects. The preferential loss of 1G starsimplies that, in general, the fraction of 2G stars in the tidaltails is smaller than the fraction of 2G stars inside the clus-ters; the two fractions are identical only for clusters in whichthe two populations are completely mixed. Fig. 6 shows thatthis is indeed the case and that the fraction of 2G stars in thepopulation of escaping stars is in general smaller than theglobal 2G fraction inside the cluster. The difference betweenthe 2G fraction in the escaping population and inside thecluster depends on the degree of spatial mixing reached bythe cluster. For example, in model sg01c20sf, the two pop-ulations are completely mixed at the end of the simulationand the two fractions are indeed identical. For all the otherclusters ( N / N tot ) esc is smaller than the fraction of 2G starsinside the cluster.The bottom panel of Fig. 6 shows ( N / N tot ) esc versusthe 2G fraction inside the cluster’s half-mass radius. Thisfigure may be more useful for a comparison with future ob-servational data for clusters that are not completely mixedsince the fraction inside the cluster half-mass is often thequantity measured in observational studies.Finally we point out that clusters that do not form 2Gstars or have a small 2G fraction (see e.g. Conroy & Spergel2011, Dondoglio et al. 2020 for some studies that have ex-plored the possible mass threshold for the 2G formation)could completely dissolve and leave a stream populated onlyor mainly by 1G stars. Previous studies exploring the formation of multiple popu-lations by means of hydro simulations (see e.g. D’Ercole etal. 2008, Bekki 2010, 2011, Calura et al. 2019) have shownthat 2G stars form in a compact subsystem in the central re-gions of the 1G system resulting in a cluster characterized bycomplex multi-scale structural and kinematical properties.In Fig. 7 we show the time evolution of the ratio of the1G half-mass radius to the 2G half-mass radius, r h , / r h , ,for a few representative models. This figure clearly showstwo distinct phases in the spatial mixing process.The first phase, occurring during the early evolutionarystages when the cluster is losing a large fraction of its 1Gpopulation, leads to a significant decrease in r h , / r h , anda rapid (in a few Gyrs) evolution away from the initial struc-tural properties set at the end of the formation process. Forall the models studied in this paper, at the end of this firstevolutionary phase r h , / r h , is significantly different fromits initial value. For most of the models presented here thevalue of r h , / r h , at the end of this phase falls approxi-mately in the range ≈ − . This early phase is followed by amore gradual and slower decrease of r h , / r h , as the clus-ter enters the long-term evolutionary stage in which mixing MNRAS , 1–16 (2020)
E. Vesperini et al. . . . . . . N N tot ( N G N t o t ) e sc l lllll l l lllll l l l l ll l l llll l l lll l l l l lll l l l l l l l llll llllllll sg01c20sfsg01c20wfsg025c10sfsg025c10wfsg025c20sfsg025c20wfsg025c20w04sfsg04c20sf . . . . . . N N tot ( r < r h ) ( N G N t o t ) e sc llllll llllll llllllllllll llllll llllll llllll llllll Figure 6.
Evolution of the number fraction of escaping stars be-longing to the 2G in time intervals of 1 Gyr between t = Gyrand t = Gyr versus the global fraction of 2G stars inside thecluster (top panel) and versus the fraction of 2G stars inside acluster’s half-mass radius (bottom panel). For each line the pointwith the smallest (largest) value of escapers ( N / N tot ) esc cor-responds to t = Gyr ( t = Gyr). For each point, the valueof ( N / N tot ) esc represents the fraction of 2G stars in the pop-ulations of stars that escaped between time t Gyr and ( t − ) Gyr and the value of ( N / N tot ) is the fraction of 2G stars in-side the cluster (top panel) or inside the cluster’s half-mass ra-dius (bottom panel) at time t Gyr. Different lines correspond tothe following models: sg01c20sf (red line), sg01c20wf (blue line),sg025c10sf (purple line), sg025c10wf (brown line), sg025c20sf (or-ange line), sg025c20wf (cyan line), sg025c20w04sf (black line),sg04c20sf (pink line). The thin dashed line represents values forwhich ( N / N tot ) esc = ( N / N tot ) . Time (Myr) r h1 G r h2 G sg01c20sfsg01c20wfsg025c20sfsg025c20wfsg025c20w04sfsg025c20w04wfsg025c10sfsg025c10wf Figure 7.
Time evolution of the ratio of the 1G half-mass radius tothe 2G half-mass radius for the following models: sg01c20sf (solidblack line), sg01c20wf (dashed black line), sg025c20sf (solid pur-ple line), sg025c20wf (dashed purple line), sg025c20w04sf (dottedpurple line), sg025c20w04wf (dot-dashed purple line), sg025c10sf(solid green line), sg025c10wf (dashed green line). A dotted hori-zontal line at r h , / r h , = is plotted to indicate the value of thisratio corresponding to the complete spatial mixing of the twopopulations. is driven by the effects of two-body relaxation. At the endof the simulations, some clusters retain some differences be-tween the spatial distributions of the 1G and the 2G spatialpopulations although in all cases these differences are muchsmaller than those in the initial conditions. Some clustersreach complete spatial mixing and, as discussed in Vesperiniet al. (2013) (see also Miholics et al. 2015, Henault-Brunetet al. 2015, Tiongco et al. 2019), complete mixing is reachedwhen clusters have lost a significant amount of mass duringits long-term evolution (in addition to any mass loss sufferedduring their early evolution). This is the case, for example,for model sg01c20sf.In order to further characterize these two phases in thespatial mixing process and their link to the different evolu-tionary stages and the mass loss occuring during a cluster’sevolution, we show in Fig. 8 the evolution of r h , / r h , as afunction of the slope the stellar MF, α . − . for a few repre-sentative models. The evolution on the α . − . - r h , / r h , plane clearly shows the two distinct mixing phases in allmodels presented. During the first phase r h , / r h , de-creases significantly while α . − . is approximately constant;this is the early evolutionary phase during which M / M tot is rapidly increasing as a result of the preferential loss of1G stars; as discussed in section 3.1, the mass loss occuringduring this early phase does not affect α . − . . The secondphase is instead characterized by mass loss and spatial mix-ing driven by the effects of two-body relaxation. The evo-lution of α . − . is the consequence of the preferential loss MNRAS000
Time evolution of the ratio of the 1G half-mass radius tothe 2G half-mass radius for the following models: sg01c20sf (solidblack line), sg01c20wf (dashed black line), sg025c20sf (solid pur-ple line), sg025c20wf (dashed purple line), sg025c20w04sf (dottedpurple line), sg025c20w04wf (dot-dashed purple line), sg025c10sf(solid green line), sg025c10wf (dashed green line). A dotted hori-zontal line at r h , / r h , = is plotted to indicate the value of thisratio corresponding to the complete spatial mixing of the twopopulations. is driven by the effects of two-body relaxation. At the endof the simulations, some clusters retain some differences be-tween the spatial distributions of the 1G and the 2G spatialpopulations although in all cases these differences are muchsmaller than those in the initial conditions. Some clustersreach complete spatial mixing and, as discussed in Vesperiniet al. (2013) (see also Miholics et al. 2015, Henault-Brunetet al. 2015, Tiongco et al. 2019), complete mixing is reachedwhen clusters have lost a significant amount of mass duringits long-term evolution (in addition to any mass loss sufferedduring their early evolution). This is the case, for example,for model sg01c20sf.In order to further characterize these two phases in thespatial mixing process and their link to the different evolu-tionary stages and the mass loss occuring during a cluster’sevolution, we show in Fig. 8 the evolution of r h , / r h , as afunction of the slope the stellar MF, α . − . for a few repre-sentative models. The evolution on the α . − . - r h , / r h , plane clearly shows the two distinct mixing phases in allmodels presented. During the first phase r h , / r h , de-creases significantly while α . − . is approximately constant;this is the early evolutionary phase during which M / M tot is rapidly increasing as a result of the preferential loss of1G stars; as discussed in section 3.1, the mass loss occuringduring this early phase does not affect α . − . . The secondphase is instead characterized by mass loss and spatial mix-ing driven by the effects of two-body relaxation. The evo-lution of α . − . is the consequence of the preferential loss MNRAS000 , 1–16 (2020) ynamics of multiple populations −1.8 −1.6 −1.4 −1.2 −1.0 −0.8 a - r h1 G r h2 G sg01c20sfsg01c20wfsg025c20sfsg025c20w04sfsg025c10sf Figure 8.
Evolution of the ratio of the 1G half-mass radius to the2G half-mass radius versus the slope of the mass function be-tween 0.3 and 0.8 m (cid:12) , α . − . , for the following models: sg01c20sf(black line), sg01c20wf (orange line), sg025c20sf (purple line),sg025c20w04sf (green line), sg025c10sf (cyan line). A dotted hor-izontal line at r h , / r h , = is plotted to indicate the value of thisratio corresponding to the complete spatial mixing of the twopopulations. l l l . . . . . . mass bin0.1−0.3 0.3−0.6 0.6−0.8 r h1 G r h2 G l l ll l ll l ll l ll l l sg01c20sfsg01c20wfsg025c20wfsg025c20w04sfsg025c20w04wfsg025c10wf Figure 9.
Ratio of the 1G half-mass radius to the 2G half-massradius at t = Gyr for main sequence stars in three differentmass bins for the following models: sg01c20sf (solid black line),sg01c20wf (dashed black line), sg025c20wf (dashed purple line),sg025c20w04sf (dotted purple line), sg025c20w04wf (dot-dashedpurple line), sg025c10wf (dashed green line). of low-mass stars during this second phase; r h , / r h , con-tinues to decrease although much more slowly. As shown inthis figure, the path followed by all clusters in the α . − . - r h , / r h , plane shows little variation among the modelsstudied; this is consistent with the results of Vesperini et al.(2013) showing that the long-term mixing of the 1G and 2Gspatial distributions are closely linked with mass loss occur-ring during the cluster’s long-term evolution.We strongly emphasize that using the slope of thepresent-day MF as a general indicator of the extent of massloss due to two-body relaxation and of the degree of spatialmixing in different clusters relies on the assumption that allclusters formed with a universal initial MF; if all clustersformed with a universal initial MF, the current value of theglobal MF’s slope is determined by the cluster’s dynamicalhistory and loss of stars during its long-term evolution andcan be used to estimate the expected extent of spatial mix-ing. On the other hand, should the stellar initial MF notbe universal (see e.g. Kroupa 2020 and references therein),caution must be used in linking the expected spatial mixingwith the slope of the mass function; for example, if a clusterformed with an initial MF flatter than a standard Kroupa(2001) mass function (see e.g. Zonoozi et al. 2011, Giersz &Heggie 2011, Webb et al. 2017, Henault-Brunet et al. 2020,Cadelano et al. 2020 for some examples of clusters for whichthis could be the case), a flat slope of the present-day MFdoes not imply that a cluster suffered strong mass loss dueto relaxation and that the different stellar populations mustbe completely mixed.As already discussed in Vesperini et al. (2018), the moreconcentrated spatial distribution of the 2G population im-plies a more rapid segregation of massive 2G stars and diffu-sion towards the outer regions of the low-mass 2G stars. Oneof the implications of this difference is that the timescale of1G-2G spatial mixing depends on the stellar mass: low-massstars tend to mix more rapidly than massive stars. Fig. 9shows r h , / r h , for stars in different mass bins for a few ofthe models studied in this paper: with the exception of themodel sg01c20sf which has reached complete spatial mixing,all the other models show that at t = Gyr r h , / r h , depends on the stellar mass although the trend we find isprobably too weak to be detected observationally.A more detailed characterization of the internal struc-tural properties of the two populations is provided by thetime evolution of the ratio of the 2G to the 1G surface num-ber density profiles in Fig. 10. We show the surface densityprofiles to establish a closer contact with quantities that canbe determined observationally; all the profiles are calculatedincluding only main sequence stars with masses between 0.1and 0.8 m (cid:12) . The radial profile in Fig. 10 and in all the subse-quent figures show the median profile calculated from thirtydifferent random 2D projections. Each profile is normalizedto the global 2G-to-1G number ratio so that complete mix-ing corresponds to a flat profile with a value equal to 1 atall radii. This figure further illustrates the details of the dy-namical path toward spatial mixing of the two populations.The evolution of the ratio of the 2G to the 1G densityprofile is driven by a combination of mixing of the two pop-ulations and the evolution of the global ratio of the numberof 2G to 1G stars. Before complete mixing is reached thecluster 2G-to-1G surface density ratio can be characterizedby a flat portion both in the innermost and in the outermost MNRAS , 1–16 (2020) E. Vesperini et al. . . . . . . . R R hl ( S G S G )( N G N G ) Figure 10.
Radial profiles of the ratio of the 2G to the 1G surfacenumber density profile (normalized to the global 2G to 1G numberratio) for the sg01c20sf model at t=0 (black line), 0.5 Gyr (redline), 2 Gyr (blue line) 4 Gyr (cyan line), 8 Gyr (purple line),and 12 Gyr (orange line). Radius is normalized to the cluster’shalf-light radius at the time the profile is calculated. regions: the two populations are not spatially mixed yet buttheir surface density profiles share a similar variation withradius in those regions leading to a flat profile of the 2G-to-1G density ratio (in some cases the outermost radial profileof the 2G profile can be slightly shallower than that of the1G leading to a slightly increasing 2G-to-1G surface densityratio). Fig. 10 also clearly shows the two phases of mixingassociated with the early and long-term evolution alreadydiscussed above.Finally in Fig. 11 we show the final density profile for afew representative models which have reached different de-grees of spatial mixing after 12 Gyr of evolution; Fig. 12shows the 2G-to-1G ratio of the surface density profiles at t = Gyr for the three models shown in Fig. 11. In thesg01c20sf model the two populations are essentially com-pletely mixed and follow similar density profiles. The othertwo models shown in Fig. 11, on the other hand, are notcompletely mixed after 12 Gyr. In the sg025c20sf model thedensity profiles of the 2G and the 1G populations have asimilar shape in the inner regions ( R (cid:46) . R hl ) (which there-fore corresponds to an approximately flat ratio of Σ G / Σ G ;see the purple line in Fig. 12); for R (cid:38) . R hl the Σ G de-creases more rapidly than Σ G and finally in the outermostregions ( R (cid:38) R hl ) the two populations follow again a similarprofile (which, again, corresponds to an approximately flatratio of Σ G / Σ G ; see the purple line in Fig. 12).The sg025c20wf model is the system that is the farthestfrom complete mixing among the three shown in Fig. 11: forthis model the 2G and the 1G populations follow differentdensity profiles and only in the outermost regions ( R > R hl ) the two profiles have similar variation with radius resultingin a flat Σ G / Σ G profile (see orange line in Fig. 12). As discussed in the Introduction, the study of the kinematicproperties of multiple stellar populations is still in its earlystages but a few investigations (see e.g. Richer et al. 2013,Bellini et al. 2015, 2018, Cordero et al. 2017, Dalessandroet al. 2018, Cordoni et al. 2020a, 2020b, Lee 2020) havestarted to address this aspect and further enriched our un-derstanding of the dynamical picture of multiple populationsin globular clusters.Here we focus our attention on the evolution of theanisotropy in the velocity distribution, and the dependenceof the velocity distribution on the stellar mass as the systemevolves toward energy equipartition.
Figs. 13-17 show the projected radial profile of theanisotropy defined as σ T / σ R − , where σ T and σ R are, re-spectively, the radial and tangential velocity dispersion de-fined on the 2-D plane of projection.For the three models shown in Fig. 13-15, both popula-tions start as isotropic King models but while the 1G is ini-tially tidally filling, the 2G is initially confined in a compactand tidally underfilling subsystem. These initial differencesin the structural properties imply that the 2G develops aradially anisotropic velocity distribution as it diffuses fromthe inner to the outer regions while the 1G remains approxi-mately isotropic (see also the discussion in Bellini et al. 2015,Tiongco et al. 2016, 2019). As discussed in Tiongco et al.(2016), in systems initially underfilling, the radial anisotropydeveloped during a cluster’s evolution eventually starts todecrease as the cluster continues its evolution and can becompletely erased for clusters that within a Hubble timereach the very advanced stages of their evolution and lose alarge fraction of their mass.The three models shown in Figs.13-15 show theanisotropy profiles of three clusters which have reached dif-ferent stages of their dynamical evolution at t = Gyr: thesg01c20sf (Fig. 13) is the model in the most advanced dy-namical phase (as illustrated also by the fact that at t = Gyr the two populations are spatially mixed; see Fig. 12):the radial anisotropy developed by the 2G during its evolu-tion has been erased by the effects of relaxation and massloss and both populations have a now similar anisotropy pro-file (isotropic in the inner regions and slightly tangentiallyanisotropic in the outer regions).The other two models shown in Fig. 14 and Fig. 15 arein a less advanced dynamical phase and in both of themthe 2G shows the signature of the anisotropy radial pro-file developed during its previous evolution and diffusiontowards the cluster’s outer regions: the 2G is isotropic inthe innermost regions ( R < R hl ) and radially anisotropic inthe intermediate/outer regions. In the outermost regions,the radial anisotropy of the 2G population decreases againand the velocity distribution turns into an approximatelyisotropic/tangentially anistropic distribution.It is important to emphasize that while in these models MNRAS000
Figs. 13-17 show the projected radial profile of theanisotropy defined as σ T / σ R − , where σ T and σ R are, re-spectively, the radial and tangential velocity dispersion de-fined on the 2-D plane of projection.For the three models shown in Fig. 13-15, both popula-tions start as isotropic King models but while the 1G is ini-tially tidally filling, the 2G is initially confined in a compactand tidally underfilling subsystem. These initial differencesin the structural properties imply that the 2G develops aradially anisotropic velocity distribution as it diffuses fromthe inner to the outer regions while the 1G remains approxi-mately isotropic (see also the discussion in Bellini et al. 2015,Tiongco et al. 2016, 2019). As discussed in Tiongco et al.(2016), in systems initially underfilling, the radial anisotropydeveloped during a cluster’s evolution eventually starts todecrease as the cluster continues its evolution and can becompletely erased for clusters that within a Hubble timereach the very advanced stages of their evolution and lose alarge fraction of their mass.The three models shown in Figs.13-15 show theanisotropy profiles of three clusters which have reached dif-ferent stages of their dynamical evolution at t = Gyr: thesg01c20sf (Fig. 13) is the model in the most advanced dy-namical phase (as illustrated also by the fact that at t = Gyr the two populations are spatially mixed; see Fig. 12):the radial anisotropy developed by the 2G during its evolu-tion has been erased by the effects of relaxation and massloss and both populations have a now similar anisotropy pro-file (isotropic in the inner regions and slightly tangentiallyanisotropic in the outer regions).The other two models shown in Fig. 14 and Fig. 15 arein a less advanced dynamical phase and in both of themthe 2G shows the signature of the anisotropy radial pro-file developed during its previous evolution and diffusiontowards the cluster’s outer regions: the 2G is isotropic inthe innermost regions ( R < R hl ) and radially anisotropic inthe intermediate/outer regions. In the outermost regions,the radial anisotropy of the 2G population decreases againand the velocity distribution turns into an approximatelyisotropic/tangentially anistropic distribution.It is important to emphasize that while in these models MNRAS000 , 1–16 (2020) ynamics of multiple populations − + + + + R R hl S ( R ) − + + + + R R hl S ( R ) − + + + + R R hl S ( R ) Figure 11.
Surface number density radial profile of the 2G (blue line), 1G (red line) and total (black line) at t = Gyr for the sg01c20sf(left panel), sg025c20sf (middle panel), and sg025c20wf (right panel) models.
R R hl ( S G S G )( N G N G ) sg01c20sfsg025c20sfsg025c20wf Figure 12.
Radial profiles of the ratio of the 2G to the 1G surfacedensity profile (normalized to the global 2G to 1G number ratio)at t = Gyr for the sg01c20sf (black line), sg025c20sf (purpleline), and sg025c20wf (orange line) models. the radial anisotropy develops during the cluster’s long-termevolution, it could also develop for both the 1G and the 2Gpopulations during the early formation and violent relax-ation phases of the cluster evolution (see e.g. Vesperini etal. 2014, Tiongco et al. 2016).In Fig. 16 and Fig. 17 we show the anisotropy profile fortwo systems starting with initial conditions characterized bya radially anisotropic velocity distribution. The two modelsstart with different degrees of initial anisotropy (one withanisotropy radius, r a , equal to r h / and the other equal to r h / ; see section 2) and reach different stages of their dy-namical evolution. For the model with r a = r h / (Fig. 16)mass loss and dynamical evolution have erased the initialanisotropy of the 1G population and the 1G is approxi- − . − . . . R R hl s T s R - Figure 13.
Radial profile of the velocity anisotropy of the 2G (blueline), 1G (red line) at t = Gyr for the sg01c20sf model. mately isotropic while the 2G is isotropic in the inner andoutermost regions and radially anisotropic in the intermedi-ate regions. In the model r a = r h / (Fig. 17), on the otherhand, the effects of dynamics and mass loss are not sufficientto completely erase the initial anisotropy of the 1G popu-lation. Also in this case the 2G is more radially anisotropicthan the 1G but the 1G itself is characterized by a moderateradial anisotropy.Finally, in Fig. 18 we show the projected radial profileof the ratio of the tangential velocity dispersions of the 1Gand the 2G populations and of the 2G to the 1G radial veloc-ity dispersions at t = Gyr for the models sg025c20sf andsg025c20wf (see Figs.14 and 15 for the anisotropy profilesof these models). This figure shows that the differences be-tween the 1G and 2G radial anisotropy is due mainly to thesmaller tangential velocity dispersion of the 2G population.
MNRAS , 1–16 (2020) E. Vesperini et al. − . − . . . R R hl s T s R - Figure 14.
Radial profile of the velocity anisotropy of the 2G (blueline), 1G (red line) at t = Gyr for the sg025c20sf model. − . − . . . R R hl s T s R - Figure 15.
Radial profile of the velocity anisotropy of the 2G (blueline), 1G (red line) at t = Gyr for the sg025c20wf model.
This is in agreement with what found in a few observationalstudies (Richer et al. 2013, Bellini et al. 2015, Milone et al.2018, Cordoni et al. 2020a; see also Bellini et al. 2015 andTiongco et al. 2019 for N-body simulations showing the sametrend).We will present in a future study a more comprehen-sive investigation of the evolution of clusters with differ-ent initial kinematic properties; the results presented here − . − . . . R R hl s T s R - Figure 16.
Radial profile of the velocity anisotropy of the 2G (blueline), 1G (red line) at t = Gyr for the model sg025c20wf withinitial anisotropy radius r a = r h / . − . − . . . R R hl s T s R - Figure 17.
Radial profile of the velocity anisotropy of the 2G (blueline), 1G (red line) at t = Gyr for the model sg025c20wf withinitial anisotropy radius r a = r h / . show that the models considered in this study and charac-terized by a tidally filling 1G and a compact 2G populationalways result in a 2G more radially anisotropic than the1G population (or with both populations characterized byan isotropic distribution if the cluster is its advanced dy-namical stages). The 1G population can be either isotropic MNRAS000
Radial profile of the velocity anisotropy of the 2G (blueline), 1G (red line) at t = Gyr for the model sg025c20wf withinitial anisotropy radius r a = r h / . show that the models considered in this study and charac-terized by a tidally filling 1G and a compact 2G populationalways result in a 2G more radially anisotropic than the1G population (or with both populations characterized byan isotropic distribution if the cluster is its advanced dy-namical stages). The 1G population can be either isotropic MNRAS000 , 1–16 (2020) ynamics of multiple populations or radially anisotropic depending on the cluster’s early andlong-term dynamical history. It is interesting to notice thatthe 1G and 2G anisotropy profiles in some of our modelsapproximately follow those observed in the few clusters inwhich the kinematics of multiple populations has been stud-ied (see e.g. Richer et al. 2013, Bellini et al. 2015, Milone etal. 2018, Cordoni et al. 2020a, 2020b): in some clusters bothpopulations are approximately isotropic or slightly tangen-tially anisotropic while in other clusters the 1G is approxi-mately isotropic at all radii explored while the 2G is radi-ally anisotropic with an anisotropy varying with the distancefrom the cluster’s centre as found in our simulations. The initial differences in the structural properties of multi-ple populations have dynamical implications for the evolu-tion of the 1G and 2G populations toward energy equiparti-tion. In order to measure the degree of energy equipartitionreached by a given population, we use the exponential func-tion introduced by Bianchini et al. (2016) to fit the varia-tion of the velocity dispersion, σ ( m ) , with the stellar mass m , σ ( m ) ∝ exp ( − . m / m eq ) . Smaller values of m eq correspondto higher degree of energy equipartition (see Bianchini etal. 2016 for further discussion). For our analysis we haveused to total velocity dispersion measured on a 2-D planeof projection and focussed our attention on main sequencestars with masses between 0.1 and 0.8 m (cid:12) . Fig. 19 showsthe radial variation of m eq at t = Gyr for two of the mod-els studied in this paper. As expected in both cases the in-ner regions are characterized by smaller values of m eq and,therefore, by a higher degree of energy equipartition. This isthe consequence of the fact that the inner regions are thosecharacterized by shorter relaxation times. The two models,however, show different behavior in the radial variation of m eq for the 1G and the 2G populations. As discussed in theprevious sections, the system sg01c20sf is in its advancedevolutionary stages: at t = Gyr the two populations arespatially mixed and characterized by an isotropic velocitydistribution. For this model the degree of energy equiparti-tion of the two populations is very similar for R (cid:46) R hl whilein the outer parts of the radial range explored the 2G ischaracterized by slightly smaller values of m eq .Model sg025c20sf, on the other hand, is in a less ad-vanced evolutionary stage, its populations are not com-pletely mixed (see Fig. 12) and are characterized by velocitydistributions with different levels of anisotropy (Fig. 14). Forthis system (see lower panel of Fig. 19), outside the cluster’sinnermost regions ( R (cid:38) . R hl ) the values of m eq for the 2Gpopulation are significantly smaller than those of the 1Gpopulation clearly showing that the initially more compact2G system is in a more advanced phase of its evolution to-ward energy equipartition than the 1G system. In this casethe difference between the values of m eq of the 1G populationand the 2G population provides an interesting signature ofthe initial structural differences of the two populations. Thevalues of m eq calculated for the 1G and 2G stars togetherare closer to those of the 2G in the inner regions where 2Gstars are the dominant population (since the two popula-tions are not mixed yet) and gradually become more similarto those of the 1G population in the outer regions which arepredominantly populated by 1G stars. . . . . . . . . R R hl s G s G . . . . . . . . R R hl s G s G Figure 18.
Radial profile of the ratio of the 2G to the 1G σ R (black line) and of the ratio of the 2G to the 1G σ T (green line) at t = Gyr for the sg025c20sf (top panel), and sg025c20wf (bottompanel) models.
In this paper, we have presented the results of a set ofMonte Carlo simulations following the dynamical evolutionof multiple-population clusters starting with different initialcluster masses, 2G mass fractions and density profiles, 2G-to-1G size ratios, and tidal radii.We explored the role of early and long-term dynami-cal processes in driving the evolution of the main internalstructural and kinematic properties of multiple populations,illustrated how the initial differences in the dynamical prop-erties of 1G and 2G stars evolve and how, in turn, they may
MNRAS , 1–16 (2020) E. Vesperini et al.
R R hl m eq R R hl m eq Figure 19.
Radial profile (with radius normalized to the half-ligh radius) of the equipartition mass, m eq , (see section 3.3.2) forthe sg01c20sf model (top panel) and the sg025c20sf model (lowerpanel) at t = Gyr. In each panel the values of m eq for the 1G(red line), the 2G (blue line) and the 1G and the 2G combined(black line) are shown. lead the 1G and 2G populations on different evolutionarypaths.The main results of this study are the following. • In our initial conditions, we have explored systems withan initial fraction of the total mass in 2G stars, M / M tot ,ranging from 0.1 to 0.4. In order for M / M tot to reach val-ues in the range of those observed ( ∼ . − . ), clustersneed to preferentially lose 1G stars. In our simulations mostof the evolution of M / M tot occurs during the first few Gyrs of a cluster’s evolution when the system is respondingto the mass loss due to stellar evolution by expanding andpreferentially losing 1G stars which are less centrally con-centrated than the 2G stars. At the end of this early phase, M / M tot has increased and reached values within the ob-served range (see Figs. 1-2). We show that this early loss ofstars does not affect the slope of the stellar mass function(see Fig. 5) and, therefore, mass loss estimates based on theslope of the present-day mass function can not be used toinfer the extent of this early episode of star loss. Additionalearly dynamical processes not included in our analysis (e.g.tidal shocks, expansion triggered by primordial gas expul-sion) could further increase the number of 1G stars escapingduring this early evolutionary phases and contribute to theevolution of M / M tot . • Once a cluster enters the long-term evolution phasedominated by the effects of two-body relaxation, therelaxation-driven loss of stars causes only a slight increasein M / M tot since in this phase clusters, while still prefer-entially losing 1G stars, start to lose also 2G stars (see Figs.1-5). • The final values of M / M tot for our models span arange between . and . and the final total number ofstars ranges from ∼ × to ∼ . × (see Fig. 4 andTable 2). • The fraction of 2G stars in the population of escapingstars, ( N / N tot ) esc , varies during a cluster’s evolution. Dur-ing the early phases of a cluster’s evolution the populationof escapers is dominated by 1G stars. Later in the cluster’sevolution ( N / N tot ) esc falls in the range 0.3-0.7 dependingon the fraction of 2G stars in the cluster as well as on thedegree of spatial mixing of the two populations. In all cases ( N / N tot ) esc is smaller than or equal to the fraction of 2Gstars inside the cluster (see Fig. 6). Clusters that do notform 2G stars or form only a small 2G population couldcompletely dissolve and leave a stream populated only ormainly by 1G stars. • As a cluster evolves, the 1G and 2G populations mixand the differences between the spatial distributions of thetwo populations decrease. We have identified two distinctphases in the spatial mixing process: an early phase whenthe ratio of the 1G to the 2G half-mass radii, r h , / r h , ,rapidly decreases followed by a more extended phase char-acterized by a slower mixing (see Figs. 7-8). We confirm theresults of our previous study and find that complete mixingis reached after the cluster has lost a significant fraction ofmass during its long-term evolution.We find that low-massstars mix more rapidly than more massive stars although theexpected differences between r h , / r h , for low-mass starsand more massive stars are small (see Fig. 9). For a fewrepresentative models we have presented the complete finalsurface density radial profile and shown in more detail thedifferences and similarities of the 1G and 2G density profiles(see Figs. 10-12). • We have studied the evolution of the cluster’s internalkinematics and found that the 2G population is character-ized by a radially anisotropic velocity distribution. The ex-tent of the 2G radial anisotropy varies with the distance fromthe cluster’s centre: the velocity distribution is isotropic inthe innermost regions, it becomes increasingly anisotropicat larger distances from the centre and finally becomesisotropic or slightly tangentially anisotropic again in the out-
MNRAS000
MNRAS000 , 1–16 (2020) ynamics of multiple populations ermost regions. The 1G population can either be isotropicat all clustercentric distances or be characterized by radialanisotropy profile similar to that of the 2G but with radialanisotropy which is in all cases weaker than that of the 2Gpopulation. The differences between the 1G and the 2G ra-dial anisotropy are due mainly to the differences betweenthe 2G and the 1G tangential velocity dispersion (see Fig.18): the 2G tangential velocity dispersion is smaller thanthat of the 1G population. The radial velocity dispersions ofthe two populations are similar. For clusters that have lost asignificant fraction of their mass due to two-body relaxationand are in the advanced stages of their evolution both the1G and the 2G populations have isotropic velocity distribu-tions at all clustercentric distances (see Figs. 13-17 for thevarious cases). • We have studied the evolution of the 1G and the2G populations toward energy equipartition and exploredthe implications of the different initial structural proper-ties of the two populations for their evolution toward en-ergy equipartition. For dynamically old clusters which havereached spatial and kinematic mixing, the degree of energyequipartition of the 1G is similar to that of the 2G popula-tions (see Fig. 19 upper panel). For clusters in less advancedstages of their evolution, however, we find that, with the ex-ception of the inner regions where the two populations havesimilar degrees of energy equipartition, the 2G population iscloser to energy equipartition than the 1G population (seeFig. 19 lower panel).In this paper we have presented an overview of themain dynamical processes driving the evolution of multiple-population clusters, explored the evolutionary paths fol-lowed by the 1G and 2G structural and kinematic proper-ties, and studied their relationship with various dynamicalparameters. In future works we will expand the range of ini-tial conditions considered and further explore the evolutionof the dynamical properties of multiple-population clustersand their dependence on the initial conditions.
ACKNOWLEDGEMENTS
EV acknowledges support from NSF grant AST-2009193.JH was supported by Basic Science Research Pro-gram through the National Research Foundation of Ko-rea (NRF) funded by the Ministry of Education (No.2020R1I1A1A01051827). MG and AH were partially sup-ported by the Polish National Science Center (NCN)through the grant UMO-2016/23/B/ST9/02732. AH is alsosupported by the Polish National Science Center grant Mae-stro 2018/30/A/ST9/00050.
DATA AVAILABILITY STATEMENT
The data presented in this article may be shared on reason-able request to the corresponding author.
REFERENCES
Bastian, N.; Lamers, H. J. G. L. M.; de Mink, S. E.; Longmore,S. N.; Goodwin, S. P.; Gieles, M., 2013, MNRAS, 436, 2398 Baumgardt H., Makino J., 2003, MNRAS, 340, 227Baumgardt H., Hilker, M., Sollima, A., Bellini, A., 2019, MNRAS,482, 5138Beccari G., Bellazzini M., Lardo C., Bragaglia A., Carretta E.,Dalessandro E., Mucciarelli A., Pancino, E., 2013, MNRAS,431, 1995Bellini A., Piotto, G.; Bedin, L. R.; King, I. R.; Anderson, J.;Milone, A. P.; Momany Y., 2009, A&A, 507, 1393Bellini A. et al., 2015, ApJ, 810, L13Bellini A. et al., 2018, ApJ, 853, 86Bekki K., 2010, ApJ, 724, L99Bekki K., 2011, MNRAS, 412, 2241Bekki K., Jerabkova T., Kroupa P., 2017, MNRAS, 471, 2242Bianchini P., van de Ven G., Norris M. A., Schinnerer E., VarriA. L., 2016, MNRAS, 458, 3644Binney J., Tremaine S., ’Galactic Dynamics’, 2008, PrincetonUniversity PressBoberg O. M., Friel E. D., Vesperini E., 2016, ApJ, 824, 5Boily C., Kroupa P., 2003, MNRAS,338, 673Bonaca A., Hogg D. W., 2018, ApJ, 867, 101Bonaca A., Hogg D. W., Price-Whelan A. M., Conroy C., 2019,ApJ, 880, 38Bonaca, A. et al. 2020a, ApJ, 889, 70Bonaca, A. et al. 2020b, submitted to AAS journals,arXiv:2012.09171Breen P. G., 2018, MNRAS, 481, L110Cadelano et al. , 2020, MNRAS, 499, 2390Calura F., D’Ercole A., Vesperini E., Vanzella E., Sollima A.,2019, MNRAS, 489, 3269Carlberg R.G., 2020, ApJ, 893, 116Carballo-Bello J. A.; Martinez-Delgado D., Navarrete C., CatelanM., Munoz R. R., Antoja T., Sollima A., 2018, MNRAS, 474,683Chernoff D.F., Shapiro S.L., 1987, ApJ, 322, 113Chernoff D.F., Weinberg M.D., 1990, ApJ, 351, 121Chun S., Lee J., Lim D., 2020, ApJ, 900, 146Conroy C., Spergel D.N., 2011, ApJ, 726, 36Contenta, F., Varri, A. L., Heggie, D. C., 2015, MNRAS, 449,L100Cordero, M. J.; Pilachowski, C. A.; Johnson, C. I.; McDonald, I.;Zijlstra, A. A.; Simmerer, J. , 2014, ApJ, 780, 94Cordero, M. J.; Pilachowski, C. A.; Johnson, C. I.; Vesperini, E.,2015, ApJ, 800, 3Cordero M. J., Henault-Brunet V., Pilachowski C. A., BalbinotE., Johnson C. I., Varri A. L., 2017, MNRAS, 465, 3515Cordoni G., et al. 2020a, ApJ, 889, 18Cordoni G., et al. 2020b, ApJ, 898, 147Dalessandro E., et al., 2014, ApJ, 791, L4Dalessandro E., et al., 2018, ApJ, 864, 33Dalessandro E., et al., 2019, ApJ, 884, L24D’Antona F., Vesperini E., D’Ercole A., Ventura P., Milone A.P., Marino A. F., Tailo M., 2016, MNRAS, 458, 2122Decressin T., Meynet G., Charbonnel C., Prantzos N., EkstroemS., 2007, A&A, 464, 1029de Mink S. E., Pols O. R., Langer N., Izzard R. G., 2009, A&A,507, L1Denissenkov P. A., Hartwick F. D. A., 2014, MNRAS, 437, L21D’Ercole A., Vesperini E., D’Antona F., McMillan S. L. W., Rec-chi S., 2008, MNRAS, 391, 825D’Ercole, A., D’Antona F., Ventura P. , Vesperini E., McMillanS.L.W., 2010, MNRAS, 407, 854D’Ercole, A., D’Antona, F., Carini, R., Vesperini, E., Ventura P.,2012, MNRAS, 423, 1521D’Ercole, A., D’Antona, F., Vesperini, E., 2016, MNRAS, 461,4088Dondoglio E., et al. , 2020, ApJ, in press (arXiv:2011.03283)Ebrahimi, H.; Sollima, A.; Haghi, H.; Baumgardt, H.; Hilker, M.,2020, MNRAS, 494, 4226MNRAS , 1–16 (2020) E. Vesperini et al.
Elmegreen B., 2017, ApJ, 836, 80Fare A., Webb J. J., Sills A., 2018, MNRAS, 481, 3027Fukushige T., Heggie, D.C., 1995, MNRAS, 276, 206Gerber J. M., Friel, E. D., Vesperini E., 2018, AJ, 156, 6Gerber J. M., Friel, E. D., Vesperini E., 2020a, AJ, 159, 50Gerber J. M., Friel, E. D., Vesperini E., 2020b, submitted to AASjournalsGieles M., et al. 2018, MNRAS, 478, 2461Giersz M., Heggie D.C., 2011, MNRAS, 410, 2698Giersz M., Heggie D. C., Hurley J. R., Hypki A., 2013, MNRAS,431, 2184Giersz M., Askar A., Wang L., Hypki A., Leveque A., SpurzemR., 2019, MNRAS, 487, 2412Gratton R., Bragaglia A., Carretta E., D’Orazi V., Lucatello S.,Sollima A., 2019,
A&ARv , 27, 8Grillmair C. J., Dionatos O., 2006, ApJ, 643, L17Grillmair C. J., 2019, ApJ, 884, 174Haghi H., Hoseini-Rad S., Zonoozi A. H., Kuepper A. H. W.,2014, MNRAS, 444, 3699Hansen, T. T., Riley, A. H., Strigari, L. E., Marshall, J. L., Fer-guson, P. S., Zepeda, J., Sneden, C., 2020, ApJ, 901, 23Henault-Brunet, V.; Gieles, M.; Agertz, O.; Read, J. I., 2015,MNRAS, 450, 1164Henault-Brunet V., Gieles M., Strader J., Peuten M., BalbinotE., Douglas K. E. K., 2020, MNRAS, 491, 113Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, MNRAS,360, 974Hong J., Vesperini E., Sollima A., McMillan S. L. W., D’AntonaF., D’Ercole A., 2015, MNRAS, 449, 629Hong J., Vesperini E., Sollima A., McMillan S. L. W., D’AntonaF., D’Ercole A., 2016, MNRAS, 457, 4507Hong J., Patel S., Vesperini E., Webb J. J., Dalessandro E., 2019,MNRAS, 483, 2592Howard C. S., Pudritz R. E., Sills A., Harris W. E., 2019, MN-RAS, 486,1146Hypki A., Giersz, M., 2013, MNRAS, 429, 1221Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897Ibata R.A., Bellazzini M., Malhan K., Martin N., Bianchini P.,2019, Nature Astronomy, 3, 667Ibata R.A., et al., 2020, arXiv:2012.05245Ji A., et al. 2020, AJ, 160, 181Kaderali S., Hunt J. A. S.; Webb J. J.; Price-Jones N., CarlbergR.” 2019, MNRAS, 484, L114Kamann S., et al., 2020, A&A, 635, 65Kin J.J., Lee Y.W, 2018, ApJ, 869, 35King I., 1966, AJ, 71, 64Krause M., Charbonnel C., Decressin T., Meynet G., PrantzosN., 2013, A&A, 552, A121Kroupa P., 2001, MNRAS, 322, 231Kroupa P., 2020, Star Clusters: From the Milky Way to theEarly Universe. Proceedings of the International Astronomi-cal Union, A.Bragaglia, M.Davies, A.Sills, E.Vesperini editors,Volume 351, pp. 117-121Lane J. M. M., Navarro J. F., Fattahi A., Oman K. A., Bovy, J.,2020, MNRAS, 492, 4164Lardo C., Bellazzini, M.; Pancino, E.; Carretta, E.; Bragaglia, A.;Dalessandro E., 2011, A&A, 525, A114Lee J.-W., 2017, ApJ, 844, 77Lee J.-W., 2020, ApJ, 888, L6Leon S., Meylan G., Combes F., 2000, A&A, 359, 907Li H., Gnedin O.Y., 2019, MNRAS, 486, 4030Libralato M., Bellini A., Piotto G., Nardiello D., van der MarelR., Anderson J., Bedin L. R., Vesperini E., 2019, ApJ, 109Lucatello, S.; Sollima, A.; Gratton, R.; Vesperini, E.; D’Orazi, V.;Carretta, E.; Bragaglia, A., 2015, A&A, 584, A52Mastrobuono-Battisti A., Perets H. B., 2013, ApJ, 779, 85Mastrobuono-Battisti A., Perets H. B., 2016, ApJ, 823, 61 McKenzie M., Bekki K., 2021, MNRAS, 500, 4578Miholics M., Webb J. J., Sills A., 2015, MNRAS, 454, 2166Milone A. P., et al., 2017, MNRAS, 464, 3636Milone A. P., Marino A. F., Mastrobuono-Battisti A., Lagioia E.P., 2018, MNRAS, 479, 5005Milone A. P., et al., 2020, MNRAS, 491, 515Nardiello D., et al. 2015 A&A, 573, 70Pancino E. et al. 2007, ApJ, 661, L155Reino S., Rossi E.M, Sanderson R.E., Sellentin E., HelmiA., Koppelman H.H., Sharma S., 2020, submitted to MN-RAS(arXiv:2007.00356)Renaud F., 2018, NewAR, 81, 1Renaud F., 2020, in Proceedings IAU Symposium 351, BragagliaA., Davies M., Sills A., Vesperini E., eds., pp. 40-46Richer H. B., Heyl J., Anderson J., Kalirai J. S., Shara M. M.,Dotter A.,Fahlman G. G., Rich R. M., 2013, ApJ, 771, L15Simioni M., Milone A. P., Bedin L. R., Aparicio A., Piotto G.,Vesperini E., Hong J., 2016, MNRAS, 463, 449Sollima A., Ferraro F. R., Bellazzini M., Origlia L., Straniero O.,Pancino E., 2007, ApJ, 654, 915Sollima A., Baumgardt H., 2017, MNRAS, 471, 3668Sollima A., 2020, MNRAS, 495, 2222Starkman N., Bovy J., Webb J. J., 2020, MNRAS, 493, 4978Thomas G. F., et al. 2020, ApJ, 902, 89Tiongco M. A., Vesperini E., Varri A. L., 2016, MNRAS, 455,3693Tiongco M. A., Vesperini E., Varri A. L., 2019, MNRAS, 487,5535Trenti M., Vesperini E., Pasquato M., 2010, ApJ, 708, 1598Ventura P., D’Antona F., Mazzitelli I., Gratton R., 2001, ApJ,550, L65Vesperini E., Heggie D.C., 1997, MNRAS, 289, 898Vesperini E., McMillan S. L. W., Portegies Zwart S., 2009, ApJ,698, 615Vesperini E., McMillan S. L. W., D’Antona F., D’Ercole A., 2011,MNRAS, 416, 355Vesperini E., McMillan S. L. W., D’Antona F., D’Ercole A., 2013,MNRAS, 429, 1913Vesperini E., Varri, A. L.; McMillan, S. L. W.; Zepf, S. E., 2014,MNRAS, 443, L79Vesperini E., Hong J., Webb J. J., D’Antona F., D’Ercole A.,2018, MNRAS, 476, 2731Wang L., Kroupa P., Takahashi K., Jerabkova T., 2020, MNRAS,491, 440Webb J. J., Leigh N. W. C., 2015, MNRAS, 453, 3278Webb J. J., Vesperini E., Dalessandro E., Beccari G., Ferraro F.R., Lanzoni B., 2017, MNRAS, 471, 3845Webb J.J., Bovy J., 2019, MNRAS, 485, 5929Whitehead A. J., McMillan S. L. W., Vesperini, E., PortegiesZwart, S., 2013, ApJ, 778, 118Zennaro M., Milone A. P., Marino A. F., Cordoni G., Lagioia, E.P., Tailo, M., 2019, MNRAS, 487, 3239Zonoozi A. H., Kuepper A. H. W., Baumgardt H., Haghi H.,Kroupa P., Hilker M., 2011, MNRAS, 411, 1989This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000