Monte Carlo simulations of multiple populations in globular clusters: constraints on the cooling flow vs. accretion scenario using million bodies simulations
aa r X i v : . [ a s t r o - ph . GA ] F e b MNRAS , 1–18 (2020) Preprint 4 February 2021 Compiled using MNRAS L A TEX style file v3.0
Monte Carlo simulations of multiple populations in globular clusters:constraints on the cooling flow vs. accretion scenario using millionbodies simulations
A. Sollima ⋆ INAF Osservatorio di Astrofisica e Scienza dello spazio di Bologna, via Gobetti 93/3, 40129 Bologna, Italy
Accepted 2021 January 15. Received 2021 January 15; in original form 2020 November 9
ABSTRACT
I simulate the evolution of a stellar system hosting two stellar populations whose initial set up is defined accordingto the two main scenarios proposed for the origin of multiple populations in Galactic globular clusters: (i) formationof a second generation from a cooling flow of pristine+polluted gas and (ii) accretion of polluted gas onto the proto-stellar disks of a fraction of low-mass stars. For this purpose, Monte Carlo simulations containing from 10 up to3 · particles have been run including the effect of stellar evolution, binary interactions, external tidal field and adetailed modelling of the proto-stellar disk structure. The early accretion of gas onto proto-stellar disks is unable toproduce discrete populations and to alter the chemical composition of a significant ( > t disk ∼ M yr ) than that predicted by models is assumed. Moreover, in this scenario themixing timescale of the two populations is too short to reproduce the observed segregation of the chemically enrichedpopulation. On the other hand, simulations run within the cooling flow scenario can evolve after a Hubble time intostellar systems with a first-to-second population mass ratio similar to that observed in globular clusters, provided thatan initial filling-factor r h /r J > .
15 is adopted. However, in the weak tidal field regime a radial segregation of thesecond population stronger than what observed in Milky Way globular clusters at large Galactocentric distances ispredicted. This discrepancy disappears in simulations following eccentric orbits in a realistic axisymmetric potential.
Key words: methods: numerical – stars: kinematics and dynamics – stars: Population II – globular clusters: general
Since the beginning of the new century, the existence ofmultiple populations in Galactic globular clusters (GCs)emerged and revolutionized our view of these stellar systems(Gratton et al. 2019). In particular, the milli-mag precisionof the photometers on-board the Hubble Space Telescope al-lowed to distinguish a series of discrete sequences in vari-ous regions of the colour-magnitude diagrams of almost allknown GCs (Piotto 2006; Piotto et al. 2015; Milone et al.2017, and references therein). Stars in different sequence arecharacterised by different abundance of light elements (He,C, N, O, Mg, Al, Na) often anticorrelated between them,while iron-peak and α -elements appears rather homogeneous(Carretta et al. 2009). This evidence suggests that GCs un-derwent a process of self-enrichment where new stellar pop-ulations originated from gas enriched in p-capture elementsby the evolution of stars of previously formed populations.The stars displaying an evolved chemical composition (withenhanced He, Na, Mg and depleted O and Al, hereafter re-ferred as ”second population”; SP) are generally more con- ⋆ E-mail: [email protected] centrated than the stars with canonical abundance patterns(”first population”; FP) (Lardo et al. 2011; Dalessandro et al.2019).Moreover, the distribution of stars in the Na-O and Mg-Alabundance planes cannot be reproduced by a simple closedbox model where all the material expelled by a given polluteris recycled in purity, regardless of the nature of the polluter.This implies that the polluted gas has been diluted with aconsistent fraction of pristine gas (D’Ercole et al. 2010).The lack of an extended Main Sequence turn-off mor-phology together with the observed abundance pattern(Langer, Hoffman, & Sneden 1993, anomalies requiring p-cycles occurring at high temperatures reachable only in mas-sive stars) suggest a fast enrichment occurring at early stages( t <
Myr ) led by the evolution of relatively massive(
M > M ⊙ ) stars (Renzini 2008). However, while in a stel-lar population characterised by a standard initial mass func-tion (Kroupa 2001) these stars release less than 30% of thetotal mass budget, the SP generally constitutes more than60% of the whole cluster population (Carretta et al. 2009;Milone et al. 2017). To overcome this problem (often referredas ”the mass budget problem”), many hypotheses have beenput forward (Bastian & Lardo 2018), which can be schemat-ically divided in two main scenarios: © Sollima • The cooling flow scenario: after the formation of the FP,the cool gas flowing from the evolution of some massive pol-luter accumulates in the central region of the cluster whereit is mixed with a fraction of pristine gas and forms theSP (Calura et al. 2019). The majority ( > • The accretion onto proto-stellar disks scenario: the FPforms with a strong mass segregation, with the most massivestars located in the central region of the cluster. The pol-luted gas released in the core by these massive stars is thenaccreted into the proto-stellar disks of those low-mass starspassing through the core. Only stars with low orbital ener-gies accrete gas and transform into SP stars, being thereforenaturally settled preferentially in the central cluster region.This scenario, originally proposed by Bastian et al. (2013),has the advantage to require only the small amount of pol-luted gas which adds up to that already available in low-massproto-stars which naturally constitute a reservoir of pristinegas. Useful polluters for this scenario can be FRMS or MIB.In recent years, many groups attempted to test theabove scenarios using chemical evolution models to repro-duce the observed abundance patterns (Lind et al. 2011;D’Ercole et al. 2012; Cassisi & Salaris 2014). Unfortunately,none of the proposed polluters is able to adequately repro-duce the distribution of stars in the Na-O and Mg-Al anticor-relations planes, unless an ad-hoc modification of the stellaryields is made (Bastian, Cabrera-Ziri, & Salaris 2015).From the dynamical point of view, the cooling flow scenario has been widely tested using N-body sim-ulations focused on the derivation of the dynami-cal mixing timescale (Decressin, Baumgardt, & Kroupa2008; Vesperini et al. 2013), on the comparison of themass function (Vesperini et al. 2018), degree of rota-tion (Tiongco, Vesperini, & Varri 2019), binary fraction(Vesperini et al. 2011; Hong et al. 2015, 2016) and relativefrequency (D’Ercole et al. 2008; Khalaj & Baumgardt 2015)of FP/SP stars. All these studies are based on simulationsrun assuming a simplified tidal field and with an initial num-ber of particles ∼ , significantly smaller than the ex-pected number of objects in real proto-GCs. On the otherhand, the predictions of the accretion onto proto-stellar disks scenario have been deduced only through analytical consid-erations in Bastian et al. (2013) and H´enault-Brunet et al.(2015). In both cases, low-mass stars passing through thecore in an isochrone potential were picked as SP regardless ofthe amount of gas that can be accreted by their proto-stellardisks during the passage and without a detailed modelling ofthe proto-stellar disk structure.All the above studies concluded that, under suitable initialconditions, both scenarios can reproduce the present-day gen- eral structural parameters (mass, half-mass radius) of GCs aswell as the main observational evidence regarding the prop-erties (relative fraction and radial segregation) of FP/SP, atleast in a qualitative way.In their analysis, H´enault-Brunet et al. (2015) run N-bodysimulations starting from initial conditions for the FP andSP specifically set to mimic the predictions of both sce-narios. From their study they found that the two scenar-ios predict slightly different residual rotation. However, theN-body simulations by Tiongco, Vesperini, & Varri (2019)showed that, at least for the cooling flow scenario, the finalrotation pattern of the two populations can vary accordingto the amount of mass loss and the dynamical age of the sys-tem (see also Bellini et al. 2015). Similarly, cluster-to-clustervariation in the differential rotation pattern of multiple popu-lations can be naturally produced by other proposed scenarios(e.g. Gieles et al. 2018). The observational evidence collectedso far for differential rotation have indeed shown different be-haviours in various GCs (Pancino et al. 2007; Cordero et al.2017; Cordoni et al. 2020; Bellini et al. 2018).One of the main difficulties in this context is due to theearly epoch during which the formation of multiple popula-tions took place. The process of formation of multiple popu-lation occurs indeed in a short timescale ( < Myr ) dur-ing which radiation, gas, and stars closely interact throughcomplex processes like e.g. SNe feedback, gas-expulsion, stel-lar evolution-driven expansion, gas cooling and mixing, com-petitive accretion, etc. Such a short period is then followedby a long period during which two-body relaxation tends toerase the structural/dynamical differences left by the forma-tion process. All these processes occur with efficiencies andon timescales depending on different powers of the mass andradius of the stellar system, so that a consistent picture ofthe dynamical evolution of GCs can be obtained only withsimulations with characteristics (number of particles, fractionof binaries, mass function) as close as possible to those of realGCs. However, proto-GCs in the early Universe could havebeen as massive as 10 M ⊙ , orders of magnitudes larger thanthe present-day capabilities of N-body simulations.A valuable alternative is provided by Monte Carlo simula-tions (H´enon 1971). In this approach, the integrals of motionof particles are perturbed and updated at odds with N-bodysimulations where the orbits of individual particles need tobe computed. This allows to use a relatively large integrationtime-step which fastens the computing time by more thanone order of magnitude. At the present day, the largest MonteCarlo simulations contain a few 10 particles (Kremer et al.2020) allowing for the first time to simulate the evolutionof GC-like objects from their formation till the present-day.Still, million-bodies Monte Carlo simulations of star clustersaccounting for the formation scenarios of multiple popula-tions are not available yet.In this paper, I present Monte Carlo simulations of GC-likestellar systems whose initial conditions are set according tothe two main scenarios proposed for the formation of multiplepopulations with the aim of investigating if they can evolvetoward a structural/kinematic configuration consistent withthe available observational evidence.In Sect. 2 the Monte Carlo code adopted for the simulationsis presented. In Sect. 3 the initial setup of simulations in theframework of the cooling flow and accretion onto proto-stellardisks scenarios are described. Sect. 4 is devoted to the presen- MNRAS , 1–18 (2020) onte Carlo simulations of GC populations tation of the results and the comparison with observationaldata. I summarize the results of the paper in Sect. 5. The simulations adopted in this paper are run us-ing a modification of the Monte Carlo code de-scribed in Sollima & Mastrobuono Battisti (2014) andSollima & Ferraro (2019). The code is based on the methodoriginally developed by H´enon (1971) and refined by severalauthors (Giersz 1998; Joshi, Rasio, & Portegies Zwart 2000;Vasiliev 2015).Schematically, the cluster is represented by a sample ofstars characterized by their integrals of motion (energy andangular momentum) and their mass. Starting from an ini-tial configuration, at each time-step the following steps areperformed sequentially: • The potential is evaluated according to the masses andpositions of the stars; • Stars are randomly distributed along their orbit accord-ing to their energy and angular momentum; • Stars with a contiguous ranking in distance to the centreare assumed to interact and the perturbations in their E andL are applied.The above steps are iteratively repeated until the end of thesimulation. The time-step is chosen to ensure small perturba-tions at any distance from the centre i.e. it must be smallerthan the local relaxation-time. A detailed description of theabove technique can be found in H´enon (1971).Although this speed has the cost of losing resolution onthose processes occurring on very short timescales and limitthe treatment of the three-dimensional structure of the clus-ter to spherical symmetry, this technique has proven to ac-curately follow the dynamical evolution depicted by N-bodysimulations at low-N in terms of both structural evolutionand kinematics (Giersz et al. 2013).During the last decades, several authors added additionallevels of complexity to this method like the presence of a massspectrum (Giersz 2001; Joshi, Nave, & Rasio 2001), the ef-fects of stellar evolution (Giersz 2001; Chatterjee et al. 2010),direct integration of binary-single and binary-binary inter-actions (Fregeau & Rasio 2007) and interaction with a tidalfield (Sollima & Mastrobuono Battisti 2014).A detailed description of the basic versionof the code adopted here can be found inSollima & Mastrobuono Battisti (2014), where the treat-ment of interactions among single mass particles and theinteraction with different kinds of tidal fields are described.In a subsequent version of the code (Sollima & Ferraro 2019),the presence of a mass spectrum and the direct integrationof single-single, binary-single and binary-binary interactionshave been accounted for.In the version of the code adopted here, two further im-provements have been developed: • Stellar evolution : the evolutionary time has beenextracted from the set of evolutionary tracks ofPietrinferni et al. (2006) with Z=0.001 and [ α /Fe]=+0.4.After this timescale the mass of the star is reduced accordingto the relations of Kruijssen (2009). A kick extracted froma Maxwellian distribution and randomly oriented has been added to the remnant of massive ( M > M ⊙ ) stars. The1D dispersion of the kick velocity has been assumed to be σ kick = 100 km/s for neutron stars (originating from the evo-lution of stars with 8 < M/M ⊙ <
30) and σ kick = 80 km/s for black holes (from stars with M > M ⊙ ). When anevolving star is part of a binary system, the mass-loss andthe kick are applied, the orbit of the system is followedfor 10 orbital periods and, if not broken, its new equilib-rium characteristics (semi-major axis and eccentricity) areupdated. • Three-bodies encounters : in the densest cluster region itis possible that three objects simultaneously cross the samesmall volume for a short interval of time. It can happen thattwo objects transfer kinetic energy to the third, emerging as abound system. This process, while generally not efficient, canbecome effective in dense regions and when a large mass con-trast among the involved stars is present and represents oneof the main channels of binary formation involving massiveremnants (Giersz 1998; Morscher et al. 2013).At each iteration, for a given pair of neighbour particlesthe closest object is chosen. The 3D velocities of the threeobjects have been chosen according to their values of E andL and randomly oriented. The probability that, during thetime-step ∆ t , the three objects cross a spherical volume ofradius r m in the same crossing time interval is P = π s r m G ( m + m ) n w w (1 + β )(1 + β )∆ t where m , m and m are the masses of the three stars, n is the number density w is the relative velocity of the firsttwo stars and w that of the third star with respect to thecenter of mass of the first two and β = 2 G ( m + m ) r m w β = 2 G ( m + m + m ) r m w I adopted a maximum impact parameter r m = G ( m + m ) E + E − Φ − Φ where E i and Φ i are the energy of the considered star andthe potential at its position. A random number uniformly dis-tributed between 0 and 1 is extracted and, if smaller than theassociated probability, the close interaction is integrated us-ing a symplectic integrator with adaptive timestep (Yoshida1990). For this purpose, the impact parameters b i (i=12,3)have been chosen as b i = r m β i "s η i (cid:18) β i + 1 β i (cid:19) − where η i are random numbers uniformly distributed between0 and 1. The initial positions of the three stars are definedsuch that they simultaneously reach the minimum distancecorresponding to the extracted impact parameters. The inter-action is followed until the total energy of the system exceedszero. The velocities of the emerging objects are then con-verted in the cluster reference frame and their correspondingvalues of E and L are assigned accordingly. MNRAS , 1–18 (2020)
Sollima
In this section, the initial setup of the Monte Carlo simula-tions presented in this paper is described. Of course, some ofthe initial conditions vary according to the adopted scenarioand are discussed in the following subsections.Unless specified otherwise, in both scenarios simulationswere run with an initial FP mass of M F = 10 M ⊙ , a frac-tion of binaries f b = 5%, and a Kroupa (2001) Initial MassFunction (IMF) defined between 0.1 and 120 M ⊙ . The char-acteristics of binary stars (semi-major axis, eccentricity, massof the components) were chosen using the technique descirbedin Sollima, Bellazzini, & Lee (2012). In particular, the massesof the components are chosen by random pairing stars froma Kroupa (2001) IMF, and the period and eccentricity havebeen extracted from the distribution of Duquennoy & Mayor(1991). From an initial library, stars with semi-major axissmaller than the stable overflow criterion by Lee & Nelson(1988) or larger than the ionization limit of Hut & Bahcall(1983) are removed. In this configuration, the total numberof particles at the beginning of the simulation is ∼ . · .To investigate the effect of the initial mass on the long-termevolution of the two populations, a few simulations were runwith a lower (10 M ⊙ ) or higher (2 · M ⊙ ) initial mass.The simulated clusters move in the tidal field generated bya point mass with M g = v circ R g /G , where v circ = 220 km/s ,and follow a circular orbit at galactocentric distance R g =14 . kpc . Such a distance has been chosen to mimic the av-erage tidal field felt by Galactic GCs. For this purpose, theorbits of Galactic GCs have been integrated within the Galac-tic potential of Johnston, Spergel, & Hernquist (1995) usingthe 3D positions and velocities listed by Baumgardt et al.(2019), and the orbit-averaged Jacobi radius at the conven-tional mass of M = 10 M ⊙ ( h r J i ) has been calculated byaveraging along the orbit the instantaneous Jacobi radius ( r J )calculated using eq. A2 of Allen, Moreno, & Pichardo (2006).Galactic GCs have values of 60 < h r J i /pc <
120 with anaverage of h r J i = 85 pc . The same value of h r J i can be ob-tained in the simplified point-mass potential assuming R g =14 . kpc . Of course, this crude approximation provides a guessvalue of R g to mimic the tidal cut imposed by the Galactictidal field but cannot account for the tidal shocks occurring ina realistic Galactic potential as a result of disk crossing andperi-Galactic passages (Ostriker, Spitzer, & Chevalier 1972;Aguilar, Hut, & Ostriker 1988). To check the effect of theadopted tidal field, some simulations have been run adopt-ing a smaller Galactocentric distance ( R g = 5 . kpc ) orthe tidal field generated by the axisymmetric potential ofJohnston, Spergel, & Hernquist (1995) and the orbit of theGC NGC1851 (see Fig. 1). This last orbit, although having arelatively large orbit-averaged Jacobi radius h r J i = 109 pc ,is extremely eccentric ( e = 0 .
9) with a peri-Galactic distanceof ∼ ∼
450 Myr, so that some 60bulge- and 120 disk- shocking occur during the entire clusterevolution.The simulations are run for 12 Gyr allowing to follow theevolution of the relative fraction and of the various structuralproperties of FP and SP.The initial and final properties of the entire set of simula-tions are listed in Table 1.
Figure 1.
Adopted NGC1851-like orbit in the X-Y (left panel), X-Z(middle panel) and R-Z (right panel) planes.
The simulations have been preformed in the framework of thescenario proposed by D’Ercole et al. (2008). In this model,the self-enrichment occurs in the central region of the clus-ter where the gas expelled by intermediate-mass AGB starscools and mix with a fraction of pristine gas, residual fromthe formation of the FP. According to D’Ercole et al. (2010),assuming ad-hoc yields for AGB and super-AGB stars, it ispossible to reproduce the distribution of stars in the Na-Oanticorrelation plane assuming a significant contribution ofpristine gas such that the emerging SP has a mass ∼
10% ofthat of the FP. The formation of the SP occurs after ∼
30 Myrwhen massive AGB stars start to evolve and ends after ∼ <
100 Myr in allthe considered scenarios, the detailed timing of this processshould not have a significant effect on the final outcome ofthe simulation after a Hubble time.According to the above picture, FP stars were extractedfrom the adopted IMF and distributed according an isotropicKing (1966) (Gunn & Griffin 1979, with W = 25 for simula-tions assuming primordial mass segregation) density profile.In one case a FP with a large degree of radial anisotropy( r a /r h,F = 1; where r a is the anisotropy radius beyondwhich stellar orbits progressively start to be radially biased;Gunn & Griffin 1979) has been simulated. The half-mass ra-dius of the FP has been chosen to ensure the large initialRoche-lobe filling factor (0 . < r h /r J < .
2) needed to en-sure a quick loss of FP stars . For reference, a King (1966) model with W = 5 and r h /r J ∼ . , 1–18 (2020) onte Carlo simulations of GC populations Between 30 Myr and 100 Myr after the beginning of thesimulation, SP stars are continuously added in the centralregion of the cluster following a King (1966) profile with cen-tral adimensional potential W = 5 and a half-mass radiuscorresponding to 10% of that of the FP ( r h,S = 0 . r h,F in allsimulations but two where different value of r h,S /r h,F = 0 . r a = r h,S . Inone case the degree of anisotropy of FP and SP have beeninverted, with a radially anisotropic FP and an isotropic SP(see above). The SP is constituted by stars extracted froma Kroupa (2001) IMF between 0.1 and 8 M ⊙ (to avoid theexplosion of SP SNe II; D’Ercole et al. 2008) and containsthe same fraction of binaries with the same characteristicsof those in the FP. The rate at which SP stars are added is dM/dt = 0 . M F Myr − , so that at the end of the SPformation episode (after 100 Myr from the beginning of thesimulation) M S = 0 . M F . Note that, given the large uncer-tainties on the AGB and super-AGB yields (Ventura et al.2013; Bastian, Cabrera-Ziri, & Salaris 2015), as well as thepoor knowledge of the origin of the dilution process, the con-straint on the initial fraction of SP/FP is extremely weak, andother fractions have been proposed (e.g. Cabrera-Ziri et al.2015). Of course, a smaller (larger) initial SP/FP fractionwould require a more (less) efficient loss of FP stars i.e. alarger (smaller) initial Roche-lobe filling factor. Given the de-generacy between these two parameters and considering theuncertainties in the details of the dilution process, it is im-possible to place any constraint on the initial SP/FP fraction.For this reason, I kept fixed this parameter. In the scenario proposed by Bastian et al. (2013), there isonly one star formation episode where the FP forms in amass segregated configuration. After a few Myr massive ( > M ⊙ ) stars evolve and release part of their gas in the centralregion of the cluster. Low-mass stars, still in their pre-MainSequence phase, cross this region and accrete this gas ontotheir proto-stellar disks, mixing it with their own gas. Toaccount for the abundance variation of elements involved inthe nuclear burning (e.g. He), it is necessary that the accretedmaterial is brought from the surface to the stellar core, anoccurrence possible only until the star is entirely convective(D’Antona et al. 2014).To calculate the amount of gas accreted by low-mass stars,I performed an analytical modeling of the disk structure.The initial population is distributed in the phase-spaceaccording to a multi-mass models of Gunn & Griffin (1979)with W = 25. For this purpose, the simulated stellar pop-ulation has been grouped in 29 mass bins, where stars withmasses M > M ⊙ are in the last 11 bins.The gas cloud density profile has been derived as ρ g = X i =19 µ i n i m i where m i is the mass, µ i is the fraction of released gas (from Figure 2.
Proto-stellar disk lifetime ( t disk ; from Li & Xiao 2016,blue line) and duration of the fully convective phase ( t conv ;from Tognelli, Prada Moroni, & Degl’Innocenti 2011, red line) asa function of the stellar mass. Kruijssen 2009), n i is the number density of the stars in thei-th mass bin.The orbits of less-massive ( M < M ⊙ ) stars within thecluster potential are followed using a 4th-order Runge-Kuttaintegrator for the time their proto-stellar disks are able toaccrete gas ( t end ). This timescale has been calculated as theminimum between the disk lifetime ( t disk ) and the time dur-ing which the star is entirely convective ( t conv ). I adoptedthe disk lifetimes calculated by Li & Xiao (2016) (LX16) asa function of the stellar mass. According to these authors,the proto-stellar disk survives longer in massive stars thanin low-mass ones, ranging from 0.8 Myr in a 0.4 M ⊙ starto 8.5 Myr in a 2 M ⊙ star. Note that these models arebuilt assuming the disk formed by gas in stable rotation.This condition can be altered in the case of a strong ac-cretion of gas at low angular momentum, so that the sta-bility of the disk is not guaranteed (Wijnen et al. 2016).The time a star maintain a fully convective structure hasbeen taken from the pre-Main Sequence evolutionary tracksof Tognelli, Prada Moroni, & Degl’Innocenti (2011) with ametal content Z=0.001. The dependence of t disk and t conv as a function of the stellar mass is shown in Fig. 2. It canbe seen that the time during which the proto-star is able toaccrete gas and use it as a fuel in its core is a peaked functionof mass with a maximum duration of ∼ M ∼ M ⊙ . As hypothesized by Bastian et al. (2013), the gasaccretion can prolong the disk lifetime. To account for this ef-fect, the orbit integration is interrupted when the integrationtime exceeds the instantaneous value of t end ( M ( t )). More-over, to account for the uncertain modeling of disk lifetimes,additional simulations were performed assuming an increasedlifetime t disk = 20 Myr regardless of the stellar mass.The gas crossing the disk during this time interval (i.e. themaximum amount of gas that can be accreted by the star) is
MNRAS , 1–18 (2020)
Sollima
Figure 3.
Orbit of a sample star in simulation R14.4lm6rh1 (left panel). The location of the gas is marked by the red area. The correspondingevolution of the distance from the cluster centre (black), the proto-stellar disk truncation radius (blue) and the accretion rate (red) areshown in the right panel.
Figure 4.
Distribution of He mass fration in the simulations runwithin the accretion onto proto-stellar disks scenario. All his-tograms are normalized to the peak value. given by∆ M = Z t end π r d ρ g v cos α dt where v is the modulus of the star velocity, r d is the proto-stellar disk radius and α is the angle between the star posi- tion and velocity vectors. All the above quantities are time-dependent and are calculated along the star orbit.In case of isolated stars, the models ofNakamoto & Nakagawa (1994) and Li & Xiao (2016)indicate that the disk radius grows with time until a gravita-tional instability occurs or X-ray photoevaporation destroythe outer regions of the disk (Owen, Clarke, & Ercolano2012), at distances as large as 50 AU. This is confirmedby the existence of extremely extended disks (up to 1000AU; Lada et al. 2000) around young stars in the Trapeziumcluster. However, within a dense cluster the main truncationmechanism for disks is the close passage of nearby stars. Toaccount for this effect I applied a Monte Carlo scheme toimplement the technique developed by de Juan Ovelar et al.(2012) locally along the orbit followed by the star inside thecluster. At each point along the orbit, 10 synthetic particleshave been extracted with the local mass function and thecorresponding velocity distribution. An impact parameter b = b max √ η has been extracted, with η being a random number uniformlydistributed between 0 and 1, n is the local number densityand b max = (cid:18) π n (cid:19) / The truncation radius implied by the encounter with a starcolliding with the extracted mass m ′ , impact parameter b ,and a relative velocity with modulus w , is given by r ′ tr ( w, m ′ , b ) = G ( m + m ′ ) w (1 + p m ′ /m ) s b w G ( m + m ′ ) − ! The total number of collisions occurring with impact param-
MNRAS , 1–18 (2020) onte Carlo simulations of GC populations Figure 5.
Initial density (left panel) profiles and mass functions (right panel) of FP (blue lines) and SP (red lines) in the R14.4lm6rh1(solid lines) and R14.4lm6rh1td20 (dashed lines) simulations. The gas density profile is overplotted in the left panel with a black line. eter b < b max before the time t is instead given by N coll = 0 . πb max n h w i t where h w i is averaged over the 10 synthetic particles andthe coefficient 0.3 is a corrective factor accounting for themild effect of encounters occurring with inclination angles > ◦ (de Juan Ovelar et al. 2012). The individual values of r ′ tr are sorted and a relation between the ranking index and r ′ tr is derived. The instantaneous truncation radius ( t tr,t ),set by the most likely closest collision, is then defined byinterpolating through the above relation at the index i ′ =10 /N coll .The disk radius at the given time and in a given point alongthe star orbit is finally given by r d = min " (cid:18) t · yr (cid:19) pc, r tr,t To compute the actual enhancement provided by the accre-tion process, the total He mass fraction has been calculatedassuming a cosmological He abundance for the FP (Y=0.245;Planck Collaboration et al. 2016) and the extreme He abun-dance predicted by de Mink et al. (2009) Y = 0 . M + 0 .
44 ∆ MM + ∆ M where M is the original stellar mass before accretion.The evolution of the accretion rate of a sample star in theR14.4lm6rh1td20 simulation during its first Myr is shown inFig. 3. This star, after moving for 20 Myr within the innerpc, reaches a final mass of 0.14 M ⊙ and a He mass fraction Y = 0 .
3. It can be noted that the maximum accretion occurs during pericenters, where the gas density and the star speedare high. The accretion rate is however modulated accordingto the size of the proto-stellar disk which increases at thebeginning of the disk evolution, reaches a maximum and thenprogressively declines as a result of the cumulative effect ofclose collisions. Consequently, most of the accretion occurs atearly epochs, while after a few Myr the disk is almost unableto accrete gas anymore.The distributions of He mass fractions among low-mass(
M < . M ⊙ ; i.e. those expected to survive after a Hubbletime) for the four considered initial conditions are shown inFig. 4. It is immediately apparent that in all simulations, thedistribution is unimodal with a peak at the cosmological value( Y = 0 . Y > . M S /M tot < ∼ > MNRAS , 1–18 (2020)
Sollima
The initial density profiles of gas, FP and SP for the simula-tion R14.4lm6rh1 and R14.4lm6rh1td20 are shown in the leftpanel of Fig. 5. It can be seen that, in both simulations, theSP naturally forms more concentrated than FP, and followsa density profile similar to that of the polluted gas. However,only a small fraction of the gas is used to form the SP. An-other peculiarity of this scenario is related to the resultingmass functions of the two populations (see the right panel ofFig. 5). Indeed, while the FP forms with the standard Kroupa(2001) IMF, the SP is able to accrete efficiently gas only ina restricted range of mass, thus following a mass functionpeaked at log ( M/M ⊙ ∼ − .
3, with a steep slope α ∼ − − . < log ( M/M ⊙ ) < . − . < log ( M/M ⊙ ) < .
2. Qualita-tively, the same evidence is apparent in simulations where thedisk lifetime has been increased to 20 Myr at all masses. Inthis case, low-mass stars have more time to accrete gas, sothat the SP is more abundant and the mass function peak isshifted at a lower mass.
In the top panel of Fig. 6 the fractions of SP stars in thelast snapshot of the 10 simulations run within the cool-ing flow scenario are plotted against the orbit-averaged Ja-cobi radius (i.e. indicating the strength of the tidal fieldthey are subject to) and compared with the fractions mea-sured in Galactic GCs by Milone et al. (2017). Among sim-ulations moving on a circular orbit, it can be noted thatboth simulations starting without primordial mass segrega-tion (R14.4lm6w5rh36dr01 and R14.4lm6w3rh36dr01) pre-dict a fraction of SP stars N S /N tot < .
4, significantly smallerthan that observed in real GCs. Between them, the simulationwith a shallower density profile (R14.4lm6w3rh36dr01) losesslightly more FP stars than the one with the steeper profile(R14.4lm6w5rh36dr01). This is expected since in the formersimulation a larger fraction of FP stars is located at largedistances from the cluster centre with respect to the lattersimulation which is more concentrated and therefore moreresistant to evaporation (see also Vesperini 1997). Amongsimulations with primordial mass segregation, the fractionof SP stars presents a strong dependence on the Roche-lobe filling factor, with the most tidally-filling simulations( r h /r J ≥ .
15) able to reach the large observational values.Both the above evidences are expected as a result of the in-creasing amount of mass loss (mainly FP stars) experiencedby simulations starting with primordial mass segregation andwith a large filling factor. Indeed, the potential energy lost inthe cluster centre is given by ∆Φ = P G ∆ m i /r i , where ∆ m i is the mass lost by the i-th star because of stellar evolutionand r i is the distance of the star from the cluster centre. Inmass segregated simulations massive stars (subject to a largeamount of mass loss) are preferentially located at small dis-tances, so that the loss of potential energy is larger than innon-segregated simulations. These simulations strongly ex-pand to maintain the virial equilibrium pushing FP stars,preferentially located in the peripheral cluster region, out-side the Jacobi radius. In the same way, the larger is the ini-tial Roche-lobe filling factor the larger is the fraction of stars Figure 6.
Fraction of SP (top panel) and A +2 parameter (bottompanel) measured in the last snapshot of all the simulations runwithin the cooling flow scenario. Grey points represents the obser-vational fractions for the Galactic GCs (from Milone et al. 2017;Dalessandro et al. 2019). The colour code is described in the lastcolumn of Table 1. Circles and triangles represent simulations withand without primordial mass segregation. The symbol size is pro-portional to the initial Roche-lobe filling factor r h /r J . exceeding the Jacobi radius during this early expansion. In-stead, SP stars segregated in the cluster centre remain at thebottom of the potential well and are more efficiently retaineduntil the end of the simulation. The subsequent evolution hasonly a minor effect in the global fraction of SP stars. This pro-cess is depicted in the top panel of Fig. 7, where the evolutionof the number ratio of FP and SP stars of the 10 simulationsare compared.To investigate the degree of mass segregation of SP stars,the parameter A +2 (Dalessandro et al. 2019) has been calcu-lated. For this purpose, the 3D positions of particles withmasses m > . M ⊙ have been projected on the plane ofthe sky assuming an isotropic spherical symmetry. The cu- MNRAS , 1–18 (2020) onte Carlo simulations of GC populations Figure 7.
Evolution of the SP number fraction (top panel) andSP/FP half-mass radii (bottom panel) for all the simulations runwithin the cooling flow scenario. The colour code is described inthe last column of Table 1. mulative distributions of FP/SP stars located within 2 r h from the cluster centre, have been computed and normal-ized to their respective total numbers. The A +2 parameteris then calculated as the difference of the integrals of thetwo distributions. In the bottom panels of Fig. 6 the valuesof A +2 computed on the last snapshot of the 10 simulationsare plotted against h r J i and compared with those estimatedby Dalessandro et al. (2019) for 18 Galactic GCs. Note that,while real GCs show values of − . < A +2 < .
08 with aweak decreasing trend with h r J i , at the end of all the simu-lations moving on circular orbits at a distance R g = 14 . kpc the SP is much more segregated in the centre ( A +2 < − . Figure 8.
Evolution of all the simulations run within the coolingflow scenario in the M − r h plane. The large dots mark the end-points of the simulations. Grey points represents the Galactic GCswith 70 < h r J i /pc < long relaxation time of these simulations which reduces theefficiency of the dynamical mixing between the two popula-tions. The half-mass relaxation time (Spitzer 1987) is indeed t rh > Gyr at the beginning of all these simulations andremains longer than 3 Gyr during their entire evolution. Notethat the dynamical mixing efficiency of the two populationsdepends on the tidal truncation and not on the initial half-mass radius. Indeed, the same discrepancy is present in theunderfilling simulations R14.4lm6w25rh18dr01ms which, inspite of its relatively short initial half-mass relaxation time,quickly expand without loosing a significant amount of mass,thus increasing their half-mass relaxation time during theirevolution. A better agreement is instead found with simula-tion R5.1lm5w25rh18dr01ms starting with a small half-massradius while being tidally-filling. In this last case, the tidal cutkeeps the simulation compact with a relaxation time whichis short enough to allow the two populations to efficientlymix. It is also interesting to note that, comparing simu-lations R14.4lm6w25rh36dr005ms, R14.4lm6w25rh36dr01msand R14.4lm6w25rh36dr02ms, the initial concentration of theSP has only a minor effect on both the SP fraction and on itsradial segregation. In fact the SP/FP ratio mainly dependson the efficiency of evaporation of the FP, while the SP is al-most unaffected by mass loss. Moreover the SP, when initiallymore concentrated, expands faster and partially compensatesfor the initial condition.In Fig. 8 the evolution of the 10 simulations in the massvs. half-mass radius plane is shown. For comparison, the loca-tion of Galactic GCs with 70 < h r J i /pc <
100 is shown (seeSect. 3). The evolution of all simulations follows in this planea characteristic path: an initial expansion (due to the stellarevolution-driven mass loss) followed by a temporary contrac-
MNRAS , 1–18 (2020) Sollima
Figure 9.
Mass function of FP (blue line) and SP (red line) measured in the last snapshot of all the simulations run within the coolingflow scenario.MNRAS , 1–18 (2020) onte Carlo simulations of GC populations Figure 10.
Anisotropy parameter of FP (blue lines) and SP (red lines) as a function of the distance from the cluster centre measured in thelast snapshot of the simulations run according to the cooling flow scenario. Distances are normalized to the final global half-mass radius.MNRAS , 1–18 (2020) Sollima tion of the half-mass (due to the formation of the SP in thecentre), and a subsequent phase where both mass and half-mass radius decrease. The length of the path depends mainlyon the initial Roche-lobe filling factor r h /r J , with tidally-filling clusters evolving toward a less massive and more con-centrated structure than underfilling ones. All the simulationstarting with a Roche-lobe filling factor r h /r J > .
15 loose >
90% of their initial mass and reach after 12 Gyr massessignificantly smaller (0 . < ∆ log ( M/M ⊙ ) < .
8) than themean present-day mass of GCs ( h log ( M/M ⊙ ) i ∼ . cooling flow scenario considering that h r J i varies by more than an order of magnitude among theGalactic GC system.A solution is provided by the simulations moving in theaxisymmetric potential of Johnston, Spergel, & Hernquist(1995) and following the eccentric orbit of NGC1851(N1851lm6.3w25rh34dr01 and N1851lm6.3w25rh17dr01ms).These simulations, although starting in underfilling condi-tions, suffer strong mass loss during the peri-Galactic pas-sages and the disk crossing. During the simulation, thehalf-mass radius is set at a reduced size with respect toa simulation moving on a circular orbit, being character-ized by a shorter half-mass relaxation time. Consequently,both the fraction of SP stars and the A +2 turn out to beconsistent with the observational values, in spite of the ap-parently large value of h r J i . Moreover, in these simula-tions primordial mass segregation is not a necessary condi-tion: also the simulation without primordial mass segregation(N1851lm6w25rh34dr01) can indeed reproduce the fraction ofSP.In Fig. 9 the mass functions of FP and SP measured in thelast snapshot of all simulations are shown. The mass functionsof both populations are very similar, with slopes comprisedin the range − . < α < − . α < .
15 (i.e. with the FP with a slightly flatter MF thanthe SP), in agreement with what found by Vesperini et al.(2018). As expected, because of the preferential loss of low-mass stars (Baumgardt & Makino 2003), simulations loosinga large fraction of their initial mass are those with the flattestmass function.The anisotropy parameter β ≡ − σ t /σ r (where σ r and σ t are the dispersions of the radial and tangential com-ponent of the velocity) is plotted as a function of thedistance from the cluster centre in Fig.10 for all the 10simulations. Apart from small differences, in all the sim-ulations the FP shows an anisotropy profile declining atlarge radii toward the regime of tangential anisotropy ( β < Figure 11.
Same as Fig. 6 but for the simulations run according tothe accretion onto proto-stellar disks scenario. nal radial anisotropy. The results presented here agree withthose of H´enault-Brunet et al. (2015). The simulations ofTiongco, Vesperini, & Varri (2019) also predict a similar be-haviour in the early stages of evolution while after many re-laxation times they predict isotropic profiles for both FP andSP. Given the extremely long relaxation time of the simula-tions presented here, the predictions of these works can beconsidered in agreement.
In Fig. 11 the fraction of SP and the A +2 value measured unthe last snapshot of the 6 simulations run according to the accretion onto proto-stellar disks scenario. All the 6 simula-tions lead to very small SP fraction ( N S /N tot < . MNRAS , 1–18 (2020) onte Carlo simulations of GC populations Figure 12.
Same as Fig. 7 but for the simulations run according tothe accretion onto proto-stellar disks scenario. to reproduce the present-day half-mass radius of GCs. More-over, the SP is only mildly more concentrated than FP atthe beginning of these simulations. So, both FP and SP arefree to expand and interact during the subsequent evolutionwithout a preferential loss of FP stars, thus leaving the frac-tion of SP stars almost unchanged. This is apparent in thetop panel of Fig. 12 where the evolution of the SP fraction isshown for all the 6 simulations.Regarding the radial distribution of SP, all the simulationswith an increased disk lifetime have SP completely mixedwith FP. As shown in the bottom panel of Fig. 12, in com-pact and underfilling conditions, two-body relaxation is effi-cient in mixing the two populations at early epochs, while thedynamical evolution at later times ( t > Gyr ; after the stel-lar evolution-driven expansion) only marginally affects theradial segregation of SP. In this respect, simulations with ex-tended disk lifetimes host a SP containing a larger fraction oflow-mass stars (see Fig. 5) and a more extended radial dis-
Figure 13.
Same as Fig. 8 but for the simulations run according tothe accretion onto proto-stellar disks scenario. tribution with respect to simulations with the disk lifetimespredicted by the LX16 model.In Fig. 13 the evolution of the 6 considered simulations inthe mass vs. half-mass radius plane is shown. Most simula-tions show the same behaviour with an expansion (driven bystellar evolution) and a modest mass loss (∆ log ( M/M ⊙ ) ∼ .
3; almost entirely addressed to stellar evolution mass loss).In this case, simulations with an initial mass of M = 10 M ⊙ and half-mass radius r h = 1 pc well reproduce after 12 Gyrthe average present-day mass and half-mass radius of Galac-tic GCs.Note also that in this scenario the simulations run in a sim-plified potential (R14.4lm6rh1td20 and R5.1lm6rh1td20) andthat run in the axisymmetric potential (N1851lm6rh1td20)produce very similar results. Indeed, in all these simulations,the entire cluster is always contained in a small volume wellinside the Jacobi radius.In Fig. 14 the mass functions of the two populations atthe end of the 6 simulations are compared. The strong massfunction difference observed at the beginning of the simula-tion (Sect. 3.2) remains apparent also after 12 Gyr in all theconsidered simulations. This is not surprising given the rel-atively small amount of mass lost by both populations. Inparticular, while the FP has a mass function consistent witha power-law, the mass function of the SP is significantly de-pleted of low-mass ( M < . M ⊙ ) stars.The anisotropy profiles measured in the last snapshot ofthe 6 simulations are shown in Fig. 15. In all simulations,both populations show a radially anisotropic profile. This isa consequence of the many close encounters occurring in thedense cluster centres which eject stars of both populations onradial orbits at large distances (Lynden-Bell & Wood 1968;Spitzer & Shull 1975). The behaviour of the anisotropy pro-files shown in Fig. 15 differs from that predicted by the sim-ulations of H´enault-Brunet et al. (2015). In the simulations MNRAS , 1–18 (2020) Sollima
Figure 14.
Same as Fig. 9 but for the simulations run according to the accretion onto proto-stellar disks scenario. of these authors indeed the FP is isotropic at the end ofthe evolution, similarly to what happens in the cooling flow scenario, while the SP is radially anisotropic at large radii.The difference between the two works likely resides in thedifferent initial structure of the simulations. The simulationsof H´enault-Brunet et al. (2015) are indeed ∼ I presented the results of a large set of Monte Carlo sim-ulations of star clusters run adopting the initial conditionspredicted by the two main scenarios proposed so far for theformation of multiple populations. These simulations, con-taining 1 . ÷ · particles, are among the largest simu-lations ever run and allow for the first time to reach finalmasses similar or only a factor ∼ accretion onto proto-stellar disks scenario, the ac-cretion rate of individual disks has been modeled for the firsttime, accounting for the evolution of the size of the proto-stellar disk. Because of the short disk-lifetimes, only a negli-gible fraction of polluted material can be accreted onto theproto-stellar disks of low-mass stars. Even adopting a mass- independent disk lifetime of 20 Myr, the SP do not exceed10% of the whole stellar content. Moreover, the continuousdistribution of orbital energies leads to a continuous distribu-tion of the time spent by stars within the polluted gas cloud,thus avoiding the formation of discrete populations. Becauseof the small disk radii at low masses and the quick disappear-ance of the convective proto-stellar core, the IMF of SP is de-pleted of low- ( M < . M ⊙ ) and high-mass ( M > . M ⊙ )stars. The subsequent dynamical evolution does not signif-icantly change any of the above initial differences. Indeed,the initial density required to constitute a significant amountof polluted gas implies an initially concentrated structure( r h (0) = 1 pc ), so that the cluster remains extremely re-sistant to tidal effects and retains a large fraction of stars ofboth populations. The final fraction of SP turns out to betherefore inconsistent with those observed today in GalacticGCs (Milone et al. 2017). Several works foresaw the difficultyof this scenario in producing discrete populations (see e.g.Renzini 2008) and in maintaining an efficient convective corefor a significant amount of time (D’Antona et al. 2014), al-though this is the first time these problems are quantified.These issues adds to the already known problems of this sce-nario in reproducing the Li abundance of SP (D’Antona et al.2014) and the stability of disks in a strong-accretion regime(Wijnen et al. 2016).The simulations run under the cooling flow scenario, given MNRAS , 1–18 (2020) onte Carlo simulations of GC populations Figure 15.
Same as Fig. 10 but for the simulations run according to the accretion onto proto-stellar disks scenario. the large number of involved free parameters, are instead ableto reproduce at the end of their evolution the present-daygeneral properties of FP/SP (in terms of their relative frac-tion and radial segregation) reaching a global structure com-parable (i.e. a similar half-mass radius and a slightly lowermass) to those of Galactic GCs. The imperative conditionsto provide a good agreement are however that the cluster i) form in a tidally-filling configuration ( r h /r J > .
15; al-though this condition depends on the initial SP/FP massratio) and ii) maintains a compact structure during its evo-lution. The first condition ensures a large fraction of SP stars,while the second allows the spatial mixing of the two popu-lations. Indeed, the larger is the filling-factor, the more ef-ficient is the loss of FP determining the final proportion ofSP/FP. On the other hand, the smaller is the cluster sizethe shorter is the relaxation time, thus increasing the ef-ficiency of two-body relaxation. The lack of a decreasingtrend between the relative fraction of SP with the Galac-tocentric distance implies that the above conditions mustbe fulfilled at all Galactocentric distances. This can repre-sent a big problem for the cooling flow scenario. Indeed,in a logarithmic potential (like the one of the Milky Wayhalo; Johnston, Spergel, & Hernquist 1995), the Jacobi ra-dius grows as r J ∝ R / g . So, to keep the filling-factor fixed,the initial half-mass radius of the cluster should grow ac- cordingly. It is therefore impossible to satisfy both the aboveconditions across the large range of Galactocentric distances(0 . < R g /kpc <
125 Harris 1996, 2010 edition) coveredby Galactic GCs. In particular, at large Galactocentric dis-tances GCs the Galactic tidal field is so weak that the clusterspend all its life with a large half-mass radius constitutinga non-collisional environments where the SP would maintaina strong radial segregation (reaching r h,S /r h,F ∼ . MNRAS , 1–18 (2020) Sollima ever, it is known that GCs orbits become radially anisotropicin the outer halo (Vasiliev 2019), and this can contribute tomaintain a significant tidal truncation also for GCs at largeGalactocentric distances. Unfortunately, only a single eccen-tric orbit has been considered here, so that it is not clearif such an effect can solve the above problem for the GCsalong the entire range of Galactocentric distances. Furthersurveys of simulations, tailored to reproduce the observedSP/FP population ratios and following the actual orbit of asignificant number of GCs are necessary to clarify this issue.It is interesting to note that the two considered scenar-ios foresee opposite predictions for the mass functions andthe anisotropy profiles of FP and SP. The determination ofthese quantities can therefore be good tools to distinguishbetween these two scenarios. In particular, while in the cool-ing flow scenario no significant differences are expected inthe mass function of FP and SP, in the accretion onto proto-stellar disks scenario the SP should show a significant de-pletion of low-mass (
M < . M ⊙ ) masses. Unfortunately,the only available study of this kind (Milone et al. 2012) isfocused on the very central region of the cluster where theeffect of mass segregation produces strong modifications tothe global mass function. Similarly, while in the cooling flow scenario the FP is characterized by a significantly smaller de-gree of radial anisotropy than the SP, in the accretion ontoproto-stellar disks scenario the two populations are expectedto share the same anisotropic profile. However, observationalstudies provide conflicting results, with some clusters hostingstellar populations with similar anisotropy profiles and otherswith significant differences (Cordoni et al. 2020).In the present paper, I do not compare the predicted be-haviour of other kinematic parameters of FP and SP likee.g. binary and remnant fraction. While these quantitiescontain important information, their present-day appearancestrongly depend on their uncertain initial conditions. Futurework will be addressed to investigate the evolution of thesequantities spanning a wide range in their initial conditions. ACKNOWLEDGMENTS
I warmly thank Enrico Vesperini for useful discussions. I alsothank the anonymous referee for his/her helpful commentsand suggestions.
DATA AVAILABILITY
The data underlying this article will be shared on reasonablerequest to the corresponding author.
MNRAS , 1–18 (2020) o n t e C a r l o s i m u l a t i o n s o f G C pop u l a t i o n s Table 1.
Initial and final properties of the performed simulations. Initial Finalname R g Mass W M F r h M S /M tot r h,S /r h,F r h /r J M r h N S /N tot A +2 colour t disk segregation /Notekpc 10 M ⊙ pc 10 M ⊙ pc Myr Cooling flow scenarioR14.4lm6w5rh36dr01 14.4 no 5 10 36 0.091 0.10 0.200 2.21 22.68 0.235 -0.546 redR14.4lm6w3rh36dr01 14.4 no 3 10 36 0.091 0.10 0.200 1.57 21.66 0.329 -0.585 magentaR14.4lm6w25rh36dr01ms 14.4 yes 25 10 36 0.091 0.10 0.200 0.57 9.08 0.812 -0.336 blackR14.4lm6w25rh18dr01ms 14.4 yes 25 10 18 0.091 0.10 0.100 4.94 32.30 0.105 -0.093 orangeR5.1lm6w25rh18dr01ms 5.1 yes 25 10 18 0.091 0.10 0.200 0.43 4.97 0.818 -0.672 greenR14.4lm6w25rh36dr005ms 14.4 yes 25 10 36 0.091 0.05 0.200 0.51 7.74 0.825 -0.365 cyanR14.4lm6w25rh36dr02ms 14.4 yes 25 10 36 0.091 0.20 0.200 0.47 12.91 0.845 -0.222 blueR14.4lm6w20rh36dr01ms anis 14.4 yes 20 10 36 0.091 0.10 0.200 1.06 18.88 0.485 -0.525 indigo anis FP/isot SPN1851lm6.3w5rh34dr01 NGC1851-like no 5 20 34 0.091 0.10 0.094 0.75 5.73 0.710 -0.111 brownN1851lm6.3w25rh17dr01ms NGC1851-like yes 25 20 17 0.091 0.10 0.047 0.79 4.06 0.824 -0.072 brown-dashed
Accretion onto proto-stellar disks scenarioR14.4lm6rh1 14.4 yes 25 10 1 0.012 0.57 0.005 4.17 4.75 0.010 0.001 red LX16R14.4lm6rh4 14.4 yes 25 10 4 0.002 0.43 0.022 4.25 13.42 0.001 -0.188 magenta LX16R14.4lm6rh1td20 14.4 yes 25 10 1 0.082 0.70 0.005 4.35 4.32 0.121 0.019 black 20R14.4lm5rh1td20 14.4 yes 25 1 1 0.052 0.73 0.005 0.35 5.52 0.075 0.114 blue 20R5.1lm6rh1td20 5.1 yes 25 10 1 0.082 0.70 0.005 4.23 4.15 0.121 0.014 green 20N1851lm6rh1td20 NGC1851-like yes 25 10 1 0.082 0.70 0.003 4.07 4.13 0.121 0.012 brown 20 M N R A S , ( ) Sollima
REFERENCES
Aguilar L., Hut P., Ostriker J. P., 1988, ApJ, 335, 720Allen C., Moreno E., Pichardo B., 2006, ApJ, 652, 1150Bastian N., Lamers H. J. G. L. M., de Mink S. E., Longmore S. N.,Goodwin S. P., Gieles M., 2013, MNRAS, 436, 2398Bastian N., Cabrera-Ziri I., Salaris M., 2015, MNRAS, 449, 3333Bastian N., Lardo C., 2018, ARA&A, 56, 83Baumgardt H., Makino J., 2003, MNRAS, 340, 227Baumgardt H., Hilker M., Sollima A., Bellini A., 2019, MNRAS,482, 5138Bellini A., Vesperini E., Piotto G., Milone A. P., Hong J., AndersonJ., van der Marel R. P., et al., 2015, ApJL, 810, L13Bellini A., Libralato M., Bedin L. R., Milone A. P., van der MarelR. P., Anderson J., Apai D., et al., 2018, ApJ, 853, 86Cabrera-Ziri I., Bastian N., Longmore S. N., Brogan C., HollyheadK., Larsen S. S., Whitmore B., et al., 2015, MNRAS, 448, 2224Calura F., D’Ercole A., Vesperini E., Vanzella E., Sollima A., 2019,MNRAS, 489, 3269Carretta E., Bragaglia A., Gratton R. G., Lucatello S., CatanzaroG., Leone F., Bellazzini M., et al., 2009, A&A, 505, 117Carretta E., Bragaglia A., Gratton R. G., Recio-Blanco A., Lu-catello S., D’Orazi V., Cassisi S., 2010, A&A, 516, A55Cassisi S., Salaris M., 2014, A&A, 563, A10Chatterjee S., Fregeau J. M., Umbreit S., Rasio F. A., 2010, ApJ,719, 915Cordero M. J., H´enault-Brunet V., Pilachowski C. A., Balbinot E.,Johnson C. I., Varri A. L., 2017, MNRAS, 465, 3515Cordoni G., Milone A. P., Mastrobuono-Battisti A., Marino A. F.,Lagioia E. P., Tailo M., Baumgardt H., et al., 2020, ApJ, 889,18Dalessandro E., Cadelano M., Vesperini E., Martocchia S., FerraroF. R., Lanzoni B., Bastian N., et al., 2019, ApJL, 884, L24D’Antona F., Ventura P., Decressin T., Vesperini E., D’Ercole A.,2014, MNRAS, 443, 3302Decressin T., Meynet G., Charbonnel C., Prantzos N., Ekstr¨om S.,2007, A&A, 464, 1029Decressin T., Baumgardt H., Kroupa P., 2008, A&A, 492, 101de Juan Ovelar M., Kruijssen J. M. D., Bressert E., Testi L., Bas-tian N., C´anovas H., 2012, A&A, 546, L1de Mink S. E., Pols O. R., Langer N., Izzard R. G., 2009, A&A,507, L1D’Ercole A., Vesperini E., D’Antona F., McMillan S. L. W., RecchiS., 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, 1521Duquennoy A., Mayor M., 1991, A&A, 500, 337Fregeau J. M., Rasio F. A., 2007, ApJ, 658, 1047Gieles M., Charbonnel C., Krause M. G. H., H´enault-Brunet V.,Agertz O., Lamers H. J. G. L. M., Bastian N., et al., 2018,MNRAS, 478, 2461Giersz M., 1998, MNRAS, 298, 1239Giersz M., 2001, MNRAS, 324, 218Giersz M., Heggie D. C., Hurley J. R., Hypki A., 2013, MNRAS,431, 2184Gratton R., Bragaglia A., Carretta E., D’Orazi V., Lucatello S.,Sollima A., 2019, A&ARv, 27, 8Gunn J. E., Griffin R. F., 1979, AJ, 84, 752Harris W. E., 1996, AJ, 112, 1487H´enault-Brunet V., Gieles M., Agertz O., Read J. I., 2015, MN-RAS, 450, 1164H´enon M. H., 1971, Ap&SS, 14, 151Hong 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, 4507 Hut P., Bahcall J. N., 1983, ApJ, 268, 319Johnston K. V., Spergel D. N., Hernquist L., 1995, ApJ, 451, 598Joshi K. J., Rasio F. A., Portegies Zwart S., 2000, ApJ, 540, 969Joshi K. J., Nave C. P., Rasio F. A., 2001, ApJ, 550, 691Khalaj P., Baumgardt H., 2015, MNRAS, 452, 924King I. R., 1966, AJ, 71, 64Krause M., Charbonnel C., Decressin T., Meynet G., Prantzos N.,2013, A&A, 552, A121Kremer K., Ye C. S., Rui N. Z., Weatherford N. C., Chatterjee S.,Fragione G., Rodriguez C. L., et al., 2020, ApJS, 247, 48Kroupa P., 2001, MNRAS, 322, 231Kruijssen J. M. D., 2009, A&A, 507, 1409Lada C. J., Muench A. A., Haisch K. E., Lada E. A., Alves J. F.,Tollestrup E. V., Willner S. P., 2000, AJ, 120, 3162Langer G. E., Hoffman R., Sneden C., 1993, PASP, 105, 301Lardo C., Bellazzini M., Pancino E., Carretta E., Bragaglia A.,Dalessandro E., 2011, A&A, 525, A114Lee H. M., Nelson L. A., 1988, ApJ, 334, 688Li M., Xiao L., 2016, ApJ, 820, 36Lind K., Charbonnel C., Decressin T., Primas F., Grundahl F.,Asplund M., 2011, A&A, 527, A148Lynden-Bell D., Wood R., 1968, MNRAS, 138, 495Merritt D., 1985, AJ, 90, 1027Milone A. P., Piotto G., Bedin L. R., Cassisi S., Anderson J.,Marino A. F., Pietrinferni A., et al., 2012, A&A, 537, A77Milone A. P., Piotto G., Renzini A., Marino A. F., Bedin L. R.,Vesperini E., D’Antona F., et al., 2017, MNRAS, 464, 3636Morscher M., Umbreit S., Farr W. M., Rasio F. A., 2013, ApJL,763, L15Nakamoto T., Nakagawa Y., 1994, ApJ, 421, 640Osipkov L. P., 1979, PAZh, 5, 77Ostriker J. P., Spitzer L., Chevalier R. A., 1972, ApJL, 176, L51Owen J. E., Clarke C. J., Ercolano B., 2012, MNRAS, 422, 1880Pancino E., Galfo A., Ferraro F. R., Bellazzini M., 2007, ApJL,661, L155Pietrinferni A., Cassisi S., Salaris M., Castelli F., 2006, ApJ, 642,797Piotto G., 2006, IAUJD, JD 14, id. 10Piotto G., Milone A. P., Bedin L. R., Anderson J., King I. R.,Marino A. F., Nardiello D., et al., 2015, AJ, 149, 91Planck Collaboration, Ade P. A. R., Aghanim N., Arnaud M., Ash-down M., Aumont J., Baccigalupi C., et al., 2016, A&A, 594,A13Renzini A., 2008, MNRAS, 391, 354Sollima A., Bellazzini M., Lee J.-W., 2012, ApJ, 755, 156Sollima A., Mastrobuono Battisti A., 2014, MNRAS, 443, 3513Sollima A., Ferraro F. R., 2019, MNRAS, 483, 1523Spitzer L., Shull J. M., 1975, ApJ, 201, 773Spitzer L., 1987, in ”Dynamical Evolution of globular clusters”,Princeton University Press, Princeton NJ”Tiongco M. A., Vesperini E., Varri A. L., 2019, MNRAS, 487, 5535Tognelli E., Prada Moroni P. G., Degl’Innocenti S., 2011, A&A,533, A109Vasiliev E., 2015, MNRAS, 446, 3150Vasiliev E., 2019, MNRAS, 484, 2832Ventura P., Di Criscienzo M., Carini R., D’Antona F., 2013, MN-RAS, 431, 3642Vesperini E., 1997, MNRAS, 287, 915Vesperini 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., Hong J., Webb J. J., D’Antona F., D’Ercole A., 2018,MNRAS, 476, 2731Wijnen T. P. G., Pols O. R., Pelupessy F. I., Portegies Zwart S.,2016, A&A, 594, A30Yoshida H., 1990, PhLA, 150, 262MNRAS000