# Testing lowered isothermal models with direct N-body simulations of globular clusters - II: Multimass models

Miklos Peuten, Alice Zocchi, Mark Gieles, Vincent Hénault-Brunet

MMNRAS , 1–28 (2017) Preprint 21 May 2018 Compiled using MNRAS L A TEX style ﬁle v3.0

Testing lowered isothermal models with direct N -bodysimulations of globular clusters - II. Multimass models M. Peuten, (cid:63) A. Zocchi, , M.Gieles, V. H´enault-Brunet Department of Physics, University of Surrey, Guildford, GU2 7XH, UK Dipartimento di Fisica e Astronomia, Universit`a degli Studi di Bologna, viale Berti Pichat 6/2, I40127, Bologna, Italy Department of Astrophysics/IMAPP, Radboud University, PO Box 9010, NL-6500 GL Nijmegen, the Netherlands

ABSTRACT

Lowered isothermal models, such as the multimass Michie-King models, have beensuccessful in describing observational data of globular clusters. In this study we assesswhether such models are able to describe the phase space properties of evolutionary N -body models. We compare the multimass models as implemented in limepy (Gieles& Zocchi) to N -body models of star clusters with diﬀerent retention fractions for theblack holes and neutron stars evolving in a tidal ﬁeld. We ﬁnd that multimass modelssuccessfully reproduce the density and velocity dispersion proﬁles of the diﬀerent masscomponents in all evolutionary phases and for diﬀerent remnants retention. We furtheruse these results to study the evolution of global model parameters. We ﬁnd that overthe lifetime of clusters, radial anisotropy gradually evolves from the low-mass to thehigh-mass components and we identify features in the properties of observable starsthat are indicative of the presence of stellar-mass black holes. We ﬁnd that the modelvelocity scale depends on mass as m − δ , with δ (cid:39) . m can be shallower, depending on thedark remnant content, and agrees well with that of the N -body models. The reportedmodel parameters, and correlations amongst them, can be used as theoretical priorswhen ﬁtting these types of mass models to observational data. Key words: methods: numerical – stars: black holes –stars: kinematics and dynamics– globular clusters: general –galaxies: star clusters: general.

The amount of available data from observations of globu-lar clusters (GCs) is steadily increasing. With the arrivalof the ESA–Gaia data (Gaia Collaboration et al. 2016), weare entering the era of high-precision kinematics, allowingus to study properties of GCs with unprecedented detail.This calls for adequate methods of analysing and describingthem in an equally detailed way. Despite the fact that GCsare thought to be free of dark matter (Baumgardt et al. 2010;Ibata et al. 2013; Baumgardt 2017), and to have evolved tospherical and isotropic conﬁgurations as the result of two-body relaxation, GCs are complex systems to model. Theyconsist of stars and stellar remnants with diﬀerent massesand luminosities and primordial and dynamically processedbinary stars (Heggie 1975; Goodman & Hut 1989; Hut et al.1992; Heggie et al. 2006; Trenti et al. 2007). The mass andluminosity functions depend on the stellar initial mass func-tion (IMF), age and metallicity. GC stellar populations dis-play chemical anomalies (Gratton et al. 2004) and broadened (cid:63)

E-mail: [email protected] (MP) main sequences (MS), possibly the result of variations in thehelium abundance (Milone et al. 2014). Furthermore GCsevolve in a galactic tidal ﬁeld that inﬂuences their evolu-tion and present-day properties (Chernoﬀ & Weinberg 1990;Johnston et al. 1999; Takahashi & Portegies Zwart 2000;Baumgardt & Makino 2003; K¨upper et al. 2010; Rieder et al.2013; Renaud et al. 2017).Modelling GCs on a star-by-star basis using direct N -body models has only become possible recently: Hurley et al.(2005) presented the ﬁrst N -body simulation of an open clus-ter, Zonoozi et al. (2011) modelled a low-mass GC and ﬁnallyHeggie (2014) and Wang et al. (2016) presented the ﬁrst N -body simulations of GCs with N ∼ . The faster MonteCarlo method allows to explore the parameter of the initialconditions to some extent (Heggie & Giersz 2008; Gierszet al. 2013). To infer properties for a large number of GCswith models with several degrees of freedom, static modelsthat are fast to calculate are required. By using relativelysimple models, that are motivated by the underlying physi-cal processes that drive their evolution, diﬀerences betweenmodels and observations can be used to increase our under-standing (Binney & McMillan 2011). c (cid:13) a r X i v : . [ a s t r o - ph . GA ] J u l M. Peuten

In the context of GCs, the King (1966) models are oftencompared to observations, although they cannot describe allGCs successfully. For example, McLaughlin & van der Marel(2005) ﬁnd that the more extended Wilson (1975) models arebetter in describing the surface brightness proﬁles of someGalactic GCs. In addition, both King and Wilson modelshave isothermal cores, which are not able to describe thelate stages of core collapse (Lynden-Bell & Eggleton 1980;Cohn 1980). The models we are going to test and discussin the context of GCs are multimass, anisotropic and spher-ical models (hereafter multimass models), which describethe properties of GCs considering their stellar mass func-tion (MF) in the form of mass bins. This formulation allowsfor diﬀerent behaviour of the diﬀerent components. Thesemodels are deﬁned by a distribution function (DF) whichis a solution of the collisionless Boltzmann equation assum-ing a Maxwellian velocity distribution that is ‘lowered’ tomimic the eﬀect of a negative escape energy as the resultof the galactic tides. The multimass formulation of a Kingmodel was ﬁrst introduced by Da Costa & Freeman (1976,we note that a formulation of a multimass model was al-ready presented in Oort & van Herk 1959). Gunn & Griﬃn(1979) extended these models including radial anisotropy asformulated by Eddington (1915) for isothermal models (seealso Michie 1963).Since their introduction, multimass models have beensuccessfully used in a multitude of studies such as in Illing-worth & King (1977), Pryor et al. (1986), Lupton et al.(1987), Meylan (1987), Richer & Fahlman (1989), Meylan& Mayor (1991), Meylan et al. (1995), Sosin (1997), Piotto& Zoccali (1999) and Richer et al. (2004) to name a few.More recently, they have been used in observational studies,such as those by Paust et al. (2010), Beccari et al. (2010),Sollima et al. (2012), Beccari et al. (2015) and Sollima et al.(2017), as well as in theoretical studies (Takahashi & Lee2000; Sollima et al. 2015).Davoust (1977) realized that the DFs of the Woolley(1954), King (1966) and Wilson (1975) models can be writ-ten as a single DF with one additional integer parame-ter. Gomez-Leyton & Velazquez (2014) further generalizedthis formulation, allowing to calculate models in betweenthe three classical models. (Gieles & Zocchi 2015, here-after GZ15) took up these formulation and added radialanisotropy as deﬁned in the Michie–King models (Michie1963) and multiple mass components as in Gunn & Griﬃn(1979). GZ15 introduced a power-law dependence betweenmass and anisotropy radius for each mass bin, while Gunn &Griﬃn (1979) argued that was not necessary because mostevents that inﬂuence the anisotropy are mass independent ornot very important. GZ15 also implemented the possibilityto change the degree of mass segregation with an additionalparameter δ that describes the relation between the velocityscale and mass, which in most models is assumed to be equalto 1 / N MB +5 parameters and2 scales, with N MB being the number of mass bins, comparedto 3 parameters and 2 scales for the single-mass model. It is therefore easier to ﬁt multimass models to the data becausethey have more degrees of freedom (McLaughlin 2003). Notonly the selection of the right number of mass bins, butalso how they are deﬁned is criticized as a ‘usual compro-mise between convenience and realism’, as Meylan & Heggie(1997) put it. Given the numerous studies successfully us-ing multimass models, this problem does not seem to be toomuch of a concern, but we nevertheless explored it in ourstudy. Another assumption for which multimass models havebeen criticized is the assumption of equipartition of energy(McLaughlin 2003; Trenti & van der Marel 2013). Indeedthe velocity scale is usually assumed to scale with the massas m − / , but we note that evolutionary multimass mod-els only achieve partial equipartition (Merritt 1981; Miocchi2006; GZ15; Bianchini et al. 2016) as the result of the escapevelocity.Several aspects of multimass models were already anal-ysed with the help of Fokker–Planck (Takahashi & Lee 2000)and N -body simulations (Sollima et al. 2015). The goal ofthis study is to compare the multimass models in the for-mulation by GZ15 to a set of N -body models to assess thequality of the former and to analyse whether some of theabove mentioned criticism is justiﬁed. In this comparisonwe do not include any source of uncertainties, such as obser-vational biases, to see how good the models are under idealconditions. Hence, we determine the MF of the multimassmodels directly from the N -body data.The comparison is done by ﬁtting the multimass mod-els to snapshots from diﬀerent N -body models using aMarkov Chain Monte Carlo (MCMC) method. Additionally,we study the new parameters which are now available in theextended formulation of the models by GZ15. In particu-lar, the continuous truncation parameter as introduced byGomez-Leyton & Velazquez (2014) and the parameter thatcontrols the mass dependence of the anisotropy for each massbin. Furthermore, we study the behaviour of the mass seg-regation parameter δ , which in previous studies was ﬁxed to δ = 1 /

2. By letting this parameter free, we can test whetherthis assumption is justiﬁed. By varying the amount of stellar-mass black holes (BHs) and neutron stars (NSs) retained inthe diﬀerent N -body models, we also study their impact onthe cluster as well as on the diﬀerent parameters of the best-ﬁtting models.In Zocchi et al. (2016), a similar analysis was presentedfor single-mass models. This comparison showed that thesingle-mass models are successful in describing the diﬀerentphases of the dynamical evolution. Zocchi et al. (2016) stud-ied the development of radial anisotropy in GCs and foundthat the models can be used to put limits on the expectedamount of radial anisotropy.This paper is structured as follows: in the next section,we give a brief overview of the multimass models. Then inSection 3, we discuss how the N -body models were gener-ated and we discuss their properties. In Section 4, we presentthe method used for the analysis and the challenges we en-countered. The radial proﬁles of density, velocity dispersionand anisotropy are discussed in Section 5. In Section 6, wediscuss the values of the best-ﬁtting model parameters andscales, and their implications. Finally, in Section 7, we dis-cuss our results and present our conclusions. MNRAS , 1–28 (2017) esting isothermal models - II. Multimass The multimass models used in this study are provided by the limepy (Lowered Isothermal Model Explorer in

PYthon ) software package (GZ15). A model has diﬀerent components,each representing stars in a mass range, characterized by amean and total mass. The DF of the j th mass component isgiven by f j (cid:0) E, J (cid:1) = A j exp (cid:18) − J r ,j s j (cid:19) E γ (cid:18) g, − E − φ ( r t ) s j (cid:19) , (1)for E < φ ( r t ) and 0 otherwise. The speciﬁc energy E = v / φ ( r ) is one of the two integrals of motion, where v isthe velocity and φ ( r ) the speciﬁc potential at a distance r from the centre. The energy E is lowered by the potential atthe truncation radius φ ( r t ). The function E γ is deﬁned as E γ ( a, x ) = (cid:40) exp ( x ) a = 0exp ( x ) P ( a, x ) a > P ( a, x ) ≡ γ ( a, x ) / Γ ( a ) the regularized lower incom-plete gamma function. The other integral of motion is thespeciﬁc angular momentum J = rv t , where v t is the tangen-tial component of the velocity vector.The anisotropy radius r a is a parameter that controlshow anisotropic the model is. The system is isotropic in thecentre, radially anisotropic in the intermediate part and near r t it is isotropic again. For small values of r a , the models arestrongly anisotropic and for values of r a lager than r t , themodels are completely isotropic. GZ15 include a power-lawdependence between mass and anisotropy radius: r a ,j = r a µ ηj (3)with µ j the dimensionless mean mass of stars in the j th masscomponent, deﬁned as: µ j = m j ¯ m (4)where m j is the mean mass of stars in the j th component and¯ m a reference mass which we set equal to the global meanmass. If η is set to zero the anisotropy radius is independentof the mass as in Gunn & Griﬃn (1979). limepy expects r a to be input in units of the King radius ( r ), ˆ r a = r a /r ,hence ˆ r a is the parameter we vary.The truncation parameter g was introduced by Gomez-Leyton & Velazquez (2014) and describes the polytropic partnear the escape energy. The polytropic index n relates to g as n = g + 3 / g = 0 and r a (cid:29) r t , theDF is identical to the one from the Woolley model (Woolley1954). A Michie–King (Michie 1963; King 1966) model isreproduced for g = 1 and for g = 2 one gets the non-rotatingWilson model (Wilson 1975). The range of possible values forthe model parameter are 0 ≤ g ≤ . g = 3 .

5. The ﬁnal parameter neededto deﬁne the models is the dimensionless central potential W (King 1966, ˆ φ in GZ15) which speciﬁes how centrallyconcentrated the model is. It is a boundary condition forsolving Poisson’s equation. https://github.com/mgieles/limepy Besides these parameters, there are also two constantswhich deﬁne the physical scales of the model: one is theglobal velocity scale s and the other is the normalizationconstant A which sets the phase space density. Instead ofthese scales, the code needs as input the total cluster mass M Cl and a radial scale r scale (which can be r , the half-massradius r h , the viral radius r v or r t ), which are internallyconverted to A and s .The velocity scale s j is deduced from s as: s j = sµ − δj (5)It is usually assumed that δ = 1 /

2, but in this study wedetermine the value of this parameter from the ﬁts to the N -body models.The constants s j and A j are connected to the mass ineach component ( M j ), which the user provides together with m j . It must be noted that M Cl is a required input parameter,independent from the M j parameters, because the latter areonly used to compute the relative masses in each component.Only after the model is solved, (cid:80) j M j = M Cl .Given these ﬁve parameters ( g , W , δ , r a and η ) andtwo scales ( M Cl , r scale ) together with the description of themass bins ( M j , m j ) limepy ﬁrst calculates the density foreach mass bin via: ρ j = (cid:90) f j (cid:0) E, J (cid:1) d v (6)Then, the dimensionless Poisson equation is solved ∇ ˆ φ = − (cid:88) j α j ˆ ρ j (7)with α j = ρ j, /ρ , ˆ ρ j = ρ j /ρ j, and the dimensionless pos-itive potential ˆ φ = ( φ ( r t ) − φ ) /s , is iteratively solved byvarying α j until the calculated M j converges to the inputvalues. After the model is solved, it is scaled to M Cl and r scale . We can then ﬁnd the likelihood for any phase spacecoordinate using the DF (equation 1).In equation (4), we set ¯ m to the mean mass of thecluster. In the formulation by Da Costa & Freeman (1976);Gunn & Griﬃn (1979) and GZ15, ¯ m is the central densityweighted mean mass. After performing several comparisons,we found that models calculated by the two diﬀerent formu-lations give the same results within the numerical uncertain-ties. Furthermore we found that using the global mean massinstead speeds up the calculation, especially for models withBHs. When using the global ¯ m the meaning of two modelparameters is modiﬁed compared to Da Costa & Freeman(1976); Gunn & Griﬃn (1979) and GZ15: W and r a bothrepresent their value for a hypothetical mass group with amass of ¯ m . Besides computational improvement, this changefrom the original formulation also allows us to compare themultimass W value with the single-mass W value, as bothrepresent the W value for the mean mass group.One can easily translate the values given in one ¯ m deﬁni-tion ( W , ˆ r a ) to another ¯ m ∗ deﬁnition ( W ∗ , ˆ r ∗ a ) by applyingthe following two equations: W ∗ = W (cid:18) ¯ m ∗ ¯ m (cid:19) δ (8)ˆ r ∗ a = ˆ r a (cid:18) ¯ m ∗ ¯ m (cid:19) ( η + δ ) (9) MNRAS , 1–28 (2017)

M. Peuten

The δ term in equation (9) comes from the r dependenceof ˆ r a .As further improvement to the original formulation ofthe limepy models we found that radially anisotropic mod-els can be constructed faster if one ﬁrst calculates the M j array of the corresponding isotropic model and then usesthis model as starting point to solve the anisotropic model.As with the previous improvement the diﬀerences are onlyof numerical nature. This procedure is now implemented inthe current distribution of limepy . N -BODY MODELS For the computation of the N -body data, we use the ap-proach presented in Trenti et al. (2010): the stellar evolutionis done ﬁrst and separately from the dynamical evolution.We do this because Galactic GCs have diﬀerent dynamicalages but have all roughly the same physical age of around12 Gyr. We consider models with diﬀerent retention fractionsof NSs and BHs and analyse them at various dynamical ages.Temporal units are always expressed in units of the initialhalf-mass relaxation time ( τ rh , ) of the N -body model. N -body models For this analysis, we run four N -body models, with diﬀer-ent amounts of NSs and BHs. Each N -body model was setup as a cluster with N = 10 stars initially following theH´enon isochrone model (H´enon 1959) with r h = 2 .

25 pc. AsIMF we adopted a Kroupa IMF (Kroupa 2001) in the massrange between 0 . (cid:12) and 100 M (cid:12) without any primordialbinaries. Then, by using the ﬁtting formula by Hurley et al.(2000) and assuming a metallicity of [Fe / H] = −

2, the starswere evolved to an age of about 12 Gyr.We mimic the eﬀect of supernova kick velocity by re-moving a certain fraction of NSs and BHs from the initialconditions described above. The retention fraction of NSsand BHs after supernova kicks is highly uncertain (Repettoet al. 2012; Mandel 2016). To bracket all possible cases, weconsider four diﬀerent values for the fraction of remnantsthat we retain in the cluster: 100% (all the remnants areretained, simulation N1), 33% (simulation N0.3), 10% (sim-ulation N0.1) and 0% (all the remnants are removed, simu-lation N0). The initial half-mass relaxation time for all fourclusters was τ rh , = 350 Myr before the stellar evolution andthe removal of the dark remnants, after these steps the τ rh , values are 412 Myr for N1, 426 Myr for N0.3, 427 Myr forN0.1 and 428 Myr for N0.The clusters are evolved on a circular orbit with a circu-lar velocity of V c = 220 km / s at a distance of R G = 4 kpc, ina singular isothermal galactic potential to mimic a galaxy.The equation of motion is solved in an inertial referenceframe centred on the cluster.These four stellar systems were then dynamicallyevolved with the state-of-the art N -body integrator nbody6 (Aarseth 2003), in the variant with GPU support (Nitadori& Aarseth 2012), until total dissolution of each cluster, i.e.until less than 100 objects are left in the cluster. Every ob-ject reaching a distance greater than twice the Jacobi ra-dius ( r J ) is considered lost and is removed from the N -body model. As the stellar evolution is done before the actual N -body simulation, binaries which formed in the course of thesimulations were also only evolved dynamically.A snapshot of each cluster is taken every Gyr, resultingin 48 snapshots (11 for model N1, 13 for model N0.3, 12 formodel N0.1 and 12 for model N0) which we ﬁt the multimassmodels to (Section 4). Because multimass models describe bound objects in a clus-ter, we removed any unbound object from the N -body mod-els. We discuss here how we selected the unbound objectsfor each N -body snapshot.First, we determine the Jacobi radius r J = (cid:18) GM Cl (cid:19) / (10)in which M Cl is the total mass within r J and Ω = V c /R G is the orbital angular velocity. As a ﬁrst guess, we set M Cl equal to the total mass of all stars in the snapshot and thendetermine r J through an iterative approach.With r J determined we are now able to calculate thespeciﬁc critical energy which is equal to the potential at r J E crit = φ ( r J ) = − GM Cl r J (11)The true critical energy is diﬀerent as equation (11) neglectsthe tides. We adopted this deﬁnition nevertheless to be con-sistent with the multimass models, which also do not accountfor the changed potential due to tidal eﬀects. We consideredan object bound if it is within r J and for its energy it holds E i < E crit , and we only used these bound objects in the restof this analysis. N -body models Fig. 1 shows how M Cl for the four diﬀerent N -body modelsevolves over the course of the simulation. It is apparent thatthe cluster with 100% initial BH and NS retention (simula-tion N1) has the highest initial mass loss, but there seems tobe no direct correlation between the number of BHs and NSsand initial mass loss as can be seen from the other three N -body models. Over the course of evolution the four modelsseem to have aligned their mass-loss rate which is in accor-dance with the ﬁndings of Lee & Ostriker (1987) and Gieleset al. (2011) that the escape rate of clusters with the samemass mainly depends on the tidal ﬁeld, which is the samefor all four models.In Fig. 2, we have plotted the evolution of r h for thefour diﬀerent N -body models. As can be seen in the ﬁgure,increasing the retention fraction of the BHs leads to an ex-pansion of the cluster: simulation N1, with 100% NS andBH retention has an r h which is on average twice as largeas the r h from simulation N0 with no NS and BH retention.The cluster in simulation N0.3 loses all of its BHs at around7 τ rh , and for the rest of the simulation its r h resembles that The snapshots can be retrieved from http://astrowiki.ph.surrey.ac.uk/dokuwiki/doku.php?id=tests:collision:mock_data:challenge_2

MNRAS , 1–28 (2017) esting isothermal models - II. Multimass M C l [ M ⊙ ] Age [ τ rh,0 ] N1N0.3N0.1N0

Figure 1.

Evolution of the cluster mass M Cl for the four N -bodymodels. Age is given in units of their τ rh , . r h [ p c ] Age [ τ rh,0 ] N1N0.3N0.1N0

Figure 2.

Evolution of the cluster half-mass radius r h for thefour N -body models. Age is given in units of their τ rh , . of N0. The eﬀect of stellar-mass BHs on the radius evolutionwas also described by Mackey et al. (2008). The global evo-lution of r h however is essentially the same, independent ofBHs and NSs retained, and follows the description in Gieleset al. (2011). In the ﬁrst half of their lifetime, the clusters arein the expansion-dominated phase while in the second halfthe clusters are in the evaporation-dominated phase duringwhich r h decreases again until total dissolution.Fig. 3 shows the relative number evolution for the BHsand NSs in the simulation N1, N0.3 and N0.1. As can be seenin the bottom ﬁgure, the cluster from N -body simulationN1 with 100% initial BHs and NSs retention is the onlycluster that retains its BHs almost until to the end of itslifetime, while the cluster from N -body simulation N0.3 losesall its BHs at around 7 τ rh , and the one from simulation N0.1already at 2 τ rh , . Looking at the NSs in the top of Fig. 3,we see their initial loss is not as strong as for the BHs. Butas soon as all BHs have left the cluster, the NSs escaperate increases such that the cluster N0.1 loses all its NS ataround ∼ τ rh , . Only the cluster from N -body simulationN1 has a population of NS left at the end of its lifetime. R e l a t i v e nu m be r o f N S [ % ] N1N0.3N0.1 R e l a t i v e nu m be r o f B H [ % ] Age [ τ rh,0 ] N1N0.3N0.1

Figure 3.

Evolution of the relative number of NS (top) and BH(bottom) as function of time in units of τ rh , for the three modelswhich initially retained BHs and NSs. After all BHs are lost, NSs, then being the most massiveobjects in the cluster, segregate to the centre and are thenejected from the cluster due to interactions they experiencewith each other.Tables A1–A4 list various properties of the diﬀerentsnapshots, such as the dynamical age, the bound mass andthe number of NSs and BHs.

As in Peuten et al. (2016), we ﬁnd that the mean massproﬁle is independent of the remnant retention fraction. InFig. 4, we plot this for all four N -body models at four diﬀer-ent times in their evolution: 2 . τ rh , , 7 . τ rh , , 16 . τ rh , and26 τ rh , . Looking at the diﬀerent times we see that the over-all behaviour is the same for all models independent of theirdark remnant population. Some divergences between the dif-ferent N -body models can be seen in the ﬁrst snapshot at2 . τ rh , but over the course of evolution these diﬀerencesdiminish. The proﬁles get ﬂatter over time. This is compa-rable to the behaviour found for a set of N -body modelswhere the dynamical and stellar evolution were done con-currently in Peuten et al. (2016). Here, the evolution overtime is less strong because the stellar evolution was done be-fore the dynamical evolution. We are not aware of a theoryproviding an explanation for this attractor solution of ¯ m ( r ),but for single-mass system it is known that after severalrelaxation times the evolution becomes self-similar (H´enon1961, 1965). Also it had been shown that the evolution of MNRAS , 1–28 (2017)

M. Peuten < m (r) > / < m > r / r h τ rh,0 N1N0.3N0.1N0 < m (r) > / < m > r / r h τ rh,0 < m (r) > / < m > r / r h τ rh,0 < m (r) > / < m > r / r h τ rh,0 Figure 4.

Relative mean mass (i.e. mean mass of stars in radialbins divided by the total mean mass) as a function of the distancefrom the cluster centre in units of r h , for all simulations at fourdiﬀerent times: 2 . τ rh , , 7 . τ rh , , 16 . τ rh , and 26 τ rh , . Triangles(red), circles (cyan), boxes (blue) and stars (black) refer to simu-lation N0, N0.1, N0.3 and N1, respectively. Error bars denote 1 σ uncertainties. radii and mass of the multimass systems is comparable tothose of single-mass systems (Lee & Goodman 1995; Gieleset al. 2010) but faster. Furthermore Giersz & Heggie (1996,1997) showed in multimass N -body models that after sometime the mean mass in Lagrangian shells stops evolving and,to a ﬁrst approximation, stays constant. Although we do notexplore this here, this result could be used as a theoreticalprior when comparing multimass models to data. To determine the best-ﬁtting multimass models for eachsnapshot we use the MCMC software package emcee (Foreman-Mackey et al. 2013), which is a pure- python im-plementation of the Goodman & Weare’s Aﬃne InvariantMCMC Ensemble sampler (Goodman & Weare 2010). The python implementation makes it straightforward to coupleit with limepy . Furthermore, we beneﬁt from the fact thatwe created a distributed grid computing version of python ’s map() function, thereby dynamically distributing eﬃcientlythe workload from several MCMC runs over all availableCPU cores at the University of Surrey Astrophysics com-puting facilities.The ﬁtting process consists of computing a multimassmodel based on the input parameters provided by theMCMC walker position in parameter space and the currentmass bin description (see the next section). Then the likeli-hood for each star in all mass bins is calculated using the DF(see equation 14) and the phase space position of the starfrom the N -body snapshot. By randomly varying the walkerpositions in parameter space, the MCMC algorithm tries toﬁnd those parameters which maximize the product of theseindividual likelihoods. The best-ﬁtting value of each param-eter is estimated as the median of the marginalized posteriordistribution using all walker positions from all chains afterremoving the initial burn-in phase. This generally coincideswith the value of the parameter providing the largest like-lihood. For the 1 σ errors, we use the values from the 16thand 84th percentiles. As mentioned in Section 1, the mass bin selection for multi-mass models is in most cases a choice of convenience (Meylan& Heggie 1997) as throughout the literature there is no gen-eral rule on how to select the best. This can be partiallyexplained by the fact that most publications consider dif-ferent data for their analyses, and have diﬀerent researchtargets, leading to diﬀerent approaches on how to set up theMF. However, we do know everything about our N -bodymodels, and this allows us to test the mass bin selectionfor multimass models. In particular, we want to understandwhat the minimum number of mass bins is to get a stableresult and how to choose the bins.For this analysis we use the N -body snapshot of simu-lation N1 at the time of 2 . τ rh , . As we wanted to trace theoverall evolution of the diﬀerent star types, we opted againstmixing them and therefore we give every star type at leastone bin. This means that we have at least ﬁve mass bins, MNRAS , 1–28 (2017) esting isothermal models - II. Multimass one each for the MS stars, the evolved stars (ES), the whitedwarfs (WD), NSs and BHs. Looking at the BHs, NSs andESs we decided to not further split them into several binsgiven the fact that they have either a rather small range ofpossible masses and/or are to low in number to justify thesplit. This leaves us with MS stars and WDs which are bothnumerous and do have a large mass range: 0 . − .

83 M (cid:12) forthe MS stars and 0 . − .

44 M (cid:12) for the WDs in our N -bodysnapshots.First, we determine how the bin selection inﬂuences theresults of the analysis. For this, we choose four diﬀerent bin-ning methods: a logarithmic binning, a linear binning, a bin-ning where in each mass bin there is an equal number of starsand a binning where there is the same amount of mass ineach mass bin. We ﬁxed the number of WD mass bins toone and then for each bin type we calculated the multimassmodel repeatedly with increasing number of MS star bins.The general idea here is that with increasing number of massbins, the overall parameters like M Cl , r h , etc., should con-verge to the value one would get for the ideal case, whereeach star has its own mass bin. We ﬁnd that the results arealmost independent of the way one chooses the binning: thenumber of MS mass bins needed to converge is the sameand the diﬀerence between the diﬀerent models are for allproperties generally less than 5%. We therefore choose forthe further analysis the logarithmic binning.Then we determine the minimum number of bins neededto get stable results, as increasing the number of bins alsoincreases the computation time of the models. We variedthe number of mass bins of the MS stars and the WDs in-dependently from each other. Here again, we see that withincreasing number of mass bins the diﬀerent quantities con-verge. We ﬁnd that we need at least two WDs mass bins andat least four MS stars mass bins for the diﬀerent quantitiesto converge. Increasing the number of bins any further doesnot improve the results (values are comparable within 5%).For our further analysis, we opt to use ﬁve MS stars massbins and three WDs mass bins. Therefore, in total we con-sider eleven mass bins (MS: 5; ES: 1; WD: 3; NS: 1; BH: 1)to set up our multimass models. Tables A5–A8 list the massbins for all N -body snapshots used. Before we present the results, we discuss a particularitywhich we encountered in our analysis. The potential of limepy models is spherical, however, the true potential ofthe cluster is triaxial because of the eﬀect of tides. Also theLagrange points of the cluster, through which stars can es-cape (Fukushige & Heggie 2000; Baumgardt 2001; K¨upperet al. 2010; Claydon et al. 2017), are not accounted for inthe multimass model. Therefore, the models are not able todescribe the objects near the critical energy correctly. Fromthis, it follows that some of the objects which are unboundin the true potential are still found bound in our deﬁnitionand cannot be described by the model correctly. These ob-jects pose a problem because they drive the ﬁt to unrealisticparameter values. In this work, every post MS star which is not a remnant iscalled an ES.

To cope with this problem, we introduced an artiﬁcialbackground population with a constant likelihood (i.e. uni-form distribution) in phase space. We added the artiﬁcialbackground population with a total mass of around 1% ofthe original cluster mass to the N -body snapshot. This back-ground population has the same MF as the cluster. The up-per limit for the maximal distance and velocity are chosento be twice the maximal values from the original snapshot( r max and v max ).We describe the likelihood function of the backgroundmodel as L B = M Back

V M tot (12)where M tot is the total mass of the snapshot including theartiﬁcial background and M Back is the mass of the back-ground only. The phase space volume V is deﬁned as V = 43 π (2 r max ) × π (2 v max ) (13)The total likelihood of an object for a given model iscalculated as L = f ( E, J ) M tot + M Back

V M tot (14)When integrated over the whole phase space volumewithin 2 v max and 2 r max the ﬁrst term equals to M Cl / M tot and the second to M Back / M tot , giving a total likelihood ofunity, as required.

We initiate the MCMC walkers in a randomly chosen spherein parameter space. For some snapshots, we run several ﬁtswith diﬀerent initial conditions to test for any divergence.We chose ﬂat priors restricted mainly by currently observedvalues for the parameters and/or by the range in which theyare considered physically valid. For the MCMC ﬁtting, westarted out with around 500 walkers and found good ﬁts forthe N -body model N0 without BHs and NSs. For the othermodels, prominently those with BHs, converging ﬁts wereonly achieved with at least 2000 walkers. On average, eachMCMC chain was run for 1000 iterations and convergencewas reached after around 300 iterations, which we trimmedfrom the MCMC chains for the calculation of the best-ﬁttingparameters. The MCMC chain took on average longer toconverge in snapshots with BHs than in snapshots without.In some cases, we also had to adjust the emcee scale pa-rameter a , which is generally set to 2, to increase the ac-ceptance rates (for details on how this aﬀects the MCMCalgorithm see the discussion in Foreman-Mackey et al. 2013,their Eq. 2).In Figs. 5 and 6, we show the marginalized posteriorprobability distribution for each parameter as well as the2D projections of the posterior probability distribution rep-resenting the covariance between the diﬀerent ﬁtting param-eters for two MCMC runs. Figs. 5 and 6 show the results ofthe ﬁtting to the N -body model N0 at 2 . τ rh , , and to the N -body model N1 at 17 . τ rh , , respectively. The obvious dif-ference between the two models is that for model N1 thereare two parameters, namely r a and η that do not converge toa single value. The stellar system in this particular N -bodysnapshot is isotropic: values of r a larger than r t generate MNRAS , 1–28 (2017)

M. Peuten

Figure 5.

Marginalized posterior probability distribution and 2D projections of the posterior probability distribution for the modelparameters and scales. This ﬁgure shows the results of the MCMC ﬁtting to the N -body model N0 at 2 . τ rh , . The dashed lines in themarginalized posterior probability distribution indicate the 16th, 50th and 84th percentiles. isotropic models, equally likely to reproduce the data, andfor this reason, the values of r a , and consequently of η , can-not be constrained.Looking at the 2D projections of the posterior proba-bility distribution in Figs. 5 and 6, one can see that they arenearly circular for most of the parameter pairs. This showsthat when using the full phase space information of eachstar, degeneracies between the diﬀerent parameters can bealleviated. Tables A9–A12 list the best-ﬁtting parameters for allthe N -body snapshots we considered. N -BODY MODELS In the ﬁrst part of our analysis, we compare the best-ﬁttingmultimass models with the results directly computed fromthe N -body snapshots for each model at all times. MNRAS , 1–28 (2017) esting isothermal models - II. Multimass Figure 6.

Marginalized posterior probability distribution and 2D projections of the posterior probability distribution for the modelparameters and scales. This ﬁgure shows the results of the MCMC ﬁtting to the N -body model N1 at 17 . τ rh , . The dashed lines inthe marginalized posterior probability distribution indicate the 16th, 50th and 84th percentiles. The best-ﬁtting values of η and r a areunconstrained because this stellar system is not radially anisotropic. First, we compare the mass density proﬁles of the best-ﬁtting multimass models and the N -body models. For this,we binned each mass bin of the N -body data such that ineach radial bin there are at least 30 objects and each ra-dial bin has a minimal radial width of 0 .

15 pc. We assumedPoisson errors for the uncertainties of the binned data anddeﬁne the position uncertainty by the 16th and 84th per- centiles of the distribution of the positions of the objects ineach bin. In Fig. 7, we compare the mass density proﬁles forthe three models N1, N0.3 and N0 at four diﬀerent times:2 . τ rh , which is the ﬁrst snapshot for each N -body model,7 . τ rh , , 16 . τ rh , and 26 τ rh , the snapshot at the end of theclusters lifetime. For clarity we only show one mass bin perstellar type. We did not include the results of model N0.1,since they are similar to the results of model N0. Togetherwith the best-ﬁtting result from the multimass models, we MNRAS , 1–28 (2017) M. Peuten

Figure 7.

Comparison of the mass density proﬁles for models N1, N0.3 and N0 at four diﬀerent ages: 2 . τ rh , , 7 . τ rh , , 16 . τ rh , and26 τ rh , . The points represent the binned N -body data, the thick lines represent the best-ﬁtting multimass models proﬁles and the thinlines represent the results from the walker positions at the last iteration. Red represents the MS stars, cyan – ESs, blue – WDs, pink –NSs and black – BHs. Error bars denote 1 σ uncertainties. MNRAS , 1–28 (2017) esting isothermal models - II. Multimass also plot the results from the walker positions at the lastiteration of the MCMC routine, reﬂecting the uncertaintiesof the results.As can be seen from Fig. 7, the best-ﬁtting multimassmodels reproduce the mass density proﬁles of the diﬀerentmass components. Diﬀerences are only found in the outer-most regions and innermost regions as well as for cases wherethe number of objects in a mass bin is low. For the outerregions, a diﬀerence is expected: as already discussed in Sec-tion 4.2, the models assume that the cluster is sphericallysymmetric, which is not the case in the N -body models asthe clusters are slightly elongated due to tidal forces.For the diﬀerences in the most central parts one seesthat the mass density is underestimated for the heaviermass bins while for the lighter mass bins it is overestimated.This could be explained by the fact that these models arepost core collapse, therefore their density proﬁles are slightlydiﬀerent from isothermal models (Lynden-Bell & Eggleton1980). Given that the diﬀerences in the centre are small, onecan see that multimass models are able to describe post-collapse models.The overall agreement between multimass and N -bodymodels for the density proﬁles is consistent with the ﬁnd-ings of Sollima et al. (2017) who ﬁtted Michie–King modelsto observational data of NGC 5466, NGC 6218 and NGC6981: for all three GCs, the multimass models reproducethe observed mass density proﬁles (see their ﬁg. 6). Diﬀer-ences are only found in the outer regions, most likely for thesame reasons as discussed above.When one compares the diﬀerent models in Fig. 7 atthe same dynamical ages it can be seen that models with-out BHs are denser and the stars are found far more con-centrated than in models with BHs. The BHs in the centre‘push’ the lower mass stars out of the core, which results in alarge core radius ( r c ) as well as a larger r h , an eﬀect alreadystudied by Peuten et al. (2016) for the cluster NGC 6101.In the evolution of model N0.3, one can see how the clus-ter changes when all BHs have been lost: the central regionsget eﬃciently populated by the next lighter objects and theresulting mass density proﬁle of the cluster looks as concen-trated as the one from model N0 at that same dynamicalage, leaving no clue about its diminished BH population. For the comparison of the velocity dispersion proﬁles inFig. 8, we used the same snapshots as for the mass den-sity proﬁles. For the calculations of the velocity dispersions,we are using a mass-weighted approach to make the valuescomparable to the values from the multimass model as theyare calculated for the mean mass of each mass bin. The ve-locity dispersion is therefore calculated as: σ k = (cid:80) Ni m i [ v k,i − (cid:104) v k (cid:105) ] (cid:80) Ni m i k = r, θ, φ (15)with (cid:104) v k (cid:105) = (cid:80) Ni m i v k,i (cid:80) Ni m i k = r, θ, φ (16)the mass-weighted mean velocity for each component. Thecalculation of the uncertainties of the binned N -body data was done using the description from Pryor & Meylan (1993)using their equation (12). Again, the results from the best-ﬁtting multimass models are in agreement with the datafrom the N -body models. As with the mass density proﬁles,small diﬀerence can be seen in the outermost regions. In theplot of model N0.3 at 7 . τ rh , , there is no value from the N -body snapshot for the BHs as there is only one BH left,in which case σ is undeﬁned.When comparing the diﬀerent models at the same dy-namical age we ﬁnd that in clusters with BHs, the veloc-ity dispersions for the diﬀerent mass bins are smaller thanin clusters without BHs (see discussion in Section 6.5). Asthe cluster is losing its BH population (see for example theevolution of model N0.3), the velocity dispersions of the dif-ferent mass bins increase to the values seen in model N0which had all its BHs removed before the actual N -bodyevolution, again leaving no hint of the lost BH population.In Section 6.5, we will look again at this relation and discussan explanation for this behaviour. In this section, we consider the anisotropy of the velocities.In Section 5.3.1, we consider the anisotropy proﬁle within thecluster and in Section 5.3.2 we consider the global anisotropyof the cluster as a whole.

The anisotropy parameter β is deﬁned as (Binney &Tremaine 1987): β ≡ − σ σ (17)with σ r the radial velocity dispersion and σ t the tangentialvelocity dispersion. For β <

0, the orbits are tangentiallybiased, for β = 0 they are isotropic, for 0 < β < β = 1 they are radial.In Fig. 9, we compare the anisotropy proﬁles from thebest-ﬁtting multimass models for a selection of mass binsto the anisotropy proﬁles from the N -body snapshots. Asthe β parameter is more aﬀected by random scatter, we hadto bin the data from the N -body snapshots diﬀerently thanin the previous two plots. We varied the number of objectsper bin, such that the average uncertainty in β is ≤ . β uncertainty was always well above 0 . N -body data, we ﬁndthat when the snapshot has some degree of radial anisotropythe multimass models qualitatively reproduce them. Thiscan be seen best with the mass bins from the MS stars.Also diﬀerences between the best-ﬁtting multimass predic-tion and the binned data can be seen at the outer regionsof the cluster. When some of the mass bins are tangentiallyanisotropic the best-ﬁtting model is isotropic as our multi-mass models cannot describe any other kind of anisotropy.Looking at the data from the snapshots itself, wesee that the heaviest mass bins become more radially MNRAS , 1–28 (2017) M. Peuten

Figure 8.

Comparison of the velocity dispersion proﬁles for models N1, N0.3 and N0 at four diﬀerent ages: 2 . τ rh , , 7 . τ rh , , 16 . τ rh , and 26 τ rh , . The points represent the binned N -body data, the thick lines represent the best-ﬁtting multimass models proﬁles and thethin lines represent the results from the walker positions at the last iteration. Red represents the MS stars, cyan – ESs, blue – WDs,pink – NSs and black – BHs. Error bars denote 1 σ uncertainties. MNRAS , 1–28 (2017) esting isothermal models - II. Multimass Figure 9.

Comparison of the anisotropy proﬁles for models N1, N0.3 and N0 at four diﬀerent ages: 2 . τ rh , , 7 . τ rh , , 16 . τ rh , and26 τ rh , . The points represent the binned N -body data, the thick lines represent the best-ﬁtting multimass models proﬁles and the thinlines represent the results from the walker positions at the last iteration. Red represents the MS stars, cyan – ESs, blue – WDs, pink –NSs and black – BHs. Error bars denote 1 σ uncertainties.MNRAS , 1–28 (2017) M. Peuten anisotropic, while the low-mass bins become ﬁrst isotropicand then tangential anisotropic. We will discuss the evolu-tion of the anisotropy further in Section 6.7 when we analysethe best-ﬁtting η parameter. To quantify the global anisotropy we use the parameter κ ,that was introduced by Polyachenko & Shukhman (1981)and is deﬁned as κ = 2 K r K t (18)where K r = 0 . (cid:80) i m i v ,i is the radial component of thekinetic energy and K t = 0 . (cid:80) i m i v ,i the tangential com-ponent. For κ = 1, the models are isotropic, for κ > κ < κ obtained for thebest-ﬁtting multimass models and for the N -body snapshots.For the uncertainties of the N -body data, we used Pois-son statistics. The best-ﬁtting multimass models are able toqualitatively reproduce the overall behaviour of the κ pa-rameter. It can be seen that at later times, the low-massstars are tangentially anisotropic, which cannot be repro-duced by limepy . This also explains why the best-ﬁttingvalue of κ (and β ) from the multimass models does not con-verge to unity immediately (or β (cid:39) κ (and β ) still indicate someradial anisotropy, resulting in a smoother transition from ra-dial anisotropy to isotropy with respect to what is observedin the N -body data. Clusters that are dominated by tangen-tial orbits are therefore, by construction, not well reproducedby the multimass models used in our study (see also Sollimaet al. 2015).Looking at the N -body data, we see that κ of the lowestmass bins typically goes down with time, while κ of theheaviest mass bins goes up. For model N1 with BHs, thischange is faster than for model N0 with no BHs and NSs.We refer the reader to Section 6.7 for a further analysis. We focus here on the best-ﬁt ting parameters resulting fromour ﬁtting procedure. The two model scale parameters canbe computed directly from the N -body data, therefore weuse them to assess the quality of the multimass models. Forthe other ﬁve model parameters, we are analysing their evo-lution to see whether they can give us some further insightsinto the clusters. Furthermore, we also discuss the evolutionof two additional quantities ( r t and κ ) to assess the qualityof the models. In Fig. 11, we plot the M CL from the best-ﬁtting multi-mass model divided by the true cluster mass as measuredin the N -body model for all four N -body models through-out their cluster lifetime. As can be seen in this ﬁgure, thebest-ﬁtting value is always within 1% of the N -body value but almost none of these are consistent within the 1 σ un-certainties: there is some systematic error in the multimassmodels which is not accounted for yet. However, given thatthe diﬀerence is always smaller than 1% this eﬀect is negli-gible. The results are comparable to the single-mass resultsfrom Zocchi et al. (2016) where an accordance within 5% ofthe true values was found.Sollima et al. (2015) found that single-mass models un-derestimated the mass of an N -body system by 50%, de-pending on its dynamical state. There are several diﬀerencesbetween their study and ours which could lead to such dif-ferent results. They are using simulated observations as in-put for their analysis, which aﬀects the recovery of the MF.Shanahan & Gieles (2015) found that approximating mass-segregated clusters by single-component models leads to anunderestimation of the mass by a factor of two or three, es-pecially for metal-rich GCs. Also the models used in Sollimaet al. (2015) have less parameters than the ones used here,as for example they do not incorporate the variable trunca-tion parameter g . Zocchi et al. (2016) showed for single-massmodels that the total mass is better recovered when allowing g to be free. We discuss the eﬀect of g in Section 6.4. The second scale parameter that can be computed from the N -body models is r h : In Fig. 12, we plot r h from the best-ﬁtting multimass models divided by the true value as com-puted from the N -body data, for all four clusters throughouttheir lifetime. Again we ﬁnd good agreement, within a fewpercent. As with M Cl , only a few of the data points showagreement within 1 σ uncertainties. These results are com-parable to the result for the single-mass case by Zocchi et al.(2016) who found an agreement within 7%. Now that we have shown that multimass models can re-produce the most important cluster properties, we focus onanalysing the other ﬁtting parameters. First, we look at thedimensionless central potential W for which the evolutionover the whole lifetime for the four N -body models is plot-ted in Fig. 13. As discussed in Section 2, the W value inmultimass models represents the dimensionless central po-tential of a hypothetical mass group with a mass equal tothe global mean mass. As the global mean mass increasesfrom (0 . ± .

01) M (cid:12) to (0 . ± .

2) M (cid:12) during the evolu-tion of the four N -body models, the W values do not referto the same stellar population and this needs to be kept inmind when comparing W values of diﬀerent N -body modelsand/or at diﬀerent times. Using the central density weightedmean mass instead of the global mean mass has other issues,as can be seen for example in the snapshots of model N0.3at 4 . τ rh , or in model N1 at 26 . τ rh , in Fig. 14. In bothcases, the number of BHs decreases to (almost) zero reduc-ing the central density weighted mean mass more than theglobal mean mass, explaining the more signiﬁcant change ofthe W values in this ﬁgure.It is also possible to use the values of W obtainedfor diﬀerent mass bins for a comparison. As an examplewe consider the ESs mass bin, because not only does the MNRAS , 1–28 (2017) esting isothermal models - II. Multimass Figure 10.

Comparison of the values of the global anisotropy parameter κ as a function of mass of the diﬀerent components for modelsN1, N0.3 and N0 at four diﬀerent ages: 2 . τ rh , , 7 . τ rh , , 16 . τ rh , and 26 τ rh , . The red points represent the N -body data, the thickblack lines represent the best-ﬁtting multimass models values and the thin grey lines represent the results from the walker positions atthe last iteration. Error bars denote 1 σ uncertainties.MNRAS , 1–28 (2017) M. Peuten r e l a t i v e M C l Age [ τ rh,0 ] N1N0.3N0.1N0

Figure 11.

Best-ﬁtting value of the total mass of the cluster, M Cl , divided by the true mass calculated from the N -body snap-shots as a function of time in units of τ rh , . The multimass modelsreproduce the true masses within 1%. Error bars denote 1 σ un-certainties. r e l a t i v e r h Age [ τ rh,0 ] N1N0.3N0.1N0

Figure 12.

Best-ﬁtting half-mass radius r h divided by the truevalue calculated from the N -body snapshots as a function of timein units of τ rh , . The multimass models reproduce the true half-mass radii within 5%. Error bars denote 1 σ uncertainties. average mass of the ESs not vary in our models [ ¯ m ES =(0 . ± . (cid:12) ], but it also represents the objects thatare easiest to observe. In Fig. 15, we plot the evolution ofthe W values of the ESs for all four N -body models overtheir entire lifetime. The uncertainties are estimated by cal-culating the W values of the ESs for the last 10 iterationsof all the walkers and then using the values from the 16thand 84th percentiles as 1 σ uncertainties.Looking at Figs 13 and 15, one can see the eﬀect ofBHs as discussed in Section 5.1: in clusters with BHs, theother stars are pushed outwards and the cluster appears lessconcentrated with a low value of W for the mean mass starsand ESs. While clusters without BHs instead appear muchmore concentrated as the normal stars can occupy the centreand therefore have a higher W value for the mean mass starsand for the ESs. Therefore, clusters with low W value forthe observable stars are much more likely to be hosting a W Age [ τ rh,0 ] N1N0.3N0.1N0

Figure 13.

Best-ﬁtting value of the central dimensionless poten-tial, W , obtained for all four N -body models, as a function oftime in units of τ rh , . Error bars denote 1 σ uncertainties. W , C MM Age [ τ rh,0 ] N1N0.3N0.1N0

Figure 14.

Central dimensionless potential obtained for the best-ﬁtting multimass models when considering the central densityweighted mean mass as the reference mass ¯ m for the four N -body models over time in units of τ rh , . Error bars denote 1 σ uncertainties. BH population than clusters with a high W value (Merrittet al. 2004; Mackey et al. 2008; Peuten et al. 2016). The truncation parameter g provides an indication of theeﬀect of external tides on the stellar system. In Fig. 16, weplot the evolution of g for the four clusters over their life-time. The evolution is similar for all four N -body models: atthe beginning g is around 1 . g = 2) and a King (1966)model ( g = 1). As the clusters ﬁll their Roche volume, thetides interact with the clusters, stripping their outermoststars and thereby making the truncation in energy spacesteeper. This evolution is reﬂected in the truncation param-eter g decreasing as the cluster evolves, converging at theend of the lifetime to a value of around g ≈ .

73 which rep-

MNRAS , 1–28 (2017) esting isothermal models - II. Multimass W , ES Age [ τ rh,0 ] N1N0.3N0.1N0

Figure 15.

Central dimensionless potential of the ESs for thefour N -body models over time in units of τ rh , . Error bars denote1 σ uncertainties. g Age [ τ rh,0 ] N1N0.3N0.1N0

Figure 16.

Best-ﬁtting values of the truncation parameter g ob-tained for the four N -body models over time in units of τ rh , .Error bars denote 1 σ uncertainties. resent a model between a King (1966) and a Woolley (1954)model ( g = 0).The results are comparable to the single-mass modelﬁndings of Zocchi et al. (2016), though they start witha higher truncation parameter, due to the fact that theystart with a smaller initial r h /r J ratio (= 0 .

01) than wedo ( r h /r J = 0 . . r J − r J ) most objects are energetically unbound(Claydon et al. 2017).As before, we see a diﬀerence between the model withBHs (N1) and the models which are mostly BH free (N0,N0.1 and N0.3): at the beginning the value of g decreasesfaster for model N1 than for the other three models andtherefore also converges quicker to its ﬁnal value. This be-haviour in the ﬁrst half of evolution is not too surprisinggiven that r h of that model is roughly twice as large as the δ Age [ τ rh,0 ] N1N0.3N0.1N0

Figure 17.

Best-ﬁtting values of mass segregation parameter δ for all four N -body models over time in units of τ rh , . Error barsdenote 1 σ uncertainties. others, thereby the impact on the steepness of the truncationis stronger (see Section 3.3).McLaughlin & van der Marel (2005) found that Wilsonmodels are equally good or better in describing a sample ofGalactic and extragalactic clusters than King models. Withour results of the evolution of g , one can interpret McLaugh-lin & van der Marel (2005) ﬁndings that clusters with large g are still dynamically expanding towards ﬁlling the Rochevolume (see also Carballo-Bello et al. 2012). In Fig. 17, we plot the evolution of the best-ﬁtting value for δ during the whole lifetime of the four N -body models. Inthe initial stages of their evolution, all four N -body modelsare still in the process of segregation, as they were set upwithout any primordial mass segregation. Over the courseof evolution of the clusters, the value converges to around δ = 0 .

5. This is in accordance with ﬁndings of Sollima et al.(2017), who study mass segregation in observations of GCs.At late stages, there are some snapshots for which the best-ﬁtting value is δ (cid:38) .

5, however the results are compatiblewith 0 . σ . Sollima et al. (2015) found that for someof their late N -body snapshots, the multimass models un-derestimate the amount mass segregation. The same is alsofound for the best-ﬁtting multimass model to the observa-tions of NGC 6218 in Sollima et al. (2017). However, we donot ﬁnd this in these models.To further analyse the behaviour of δ over the courseof the cluster evolution, we additionally plot in Fig. 18 thecentral velocity dispersion for the diﬀerent mass bins fromthe N -body data together with the predicted central velocitydispersion from the best-ﬁtting multimass model (see alsoFig. 8). Given the results from Section 5.2, it is no surprisethat the best-ﬁtting multimass models are able to reproducethe true values.Despite the fact that the best-ﬁtting models have δ ≈ .

5, this does not mean that the multimass models are ina state of energy equipartition, as can be seen in Fig. 18.This was already pointed out by Merritt (1981) and Miocchi

MNRAS , 1–28 (2017) M. Peuten

Figure 18.

Comparison of the central velocity dispersion for the diﬀerent mass bins for Models N1, N0.3 and N0 at four diﬀerentdynamical ages: 2 . τ rh , , 7 . τ rh , , 16 . τ rh , and 26 τ rh , . The red points represent the binned N -body data, the black lines represent thebest-ﬁtting multimass models central velocity dispersion and the thin grey lines represent the results from the walker positions at thelast iteration. The dashed black line shows a σ d (0) ∝ m − / j reference line. For the snapshots with a BH population several additionalmass bins in the high-mass end were included to better show the relation. Error bars denote 1 σ uncertainties.MNRAS , 1–28 (2017) esting isothermal models - II. Multimass (2006) as well as by GZ15. In a mass-segregated multimassmodel the following relation m j s j = m i s i ∀ j, i j (cid:54) = i (19)holds true with s j and s i being the velocity scale of twodiﬀerent mass bins. For the 1D velocity dispersion at thecentre the relation: m j s j = m j σ d,j ∀ j (20)only holds true for W → ∞ . In the N -body models, we arestudying here W (cid:28) ∞ and therefore σ d,j < s j , whichmeans that our multimass models are never in a state ofenergy equipartition despite having δ = 0 . m eq can be in energy equipartition.The value of m eq depends on the mass spectrum of thecluster, and is larger for cluster with a wider mass spec-trum, i.e. when BHs are retained. This trend can also beseen in Fig. 18, where the number of mass bins followingthe σ d (0) ∝ m − / j relation, which are therefore in energyequipartition, is only greater than one for the models with-out BHs. For models with BHs, only the BHs can be inenergy equipartition as they are the only ones which havea mass greater than m eq . This leads to the largest part ofthe other objects in these clusters having a smaller spreadin the velocity dispersions, which is the reason for the re-duced mass segregation in the observable stars in clusterswith BHs.Looking at the results in Fig. 17, we can conclude thatsetting the mass segregation parameter to a ﬁxed value of δ = 0 . δ must therefore be smaller, some-thing also found by Sollima et al. (2015, 2017) for youngclusters. The last two ﬁtting parameters are coupled together as theyboth determine r a ,j for each mass bin as can be seen inequation (3).First we focus on r a . In Fig. 19 we plot r a /r t of thebest-ﬁtting model, for all four N -body models during theirlifetime. If r a /r t (cid:38)

1, the cluster is isotropic (see Section 2).Considering the deﬁnition of r a ,j (equation 3), it follows thateven if r a is well above r t for some mass bins r a ,j can stillbe below r t and therefore these mass bins still show somedegree of radial anisotropy.Looking at Fig. 19, one can see that the model whichretained all its BHs (N1) behaves diﬀerently from the othermodels. Model N1 is only radially anisotropic at the be-ginning of the lifetime and quickly becomes isotropic. Themodel without BHs (N0) loses its radial anisotropy moreslowly and only at the end of its lifetime it becomesisotropic/tangentially anisotropic. The two models in be-tween (N0.3 and N0.1) behave at the beginning diﬀerently:as long as they still have some BHs left their r a value evolvesindependently, but after all BHs are lost their r a value dropsto the r a value of N0 at the same dynamical age and fromthen on essentially follows the r a evolution of model N0. r a / r t Age [ τ rh,0 ] N1N0.3N0.1N0

Figure 19.

Anisotropy radius r a in units of the truncation radius r t as obtained for the best-ﬁtting multimass models for all four N -body models over time in units of τ rh , . Error bars denote 1 σ uncertainties. -4-3-2-1 0 1 0 5 10 15 20 25 30 η Age [ τ rh,0 ] N1N0.3N0.1N0

Figure 20.

Best-ﬁtting anisotropy parameter η for all four N -body models over time in units of τ rh , . Error bars denote 1 σ uncertainties. These results obtained for the BH-free clusters are com-parable, albeit with smaller magnitude, to the results foundfor the single-mass case by Zocchi et al. (2016). They showedthat the anisotropy radius is monotonically decreasing tillthe cluster reaches core collapse after which it is monotoni-cally increasing, and it eventually becomes so large that thecorresponding model is isotropic.Before we discuss the possible physical reasons behindthe evolution of r a , we ﬁrst have a look at the ﬁtting pa-rameter η , which is needed to include a dependence of theanisotropy radius on the mass. The anisotropy parameter η is a novel ﬁtting parameter inmultimass models. In Fig. 20, we plot the values of η for thefour clusters over time. We only plot this for the snapshotsshowing some degree of anisotropy. The most important fea- MNRAS , 1–28 (2017) M. Peuten r t / r J Age [ τ rh,0 ] N1N0.3N0.1N0

Figure 21.

Ratio of the truncation radius r t obtained for thebest-ﬁtting models to the Jacobi radius r J determined from the N -body snapshots for all four models as a function of time inunits of τ rh , . Error bars denote 1 σ uncertainties. ture is that η evolves for all clusters from a value of η ≈ . ≈ − . η ≈ − . η is not surprising given what we al-ready saw in Section 5.3, where we showed that the amountof radial anisotropy is decreasing in the low-mass bins andis increasing in the high-mass bins. This is reﬂected in thedevelopment of η changing from a positive value to a nega-tive one over time. This trend is comparable to what Sollimaet al. (2015) found in their analysis of the W5rh1R8.5 N -body model from Baumgardt & Makino (2003): they foundthat the low-mass stars, which are preferentially located inthe cluster outer regions due to mass segregation, becometangentially anisotropic. The reason for this behaviour isthat interactions occurring in the cluster centre kick starsinto the cluster halo on to radial orbits (Lynden-Bell &Wood 1968; Spitzer & Shull 1975). As stars on radial or-bits reach the cluster boundary with positive velocity, theycan escape the cluster more eﬃciently, thereby depleting thelow-mass population from stars with radial orbits, leavingonly the stars with tangential orbits in the cluster.To test the relevance of η , we rerun ﬁts to model N0 butthis time ﬁxing η = 0 comparable to the original formulationby Gunn & Griﬃn (1979). In these ﬁts, the most obviousdiﬀerence is that with increasing time, and therefore alsowith a higher absolute value of η , the uncertainties of r a increase up to ﬁve times the value recovered in the ﬁts with anon-ﬁxed η value. Therefore, the introduction of η improvesthe ability to describe the data. At the end of our comparison, we look at two quantities onwhich we do not ﬁt but which get computed by the multi-mass models and which can be calculated for the N -bodymodels. In Fig. 21, we plot r t divided by r J as determinedin Section 3.2 for the four N -body models over their whole lifetime. The values of r t and their uncertainties were com-puted using all the walker positions of the last ten iterationsof the MCMC runs. This ﬁgure shows that r t stays within3% of the computed r J and in all but two cases the results areconsistent within their 1 σ uncertainties. The largest discrep-ancies can be seen at the end of the lifetime of the clusters.Compared to the single-mass models in Zocchi et al. (2016)which showed divergence of a factor of two, the multimassmodels are able to reproduce r t accurately. In Fig. 22, we look at the evolution of the global value of κ . Here, we have plotted the comparison between the best-ﬁtting value inferred using the walker positions of the last10 iterations of the MCMC runs from each snapshots to theone calculated from the N -body snapshots directly. Theseare calculated applying equation (18) to all objects in the N -body model. For the uncertainties of the N -body data, weused Poisson statistics. As before the best-ﬁtting multimassmodels qualitatively reproduce the overall trend as long asthe N -body model is radially anisotropic. When the N -bodysnapshots becomes tangentially anisotropic, the best-ﬁttingmultimass models are isotropic. Hence, the multimass κ val-ues still shows some radial anisotropy where the true clusteris already dominated by tangentially anisotropic orbits. Forall models κ < . ± .

25, from which it follows that all mod-els are stable against radial orbit instability as discussed inPolyachenko & Shukhman (1981).

In this study, we assessed the validity of the multimassanisotropic models provided by the limepy software (GZ15)ﬁtting them to N -body models. We ﬁnd that the N -bodymodels are well described by multimass models, a resultwhich is fortunate given the long list of observational stud-ies using multimass models to analyse GCs (see Section 1).Zocchi et al. (2016) showed for the single-mass case thatthe limepy models are able to describe clusters at all evolu-tionary phases. Although the agreement is not perfect, thesystematic diﬀerences are negligible for most applicationsand parameters of interest (see also the discussion in Sec-tion 6.1).Our comparison shows that the best-ﬁtting total clus-ter masses are oﬀ by no more than 1% from the true valueas computed from the N -body snapshots. The best-ﬁttingcluster half-mass radius is reproduced within 5% and thetruncation radius is reproduced within 3%.We ﬁnd that the mass density and velocity dispersionproﬁles of the diﬀerent mass bins are well reproduced bythe multimass models. If the N -body snapshot is radiallyanisotropic then the multimass models are generally able toreproduce it.We show that in the N -body models, regardless of initialBHs and NSs retention, the truncation parameter g evolvesfrom roughly 1 . .

7. The general trend can be ex-plained by the tidal eﬀects stripping the loosely bound stars.We ﬁnd that the best-ﬁtting mass segregation parameter δ converges to a value close to 0 . N -body models,which is the value used in the original formulation by Da MNRAS , 1–28 (2017) esting isothermal models - II. Multimass κ Age [ τ rh,0 ] N1 Best-fit modelsN-body snapshots κ Age [ τ rh,0 ] N0.3 κ Age [ τ rh,0 ] N0.1 κ Age [ τ rh,0 ] N0 Figure 22.

Comparison of the values of the global anisotropy pa-rameter κ for the four diﬀerent N -body models over their wholelifetime. The black dashed lines represent the best-ﬁtting multi-mass estimate, and the red solid lines represent the true valuesdirectly calculated from the N -body snapshots. Costa & Freeman (1976). Only for young clusters which arenot yet mass segregated is the best-ﬁtting value smaller andfor models with BHs it is ∼ . η parameter shows that theanisotropy radius is mass dependent and that this massdependence changes in our N -body models over time from η = 0 . η = − . η = − . η is another relevant parame-ter when analysing radial anisotropy with multimass models.Furthermore, we ﬁnd that clusters with a BH population canbe tangentially anisotropic for most of their lifetime.The W parameter for the observable ESs is lower forthe clusters with BHs than for the clusters without BHs.Therefore, clusters which still harbour a stellar-mass BHpopulation should appear less dense when looking at theobservable stars. N -body simulation N0.1, which loses itsBHs within its ﬁrst τ rh , , does not show any strong diﬀer-ences to simulation N0, despite having a population of NSs.The inﬂuence of the NSs on a cluster is therefore negligible,compared to the impact stellar-mass BHs have on a cluster.We conclude that the limepy multimass models are anadequate tool to study the global properties of GCs, as theresults from the comparison with N -body models show agood agreement with their properties inferred from multi-mass models. ACKNOWLEDGEMENTS

This research was done as part of the

Gaia

Challenge ,whose goal is to compare and improve mass modelling tech-niques in preparation of data releases of the ESA Gaia mis-sion. We thank the organizers of the Third Gaia ChallengeWorkshop (Barcelona, 2015) for a very productive meeting.We thank the anonymous referee for constructive sugges-tions. We are grateful to Anna Lisa Varri and Antonio Sol-lima for interesting discussions. We are grateful to SverreAarseth and Keigo Nitadori for making nbody6 publiclyavailable, and to Dan Foreman-Mackey for providing the emcee software and for maintaining the online documenta-tion; we also thank Mr. Dave Munro of the University of Sur-rey for hardware and software support. MG acknowledges ﬁ-nancial support from the Royal Society (University ResearchFellowship) and AZ acknowledges ﬁnancial support from theRoyal Society (Newton International Fellowship). MP, MGand AZ acknowledge the European Research Council (ERC-StG-335936, CLUSTERS). VHB acknowledges support fromthe Radboud Excellence Initiative Fellowship.

REFERENCES

Aarseth S. J., 2003, Gravitational N-Body Simulations. Cam-bridge University Press, CambridgeBaumgardt H., 2001, MNRAS, 325, 1323 http://astrowiki.ph.surrey.ac.uk/dokuwiki/doku.php?id=start MNRAS , 1–28 (2017) M. Peuten

Baumgardt H., 2017, MNRAS, 464, 2174Baumgardt H., Makino J., 2003, MNRAS, 340, 227Baumgardt H., Cˆot´e P., Hilker M., Rejkuba M., Mieske S., Djor-govski S. G., Stetson P., 2010, in de Grijs R., L´epine J. R. D.,eds, IAU Symposium Vol. 266, IAU Symposium. pp 365–365,doi:10.1017/S174392130999130XBeccari G., Pasquato M., De Marchi G., Dalessandro E., TrentiM., Gill M., 2010, ApJ, 713, 194Beccari G., Dalessandro E., Lanzoni B., Ferraro F. R., BellazziniM., Sollima A., 2015, ApJ, 814, 144Bianchini P., van de Ven G., Norris M. A., Schinnerer E., VarriA. L., 2016, MNRAS, 458, 3644Binney J., McMillan P., 2011, MNRAS, 413, 1889Binney J., Tremaine S., 1987, Galactic dynamics. Princeton Univ.Press, Princeton, NJCarballo-Bello J. A., Gieles M., Sollima A., Koposov S., Mart´ınez-Delgado D., Pe˜narrubia J., 2012, MNRAS, 419, 14Chernoﬀ D. F., Weinberg M. D., 1990, ApJ, 351, 121Claydon I., Gieles M., Zocchi A., 2017, MNRAS, 466, 3937Cohn H., 1980, ApJ, 242, 765Da Costa G. S., Freeman K. C., 1976, ApJ, 206, 128Davoust E., 1977, A&A, 61, 391Eddington A. S., 1915, MNRAS, 75, 366Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013,PASP, 125, 306Fukushige T., Heggie D. C., 2000, MNRAS, 318, 753Gaia Collaboration et al., 2016, A&A, 595, A1Gieles M., Zocchi A., 2015, MNRAS, 454, 576Gieles M., Baumgardt H., Heggie D. C., Lamers H. J. G. L. M.,2010, MNRAS, 408, L16Gieles M., Heggie D. C., Zhao H., 2011, MNRAS, 413, 2509Giersz M., Heggie D. C., 1996, MNRAS, 279, 1037Giersz M., Heggie D. C., 1997, MNRAS, 286, 709Giersz M., Heggie D. C., Hurley J. R., Hypki A., 2013, MNRAS,431, 2184Gomez-Leyton Y. J., Velazquez L., 2014, Journal of StatisticalMechanics: Theory and Experiment, 4, 6Goodman J., Hut P., 1989, Nature, 339, 40Goodman J., Weare J., 2010, Communications in Applied Math-ematics and Computational science, 5, 65Gratton R., Sneden C., Carretta E., 2004, ARA&A, 42, 385Gunn J. E., Griﬃn R. F., 1979, AJ, 84, 752Heggie D. C., 1975, MNRAS, 173, 729Heggie D. C., 2014, MNRAS, 445, 3435Heggie D. C., Giersz M., 2008, MNRAS, 389, 1858Heggie D. C., Trenti M., Hut P., 2006, MNRAS, 368, 677H´enon M., 1959, Annales d’Astrophysique, 22, 126H´enon M., 1961, Annales d’Astrophysique, 24, 369H´enon M., 1965, Annales d’Astrophysique, 28, 62Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543Hurley J. R., Pols O. R., Aarseth S. J., Tout C. A., 2005, MNRAS,363, 293Hut P., et al., 1992, PASP, 104, 981Ibata R., Nipoti C., Sollima A., Bellazzini M., Chapman S. C.,Dalessandro E., 2013, MNRAS, 428, 3648Illingworth G., King I. R., 1977, ApJ, 218, L109Johnston K. V., Sigurdsson S., Hernquist L., 1999, MNRAS, 302,771King I. R., 1966, AJ, 71, 64Kroupa P., 2001, MNRAS, 322, 231K¨upper A. H. W., Kroupa P., Baumgardt H., Heggie D. C., 2010,MNRAS, 401, 105Lee H. M., Goodman J., 1995, ApJ, 443, 109Lee H. M., Ostriker J. P., 1987, ApJ, 322, 123Lupton R. H., Gunn J. E., Griﬃn R. F., 1987, AJ, 93, 1114Lynden-Bell D., Eggleton P. P., 1980, MNRAS, 191, 483Lynden-Bell D., Wood R., 1968, MNRAS, 138, 495 Mackey A. D., Wilkinson M. I., Davies M. B., Gilmore G. F.,2008, MNRAS, 386, 65Mandel I., 2016, MNRAS, 456, 578McLaughlin D. E., 2003, in Piotto G., Meylan G., DjorgovskiS. G., Riello M., eds, Astronomical Society of the Paciﬁc Con-ference Series Vol. 296, New Horizons in Globular Cluster As-tronomy. p. 101McLaughlin D. E., van der Marel R. P., 2005, ApJS, 161, 304Merritt D., 1981, AJ, 86, 318Merritt D., Piatek S., Portegies Zwart S., Hemsendorf M., 2004,ApJ, 608, L25Meylan G., 1987, A&A, 184, 144Meylan G., Heggie D. C., 1997, A&ARv, 8, 1Meylan G., Mayor M., 1991, A&A, 250, 113Meylan G., Mayor M., Duquennoy A., Dubath P., 1995, A&A,303, 761Michie R. W., 1963, MNRAS, 125, 127Milone A. P., et al., 2014, ApJ, 785, 21Miocchi P., 2006, MNRAS, 366, 227Nitadori K., Aarseth S. J., 2012, MNRAS, 424, 545Oort J. H., van Herk G., 1959, Bull. Astron. Inst. Netherlands,14, 299Paust N. E. Q., et al., 2010, AJ, 139, 476Peuten M., Zocchi A., Gieles M., Gualandris A., H´enault-BrunetV., 2016, MNRAS, 462, 2333Piotto G., Zoccali M., 1999, A&A, 345, 485Polyachenko V. L., Shukhman I. G., 1981, Soviet Ast., 25, 533Pryor C., Meylan G., 1993, in Djorgovski S. G., Meylan G., eds,Astronomical Society of the Paciﬁc Conference Series Vol. 50,Structure and Dynamics of Globular Clusters. p. 357Pryor C., Hartwick F. D. A., McClure R. D., Fletcher J. M.,Kormendy J., 1986, AJ, 91, 546Renaud F., Agertz O., Gieles M., 2017, MNRAS, 465, 3622Repetto S., Davies M. B., Sigurdsson S., 2012, MNRAS, 425, 2799Richer H. B., Fahlman G. G., 1989, ApJ, 339, 178Richer H. B., et al., 2004, AJ, 127, 2771Rieder S., Ishiyama T., Langelaan P., Makino J., McMillanS. L. W., Portegies Zwart S., 2013, MNRAS, 436, 3695Shanahan R. L., Gieles M., 2015, MNRAS, 448, L94Sollima A., Bellazzini M., Lee J.-W., 2012, ApJ, 755, 156Sollima A., Baumgardt H., Zocchi A., Balbinot E., Gieles M.,H´enault-Brunet V., Varri A. L., 2015, MNRAS, 451, 2185Sollima A., Dalessandro E., Beccari G., Pallanca C., 2017, MN-RAS, 464, 3871Sosin C., 1997, AJ, 114, 1517Spitzer Jr. L., Shull J. M., 1975, ApJ, 201, 773Takahashi K., Lee H. M., 2000, MNRAS, 316, 671Takahashi K., Portegies Zwart S. F., 2000, ApJ, 535, 759Trenti M., van der Marel R., 2013, MNRAS, 435, 3272Trenti M., Heggie D. C., Hut P., 2007, MNRAS, 374, 344Trenti M., Vesperini E., Pasquato M., 2010, ApJ, 708, 1598Wang L., Spurzem R., Aarseth S., Giersz M., Askar A., BerczikP., Naab T., Kouwenhoven R. S. M. B. N., 2016, MNRAS,458, 1450Wilson C. P., 1975, AJ, 80, 175Woolley R. V. D. R., 1954, MNRAS, 114, 191Zocchi A., Gieles M., H´enault-Brunet V., Varri A. L., 2016, MN-RAS, 462, 696Zonoozi A. H., K¨upper A. H. W., Baumgardt H., Haghi H.,Kroupa P., Hilker M., 2011, MNRAS, 411, 1989

APPENDIX A: N -BODY DATA AND MCMCRESULTS MNRAS , 1–28 (2017) esting isothermal models - II. Multimass Table A1.

Properties of the snapshots from the N -body model N1 with 100% initial BH and NS retention. We list the age in units ofthe initial half-mass relaxation time, the total bound mass in M (cid:12) , the half-mass radius in pc, the number of BHs in the cluster, thenumber of NSs in the cluster and the Jacobi radius r J in pc calculated as in Section 3.2. For the bound mass and the number of BHsand NSs, we also give the percentage relative to the initial values in brackets.Age M CL Half-mass radius Number of BHs Number of NSs r J ( τ rh , ) ( M (cid:12) ) (pc) (pc)0 37353 (100%) 2 .

06 121 (100%) 632 (100%) 30 . . .

03 119 (56%) 551 (87%) 27 . . .

03 85 (40%) 518 (82%) 26 . . .

56 59 (28%) 503 (80%) 25 . . .

58 55 (26%) 490 (78%) 24 . . .

54 43 (20%) 474 (75%) 22 . . .

27 35 (17%) 461 (73%) 21 . . .

71 26 (12%) 439 (69%) 19 . . .

98 18 (8 . . . .

09 12 (5 . . . . .

92 6 (2 . . . . .

58 0 (0 . . Table A2.

Properties of the snapshots from the N -body model N0.3 with 33% initial BH and NS retention. We list the age in unitsof the initial half-mass relaxation time, the total bound mass in M (cid:12) , the half-mass radius in pc, the number of BHs in the cluster, thenumber of NSs in the cluster and the Jacobi radius r J in pc calculated as in Section 3.2. For the bound mass and the number of BHsand NSs, we also give the percentage relative to the initial values in brackets.Age M CL Half-mass radius Number of BHs Number of NSs r J ( τ rh , ) ( M (cid:12) ) (pc) (pc)0 34404 (100%) 2 .

06 66 (100%) 211 (100%) 29 . . .

07 22 (33%) 199 (94%) 28 . . .

13 8 (12%) 189 (90%) 27 . . .

12 1 (1 . . . .

27 0 (0%) 150 (71%) 25 . . .

57 0 (0%) 119 (56%) 24 . . .

57 0 (0%) 94 (45%) 23 . . .

61 0 (0%) 75 (36%) 22 . . .

43 0 (0%) 59 (28%) 20 . . .

28 0 (0%) 46 (22%) 18 . . .

88 0 (0%) 37 (18%) 16 . . .

61 0 (0%) 26 (12%) 13 . . . .

89 0 (0%) 15 (7 . . . . .

08 0 (0%) 7 (3 . . , 1–28 (2017) M. Peuten

Table A3.

Properties of the snapshots from the N -body model N0.1 with 10% initial BH and NS retention. We list the age in unitsof the initial half-mass relaxation time, the total bound mass in M (cid:12) , the half-mass radius in pc, the number of BHs in the cluster, thenumber of NSs in the cluster and the Jacobi radius r J in pc calculated as in Section 3.2. For the bound mass and the number of BHsand NSs, we also give the percentage relative to the initial values in brackets.Age M CL Half-mass radius Number of BHs Number of NSs r J ( τ rh , ) ( M (cid:12) ) (pc) (pc)0 33476 (100%) 2 .

04 22 (100%) 56 (100%) 29 . . .

46 5 (23%) 54 (96%) 28 . . .

76 0 (0%) 50 (89%) 27 . . .

27 0 (0%) 24 (43%) 26 . . .

52 0 (0%) 19 (34%) 25 . . .

60 0 (0%) 12 (21%) 23 . . .

65 0 (0%) 7 (13%) 22 . . .

52 0 (0%) 5 (8 . . . .

42 0 (0%) 3 (5 . . . .

11 0 (0%) 1 (1 . . . .

75 0 (0%) 0 (0%) 14 . . . .

29 0 (0%) 0 (0%) 11 . . . .

48 0 (0%) 0 (0%) 8 . Table A4.

Properties of the snapshots from the N -body model N0 with no initial BH and NS retention. We list the age in units of theinitial half-mass relaxation time, the total bound mass in M (cid:12) , the half-mass radius in pc and the Jacobi radius r J in pc calculated asin Section 3.2. For the bound mass we also give the percentage relative to the initial value in brackets.Age M Cl Half-mass radius r J ( τ rh , ) ( M (cid:12) ) (pc) (pc)0 33042 (100%) 2 .

04 28 . . .

26 28 . . .

93 27 . . .

25 25 . . .

46 24 . . .

53 23 . . .

56 22 . . .

47 20 . . .

18 18 . . .

89 16 . . .

55 14 . . . .

12 11 . . . .

06 7 .

18 MNRAS , 1–28 (2017) e s t i n g i s o t h e r m a l m od e l s - II . M u l t i m a ss Table A5.

Mass bins of the diﬀerent snapshots of N -body model N1. We list the age in units of the initial half-mass relaxation time and for each mass bin the total mass M j and themean mass m j in units of M (cid:12) . There are a total of 11 mass bins: ﬁve for the MSs, one for the ESs, three for the WDs and one each for the NS and BHs. Age MS1 MS2 MS3 MS4 MS5 ES WD1 WD2 WD3 NS BH

Mj mj Mj mj Mj mj Mj mj Mj mj Mj mj Mj mj Mj mj Mj mj Mj mj Mj mj ( τ rh ,

0) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) )2 . .

13 3288 0 . .

32 6155 0 . .

72 185 0 .

82 3650 0 . .

78 1521 1 .

08 841 1 .

53 1504 12 . . .

13 2702 0 . .

32 5287 0 . .

72 164 0 .

82 3152 0 . .

78 1391 1 .

08 789 1 .

52 955 11 . . .

13 2182 0 . .

32 4572 0 . .

72 152 0 .

82 2775 0 . .

78 1299 1 .

09 766 1 .

52 596 10 . . .

13 1722 0 . .

32 3943 0 . .

72 138 0 .

82 2475 0 . .

78 1230 1 .

09 748 1 .

53 538 9 . . .

13 1280 0 . .

32 3283 0 . .

72 118 0 .

82 2104 0 . .

78 1154 1 .

09 725 1 .

53 379 8 . . .

13 874 0 . .

32 2549 0 . .

72 104 0 .

83 1728 0 . .

79 1055 1 .

09 706 1 .

53 280 7 . . .

13 528 0 . .

32 1818 0 .

51 1519 0 .

72 82 0 .

83 1299 0 . .

79 957 1 . .

54 192 7 . . .

13 269 0 . .

33 1172 0 .

51 1119 0 .

72 62 0 .

83 881 0 . .

79 831 1 . .

54 119 6 . . . .

13 98 0 . .

33 611 0 .

51 728 0 .

73 41 . .

83 496 0 .

61 541 0 .

79 647 1 .

11 574 1 .

56 80 6 . . .

99 0 .

13 22 . .

21 63 0 .

33 213 0 .

52 325 0 .

73 15 . .

83 204 0 .

61 272 0 . .

13 474 1 .

57 30 . . . .

12 0 .

12 0 .

24 0 .

24 4 .

52 0 .

35 14 . .

53 50 0 .

74 2 .

51 0 .

84 15 . .

62 47 . .

81 119 1 .

18 243 1 .

63 0 0

Table A6.

Mass bins of the diﬀerent snapshots of N -body model N0.3. We list the age in units of the initial half-mass relaxation time and for each mass bin the total mass M j and themean mass m j in units of M (cid:12) . There are a total of 11 mass bins: ﬁve for the MSs, one for the ESs, three for the WDs and one each for the NS and BHs. Age MS1 MS2 MS3 MS4 MS5 ES WD1 WD2 WD3 NS BH

Mj m j Mj m j Mj m j Mj m j Mj m j Mj m j Mj m j Mj m j Mj m j Mj m j Mj m j( τ rh ,

0) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) )2 . .

13 3609 0 . .

32 6732 0 . .

72 200 0 .

81 3975 0 . .

78 1622 1 .

08 302 1 .

52 224 10 . . .

13 3276 0 . .

32 6340 0 . .

72 196 0 .

82 3805 0 . .

78 1588 1 .

08 285 1 .

51 53 6 . . .

13 2869 0 . .

32 5889 0 . .

72 184 0 .

82 3603 0 . .

78 1546 1 .

08 269 1 .

50 6 .

31 6 . . .

13 2386 0 . .

32 5356 0 . .

72 179 0 .

82 3329 0 . .

78 1478 1 .

08 218 1 .

45 0 011 . .

13 1782 0 . .

32 4691 0 . .

72 164 0 .

82 3005 0 . .

78 1354 1 .

07 168 1 .

41 0 014 . .

13 1406 0 . .

32 4003 0 . .

72 155 0 .

82 2675 0 . .

78 1244 1 .

06 130 1 .

38 0 016 . .

13 990 0 . .

32 3308 0 . .

72 143 0 .

82 2274 0 .

61 1772 0 .

79 1140 1 .

05 102 1 .

36 0 018 . .

13 622 0 . .

32 2595 0 .

51 2263 0 .

72 126 0 .

82 1859 0 .

61 1559 0 .

79 1024 1 .

05 80 1 .

35 0 021 . .

13 341 0 . .

33 1862 0 .

51 1847 0 .

73 104 0 .

82 1390 0 .

61 1329 0 .

79 904 1 .

04 61 1 .

34 0 023 . .

13 144 0 . .

33 1169 0 .

52 1394 0 .

73 83 0 .

82 949 0 .

61 1055 0 . .

04 49 . .

33 0 025 . . .

13 45 . . .

33 587 0 .

52 892 0 .

73 55 0 .

81 533 0 .

61 729 0 . .

04 43 . .

32 0 028 . .

92 0 .

13 5 .

93 0 . . .

35 148 0 .

53 367 0 .

74 20 . .

83 161 0 .

61 349 0 .

81 393 1 .

04 19 . .

31 0 030 . .

17 0 .

17 0 .

39 0 .

39 0 0 6 . .

55 38 . .

74 3 .

32 0 .

83 8 .

71 0 .

62 69 0 .

85 133 1 .

06 9 .

18 1 .

31 0 0 M N R A S , ( ) M . P e u t e n Table A7.

Mass bins of the diﬀerent snapshots of N -body model N0.1. We list the age in units of the initial half-mass relaxation time and for each mass bin the total mass M j and themean mass m j in units of M (cid:12) . There are a total of 11 mass bins: ﬁve for the MSs, one for the ESs, three for the WDs and one each for the NS and BHs. Age MS1 MS2 MS3 MS4 MS5 ES WD1 WD2 WD3 NS BH

Mj mj Mj mj Mj mj Mj mj Mj mj Mj mj Mj mj Mj mj Mj mj Mj mj Mj mj ( τ rh ,

0) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) )2 . .

13 3658 0 . .

32 6844 0 . .

72 204 0 .

81 4039 0 . .

78 1644 1 .

08 82 1 .

52 53 10 . . .

13 3220 0 . .

32 6374 0 . .

72 200 0 .

82 3827 0 . .

78 1577 1 .

08 75 1 .

50 0 07 . .

13 2683 0 . .

32 5703 0 . .

72 186 0 .

81 3515 0 . .

78 1334 1 .

06 33 1 .

36 0 09 . .

13 2142 0 . .

32 5015 0 . .

72 180 0 .

82 3183 0 . .

78 1185 1 .

05 26 1 .

35 0 011 . .

13 1652 0 . .

32 4323 0 . .

72 168 0 .

82 2824 0 . .

78 1065 1 .

04 16 . .

36 0 014 . .

13 1207 0 . .

32 3606 0 . .

72 158 0 .

82 2455 0 . .

78 935 1 .

03 9 .

26 1 .

33 0 016 . .

13 787 0 . .

32 2845 0 .

51 2345 0 .

72 147 0 .

82 2064 0 .

61 1635 0 .

79 811 1 .

02 6 .

51 1 . . .

13 456 0 . .

32 2148 0 .

51 2025 0 .

73 130 0 .

82 1651 0 .

61 1414 0 .

79 711 1 .

01 3 .

91 1 . . .

13 227 0 . .

33 1411 0 .

51 1600 0 .

73 113 0 .

82 1179 0 .

61 1169 0 .

79 592 1 .

01 1 . . . . .

13 81 0 . .

33 745 0 .

52 1121 0 .

73 89 0 .

82 724 0 .

61 871 0 .

79 481 1 . . .

55 0 .

13 14 . . .

33 291 0 .

53 596 0 .

74 55 0 .

83 305 0 .

61 511 0 . .

01 0 0 0 028 . .

14 0 .

14 0 .

48 0 .

24 3 .

68 0 .

33 53 0 .

55 196 0 .

75 21 . .

83 57 0 .

61 192 0 .

82 155 1 .

01 0 0 0 0

Table A8.

Mass bins of the diﬀerent snapshots of N -body model N0. We list the age in units of the initial half-mass relaxation time and for each mass bin the total mass M j and themean mass m j in units of M (cid:12) . There are a total of nine mass bins: ﬁve for the MSs, one for the ESs and three for the WDs. Age MS1 MS2 MS3 MS4 MS5 ES WD1 WD2 WD3

Mj mj Mj mj Mj mj Mj mj Mj mj Mj mj Mj mj Mj mj Mj mj ( τ rh ,

0) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) ) ( M (cid:12) )2 . .

13 3613 0 . .

32 6817 0 . .

72 210 0 .

82 4016 0 . .

78 1659 1 . . .

13 3091 0 . .

32 6174 0 . .

72 197 0 .

81 3703 0 . .

78 1356 1 . . .

13 2543 0 . .

32 5515 0 . .

72 185 0 .

81 3385 0 . .

78 1135 1 . . .

13 2055 0 . .

32 4874 0 . .

72 174 0 .

81 3074 0 . .

78 987 1 . . .

13 1589 0 . .

32 4210 0 . .

72 159 0 .

81 2708 0 . .

78 867 1 . . .

13 1147 0 . .

32 3495 0 .

51 2667 0 .

72 140 0 .

81 2329 0 . .

78 745 1 . . .

13 747 0 . .

32 2763 0 .

51 2336 0 .

72 123 0 .

82 1945 0 .

61 1557 0 .

78 636 1 . . .

13 432 0 . .

32 2045 0 .

51 1961 0 .

73 107 0 .

82 1550 0 .

61 1362 0 .

79 555 0 . . .

13 196 0 . .

33 1316 0 .

51 1518 0 .

73 87 0 .

82 1098 0 .

61 1086 0 .

79 463 0 . . . .

13 68 0 . .

33 707 0 .

52 1045 0 .

73 66 0 .

82 681 0 .

61 795 0 .

79 353 0 . . .

32 0 .

13 13 0 .

21 52 0 .

34 248 0 .

53 578 0 .

74 43 0 .

83 287 0 .

61 473 0 . . . .

11 0 .

11 0 .

61 0 . .

31 0 .

33 25 . .

54 139 0 .

76 10 . .

84 29 0 .

62 138 0 .

82 111 0 . M N R A S , (2017

82 111 0 . M N R A S , (2017 ) e s t i n g i s o t h e r m a l m od e l s - II . M u l t i m a ss Table A9.

Results from the MCMC ﬁtting process for the snapshots from the N -body model N1 with initial 100% BH and NS retention. We list here the age of each snapshot in unitsof the initial half-mass relaxation time, the total cluster mass in M (cid:12) , the half-mass radius in pc, the dimensionless central concentration for the global mean mass, the central meanmass, and the ES stars, the truncation parameter g , the mass segregation parameter δ , the anisotropy radius for the global mean mass and the central mean mass in pc, the anisotropyparameter η and the truncation radius r t in pc. All parameters except the dimensionless central concentration for the global mean mass, and the ES stars, the anisotropy radius for theglobal mean mass and the truncation radius r t , are ﬁtting parameters; the other values were obtained from the multimass models of the 10 last walker positions of each MCMC chain.The median of the marginalized posterior distribution of each parameter is used to estimate its best-ﬁtting value and the 16th and 84th percentiles as proxy for the 1 σ uncertainties.Age M Cl r h W W , CMM W , ES g δ r a r a , CMM η r t ( τ rh , ) ( M (cid:12) ) (pc) (pc) (pc) (pc)2 . +14 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +13 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . − . +0 . − . . +0 . − . . +11 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . +2 . − . . +0 . − . − . +0 . − . . +0 . − . . +9 . − . . +0 . − . . +0 . − . +1 . − . . +0 . − . . +0 . − . . +0 . − . + ∞− +4770 − − . +0 . − . . +0 . − . . +8 . − . . +0 . − . . +0 . − . +1 . − . . +0 . − . . +0 . − . . +0 . − . + ∞− +208 − − . +0 . − . . +0 . − . . +7 . − . . +0 . − . . +0 . − . +1 . − . . +0 . − . . +0 . − . . +0 . − . – 4416 +2814 − − +4 − . +0 . − . . +6 . − . . +0 . − . . +0 . − . +1 . − . . +0 . − . . +0 . − . . +0 . − . – 3740 +2365 − − +4 − . +0 . − . . +6 . − . . +0 . − . . +0 . − . +2 . − . . +0 . − . . +0 . − . . +0 . − . – 1017 +2149 − − +4 − . +0 . − . . +5 . − . . +0 . − . . +0 . − . +3 . − . . +0 . − . . +0 . − . . +0 . − . – 3041 +1843 − − +4 − . +0 . − . . +4 . − . . +0 . − . . +0 . − . +4 . − . . +0 . − . . +0 . − . . +0 . − . – 874 +600 − − +4 − . +0 . − . . +2 . − . . +0 . − . +1 . − . . +1 . − . . +1 . − . . +0 . − . . +0 . − . +266 − +12 − − +4 − . +0 . − . Table A10.

Results from the MCMC ﬁtting process for the snapshots from the N -body model N0.3 with initial 33% BH and NS retention. We list here the age of each snapshot in unitsof the initial half-mass relaxation time, the total cluster mass in M (cid:12) , the half-mass radius in pc, the dimensionless central concentration for the global mean mass, the central meanmass, and the ES stars, the truncation parameter g , the mass segregation parameter δ , the anisotropy radius for the global mean mass and the central mean mass in pc, the anisotropyparameter η and the truncation radius r t in pc. All parameters except the dimensionless central concentration for the global mean mass, and the ES stars, the anisotropy radius for theglobal mean mass and the truncation radius r t , are ﬁtting parameters; the other values were obtained from the multimass models of the 10 last walker positions of each MCMC chain.The median of the marginalized posterior distribution of each parameter is used to estimate its best-ﬁtting value and the 16th and 84th percentiles as proxy for the 1 σ uncertainties.Age M Cl r h W W , CMM W , ES g δ r a r a , CMM η r t ( τ rh , ) ( M (cid:12) ) (pc) (pc) (pc) (pc)2 . +13 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +2 . − . . +0 . − . . +0 . − . . +14 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +14 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . − . +0 . − . . +0 . − . . +12 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . − . +0 . − . . +0 . − . . +11 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . − . +0 . − . . +0 . − . . +9 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . − . +0 . − . . +0 . − . . +9 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . − . +0 . − . . +0 . − . . +8 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . − . +0 . − . . +0 . − . . +7 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . +2 − . +0 . − . − . +0 . − . . +0 . − . . +7 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . +4 − . +0 . − . − . +0 . − . . +0 . − . . +6 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . +7 − . +1 . − . − . +0 . − . . +0 . − . . +3 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . +8346 − +696 − − . +4 . − . . +0 . − . . +2 − . +0 . − . . +0 . − . . +0 . − . . +1 . − . . +0 . − . . +0 . − . +1444 − +504 − − . +4 . − . . +0 . − . M N R A S , ( ) M . P e u t e n Table A11.

Results from the MCMC ﬁtting process for the snapshots from the N -body model N0.1 with initial 10% BH and NS retention. We list here the age of each snapshot in unitsof the initial half-mass relaxation time, the total cluster mass in M (cid:12) , the half-mass radius in pc, the dimensionless central concentration for the global mean mass, the central meanmass, and the ES stars, the truncation parameter g , the mass segregation parameter δ , the anisotropy radius for the global mean mass and the central mean mass in pc, the anisotropyparameter η and the truncation radius r t in pc. All parameters except the dimensionless central concentration for the global mean mass, and the ES stars, the anisotropy radius for theglobal mean mass and the truncation radius r t , are ﬁtting parameters; the other values were obtained from the multimass models of the 10 last walker positions of each MCMC chain.The median of the marginalized posterior distribution of each parameter is used to estimate its best-ﬁtting value and the 16th and 84th percentiles as proxy for the 1 σ uncertainties.Age M Cl r h W W , CMM W , ES g δ r a r a , CMM η r t ( τ rh , ) ( M (cid:12) ) (pc) (pc) (pc) (pc)2 . +15 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . +1 − . +0 . − . . +0 . − . . +14 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +11 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . − . +0 . − . . +0 . − . . +11 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . − . +0 . − . . +0 . − . . +10 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . − . +0 . − . . +0 . − . . +9 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . − . +0 . − . . +0 . − . . +8 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . − . +0 . − . . +0 . − . . +8 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . − . +0 . − . . +0 . − . . +7 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . +3 − . +0 . − . − . +0 . − . . +0 . − . . +6 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . +8138 − +491 − − . +4 − . +0 . − . . +4 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . +2822 − +305 − − . +4 − . +0 . − . . +2 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . +1719 − +589 − − . +4 − . +0 . − . Table A12.

Results from the MCMC ﬁtting process for the snapshots from the N -body model N0 with no initial BH and NS retention. We list here the age of each snapshot in unitsof the initial half-mass relaxation time, the total cluster mass in M (cid:12) , the half-mass radius in pc, the dimensionless central concentration for the global mean mass, the central meanmass, and the ES stars, the truncation parameter g , the mass segregation parameter δ , the anisotropy radius for the global mean mass and the central mean mass in pc, the anisotropyparameter η and the truncation radius r t in pc. All parameters except the dimensionless central concentration for the global mean mass, and the ES stars, the anisotropy radius for theglobal mean mass and the truncation radius r t , are ﬁtting parameters; the other values were obtained from the multimass models of the 10 last walker positions of each MCMC chain.The median of the marginalized posterior distribution of each parameter is used to estimate its best-ﬁtting value and the 16th and 84th percentiles as proxy for the 1 σ uncertainties.Age M Cl r h W W , CMM W , ES g δ r a r a , CMM η r t ( τ rh , ) ( M (cid:12) ) (pc) (pc) (pc) (pc)2 . +13 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +13 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +12 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . − . +0 . − . . +0 . − . . +11 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . − . +0 . − . . +0 . − . . +10 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . − . +0 . − . . +0 . − . . +9 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . − . +0 . − . . +0 . − . . +8 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . − . +0 . − . . +0 . − . . +7 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . +1 − . +0 . − . − . +0 . − . . +0 . − . . +6 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . +2 − . +0 . − . − . +0 . − . . +0 . − . . +5 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . +11 − +2 − − . +0 . − . . +0 . − . . +4 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . +3098 − +390 − − +4 − . +0 . − . . +2 − . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . +17 − +6 − − +4 − . +0 . − . M N R A S , (2017