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 file 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 different retention fractions for theblack holes and neutron stars evolving in a tidal field. We find that multimass modelssuccessfully reproduce the density and velocity dispersion profiles of the different masscomponents in all evolutionary phases and for different remnants retention. We furtheruse these results to study the evolution of global model parameters. We find 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 find 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 fitting 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 configurations as the result of two-body relaxation, GCs are complex systems to model. Theyconsist of stars and stellar remnants with different 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 field that influences their evolu-tion and present-day properties (Chernoff & 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 first N -body simulation of an open clus-ter, Zonoozi et al. (2011) modelled a low-mass GC and finallyHeggie (2014) and Wang et al. (2016) presented the first 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, differences 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) find that the more extended Wilson (1975) models arebetter in describing the surface brightness profiles 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 different behaviour of the different components. Thesemodels are defined by a distribution function (DF) whichis a solution of the collisionless Boltzmann equation assum-ing a Maxwellian velocity distribution that is ‘lowered’ tomimic the effect of a negative escape energy as the resultof the galactic tides. The multimass formulation of a Kingmodel was first 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 & Griffin(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 defined in the Michie–King models (Michie1963) and multiple mass components as in Gunn & Griffin(1979). GZ15 introduced a power-law dependence betweenmass and anisotropy radius for each mass bin, while Gunn &Griffin (1979) argued that was not necessary because mostevents that influence 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 fit 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 defined 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 justified. 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 fitting the multimass mod-els to snapshots from different 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 fixed to δ = 1 /
2. By letting this parameter free, we can test whetherthis assumption is justified. By varying the amount of stellar-mass black holes (BHs) and neutron stars (NSs) retained inthe different N -body models, we also study their impact onthe cluster as well as on the different parameters of the best-fitting 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 differentphases 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 profiles of density, velocity dispersionand anisotropy are discussed in Section 5. In Section 6, wediscuss the values of the best-fitting 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 different 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 specific energy E = v / φ ( r ) is one of the two integrals of motion, where v isthe velocity and φ ( r ) the specific 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 defined 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 thespecific 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, defined 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 & Griffin (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 final parameter neededto define the models is the dimensionless central potential W (King 1966, ˆ φ in GZ15) which specifies 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 define 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 fits 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 five parameters ( g , W , δ , r a and η ) andtwo scales ( M Cl , r scale ) together with the description of themass bins ( M j , m j ) limepy first 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 find 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 & Griffin (1979) and GZ15, ¯ m is the central densityweighted mean mass. After performing several comparisons,we found that models calculated by the two different 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 modified compared to Da Costa & Freeman(1976); Gunn & Griffin (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 defini-tion ( W , ˆ r a ) to another ¯ m ∗ definition ( 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 first 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 differences 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 first and separately from the dynamical evolution.We do this because Galactic GCs have different dynamicalages but have all roughly the same physical age of around12 Gyr. We consider models with different 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 differ-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 fitting 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 effect 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 different 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 fit 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 first 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 thespecific 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 different as equation (11) neglectsthe tides. We adopted this definition nevertheless to be con-sistent with the multimass models, which also do not accountfor the changed potential due to tidal effects. 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 different 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 findings of Lee & Ostriker (1987) and Gieleset al. (2011) that the escape rate of clusters with the samemass mainly depends on the tidal field, which is the samefor all four models.In Fig. 2, we have plotted the evolution of r h for thefour different N -body models. As can be seen in the figure,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 effect 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 first 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 figure, 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 differentsnapshots, such as the dynamical age, the bound mass andthe number of NSs and BHs.
As in Peuten et al. (2016), we find that the mean massprofile is independent of the remnant retention fraction. InFig. 4, we plot this for all four N -body models at four differ-ent times in their evolution: 2 . τ rh , , 7 . τ rh , , 16 . τ rh , and26 τ rh , . Looking at the different 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 first snapshot at2 . τ rh , but over the course of evolution these differencesdiminish. The profiles get flatter 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 fourdifferent 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 first 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-fitting 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 Affine InvariantMCMC Ensemble sampler (Goodman & Weare 2010). The python implementation makes it straightforward to coupleit with limepy . Furthermore, we benefit from the fact thatwe created a distributed grid computing version of python ’s map() function, thereby dynamically distributing efficientlythe workload from several MCMC runs over all availableCPU cores at the University of Surrey Astrophysics com-puting facilities.The fitting 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 tofind those parameters which maximize the product of theseindividual likelihoods. The best-fitting 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 different researchtargets, leading to different 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 different star types, we opted againstmixing them and therefore we give every star type at leastone bin. This means that we have at least five 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 influences theresults of the analysis. For this, we choose four different 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 fixed 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 find that the results arealmost independent of the way one chooses the binning: thenumber of MS mass bins needed to converge is the sameand the difference between the different 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 different quantities con-verge. We find that we need at least two WDs mass bins andat least four MS stars mass bins for the different 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 five 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 effect 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 definitionand cannot be described by the model correctly. These ob-jects pose a problem because they drive the fit 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 artificialbackground population with a constant likelihood (i.e. uni-form distribution) in phase space. We added the artificialbackground 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 theartificial background and M Back is the mass of the back-ground only. The phase space volume V is defined 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 first 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 fitswith different initial conditions to test for any divergence.We chose flat priors restricted mainly by currently observedvalues for the parameters and/or by the range in which theyare considered physically valid. For the MCMC fitting, westarted out with around 500 walkers and found good fits forthe N -body model N0 without BHs and NSs. For the othermodels, prominently those with BHs, converging fits 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-fittingparameters. 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 affects 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 different fitting param-eters for two MCMC runs. Figs. 5 and 6 show the results ofthe fitting 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 figure shows the results of the MCMC fitting 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 different parameters can bealleviated. Tables A9–A12 list the best-fitting parameters for allthe N -body snapshots we considered. N -BODY MODELS In the first part of our analysis, we compare the best-fittingmultimass 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 figure shows the results of the MCMC fitting 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-fitting values of η and r a areunconstrained because this stellar system is not radially anisotropic. First, we compare the mass density profiles of the best-fitting 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 anddefine 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 profiles forthe three models N1, N0.3 and N0 at four different times:2 . τ rh , which is the first 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-fitting result from the multimass models, we MNRAS , 1–28 (2017) M. Peuten
Figure 7.
Comparison of the mass density profiles for models N1, N0.3 and N0 at four different ages: 2 . τ rh , , 7 . τ rh , , 16 . τ rh , and26 τ rh , . The points represent the binned N -body data, the thick lines represent the best-fitting multimass models profiles 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, reflecting the uncertaintiesof the results.As can be seen from Fig. 7, the best-fitting multimassmodels reproduce the mass density profiles of the differentmass components. Differences 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 difference 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 differences 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 profiles are slightlydifferent from isothermal models (Lynden-Bell & Eggleton1980). Given that the differences 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 profiles is consistent with the find-ings of Sollima et al. (2017) who fitted Michie–King modelsto observational data of NGC 5466, NGC 6218 and NGC6981: for all three GCs, the multimass models reproducethe observed mass density profiles (see their fig. 6). Differ-ences are only found in the outer regions, most likely for thesame reasons as discussed above.When one compares the different 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 effect 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 efficiently populated by the next lighter objects and theresulting mass density profile 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 profiles inFig. 8, we used the same snapshots as for the mass den-sity profiles. 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-fitting multimass models are in agreement with the datafrom the N -body models. As with the mass density profiles,small difference 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 undefined.When comparing the different models at the same dy-namical age we find that in clusters with BHs, the veloc-ity dispersions for the different 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 profile within thecluster and in Section 5.3.2 we consider the global anisotropyof the cluster as a whole.
The anisotropy parameter β is defined 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 profiles from thebest-fitting multimass models for a selection of mass binsto the anisotropy profiles from the N -body snapshots. Asthe β parameter is more affected by random scatter, we hadto bin the data from the N -body snapshots differently 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 findthat 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 differences between the best-fitting 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-fitting 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 profiles for models N1, N0.3 and N0 at four different ages: 2 . τ rh , , 7 . τ rh , , 16 . τ rh , and 26 τ rh , . The points represent the binned N -body data, the thick lines represent the best-fitting multimass models profiles 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 profiles for models N1, N0.3 and N0 at four different ages: 2 . τ rh , , 7 . τ rh , , 16 . τ rh , and26 τ rh , . The points represent the binned N -body data, the thick lines represent the best-fitting multimass models profiles 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 first isotropicand then tangential anisotropic. We will discuss the evolu-tion of the anisotropy further in Section 6.7 when we analysethe best-fitting η parameter. To quantify the global anisotropy we use the parameter κ ,that was introduced by Polyachenko & Shukhman (1981)and is defined 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-fitting multimass models and for the N -body snapshots.For the uncertainties of the N -body data, we used Pois-son statistics. The best-fitting 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-fittingvalue 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-fit ting parameters resulting fromour fitting 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 five 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-fitting 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 figure, thebest-fitting 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 difference is always smaller than 1% this effect 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 differencesbetween their study and ours which could lead to such dif-ferent results. They are using simulated observations as in-put for their analysis, which affects 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 effect 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-fitting multimass models divided by the true value as com-puted from the N -body data, for all four clusters throughouttheir lifetime. Again we find 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 fitting 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 different N -body modelsand/or at different 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 significant change ofthe W values in this figure.It is also possible to use the values of W obtainedfor different 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 different components for modelsN1, N0.3 and N0 at four different ages: 2 . τ rh , , 7 . τ rh , , 16 . τ rh , and 26 τ rh , . The red points represent the N -body data, the thickblack lines represent the best-fitting 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-fitting 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-fitting 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 effect 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-fitting 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-fitting 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 theeffect 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 fill their Roche volume, thetides interact with the clusters, stripping their outermoststars and thereby making the truncation in energy spacesteeper. This evolution is reflected 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-fitting 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 modelfindings 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 difference 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 final value. This be-haviour in the first 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-fitting 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) findings that clusters with large g are still dynamically expanding towards filling the Rochevolume (see also Carballo-Bello et al. 2012). In Fig. 17, we plot the evolution of the best-fitting 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 findings of Sollima et al.(2017), who study mass segregation in observations of GCs.At late stages, there are some snapshots for which the best-fitting 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-fitting multimass model to the observa-tions of NGC 6218 in Sollima et al. (2017). However, we donot find 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 different mass bins fromthe N -body data together with the predicted central velocitydispersion from the best-fitting multimass model (see alsoFig. 8). Given the results from Section 5.2, it is no surprisethat the best-fitting multimass models are able to reproducethe true values.Despite the fact that the best-fitting 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 different mass bins for Models N1, N0.3 and N0 at four differentdynamical ages: 2 . τ rh , , 7 . τ rh , , 16 . τ rh , and 26 τ rh , . The red points represent the binned N -body data, the black lines represent thebest-fitting 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 twodifferent 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 fixed value of δ = 0 . δ must therefore be smaller, some-thing also found by Sollima et al. (2015, 2017) for youngclusters. The last two fitting 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-fitting 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 definition 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 differently 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 differently: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-fitting 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-fitting 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 first have a look at the fitting pa-rameter η , which is needed to include a dependence of theanisotropy radius on the mass. The anisotropy parameter η is a novel fitting 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-fitting 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 reflected 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 efficiently, 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 fits to model N0 butthis time fixing η = 0 comparable to the original formulationby Gunn & Griffin (1979). In these fits, the most obviousdifference is that with increasing time, and therefore alsowith a higher absolute value of η , the uncertainties of r a increase up to five times the value recovered in the fits with anon-fixed η 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 fit 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 figure 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-fitting 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-fitting multimassmodels qualitatively reproduce the overall trend as long asthe N -body model is radially anisotropic. When the N -bodysnapshots becomes tangentially anisotropic, the best-fittingmultimass 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)fitting them to N -body models. We find 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 differences are negligible for most applicationsand parameters of interest (see also the discussion in Sec-tion 6.1).Our comparison shows that the best-fitting total clus-ter masses are off by no more than 1% from the true valueas computed from the N -body snapshots. The best-fittingcluster half-mass radius is reproduced within 5% and thetruncation radius is reproduced within 3%.We find that the mass density and velocity dispersionprofiles of the different 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 effects stripping the loosely bound stars.We find that the best-fitting 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 different N -body models over their wholelifetime. The black dashed lines represent the best-fitting 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-fitting 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 find 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 first τ rh , , does not show any strong differ-ences to simulation N0, despite having a population of NSs.The influence 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 fi-nancial support from the Royal Society (University ResearchFellowship) and AZ acknowledges financial 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, 14Chernoff 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., Griffin 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., Griffin 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 Pacific 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 Pacific 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 different 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: five 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 different 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: five 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 different 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: five 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 different 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: five 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 fitting 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 fitting 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-fitting 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 fitting 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 fitting 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-fitting 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 fitting 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 fitting 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-fitting 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 fitting 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 fitting 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-fitting 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