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

TL;DR

This study evaluates the effectiveness of multimass isothermal models in replicating the phase space properties of globular clusters as simulated by direct N-body methods, across various evolutionary stages and remnant retention scenarios.

Contribution

It demonstrates that multimass models accurately reproduce density and velocity profiles of star clusters in different evolutionary phases and provides insights into the evolution of global parameters and observable features.

Findings

Multimass models successfully match N-body simulation profiles.

Radial anisotropy develops from low- to high-mass components over time.

Velocity scale depends on mass as m^{-0.5}, with variations influenced by dark remnants.

Abstract

Lowered isothermal models, such as the multimass Michie-King models, have been successful in describing observational data of globular clusters. In this study we assess whether such models are able to describe the phase space properties of evolutionary $N$-body models. We compare the multimass models as implemented in (Gieles \& Zocchi) to $N$-body models of star clusters with different retention fractions for the black holes and neutron stars evolving in a tidal field. We find that multimass models successfully reproduce the density and velocity dispersion profiles of the different mass components in all evolutionary phases and for different remnants retention. We further use these results to study the evolution of global model parameters. We find that over the lifetime of clusters, radial anisotropy gradually evolves from the low-mass to the high-mass components and we identify features in the properties of observable stars that are indicative of the presence of stellar-mass black holes. We find that the model velocity scale depends on mass as $m^{-\delta}$, with $\delta\simeq0.5$ for almost all models, but the dependence of central velocity dispersion on $m$ can be shallower, depending on the dark remnant content, and agrees well with that of the $N$-body models. The reported model parameters, and correlations amongst them, can be used as theoretical priors when fitting these types of mass models to observational data.

Tables

Table 1. Table 1: Properties of the snapshots from the N 𝑁 N -body model N1 with 100% initial BH and NS retention. We list the age in units of the initial half-mass relaxation time, the total bound mass in M ⊙ subscript M direct-product \,{\rm M}_{\odot} , the half-mass radius in pc, the number of BHs in the cluster, the number of NSs in the cluster and the Jacobi radius r J subscript 𝑟 J r_{\rm J} in pc calculated as in Section  3.2 . For the bound mass and the number of BHs and NSs, we also give the percentage relative to the initial values in brackets.

Age	$M_{\mathrm{CL}}$	Half-mass radius	Number of BHs	Number of NSs	$r_J$
(${\tau }_{\mathrm{rh},0}$)	($M_{\odot }$)	(pc)			(pc)
$0$	$37353$ ($100\%$)	$2.06$	$121$ ($100\%$)	$632$ ($100\%$)	$30.12$
$2.4$	$30408$ ($81\%$)	$5.03$	$119$ ($56\%$)	$551$ ($87\%$)	$27.98$
$4.9$	$25737$ ($69\%$)	$6.03$	$85$ ($40\%$)	$518$ ($82\%$)	$26.51$
$7.3$	$22077$ ($59\%$)	$6.56$	$59$ ($28\%$)	$503$ ($80\%$)	$25.24$
$9.7$	$19109$ ($51\%$)	$6.58$	$55$ ($26\%$)	$490$ ($78\%$)	$24.08$
$12.1$	$15827$ ($42\%$)	$6.54$	$43$ ($20\%$)	$474$ ($75\%$)	$22.68$
$14.6$	$12618$ ($34\%$)	$6.27$	$35$ ($17\%$)	$461$ ($73\%$)	$21.10$
$17.0$	$9434$ ($25\%$)	$5.71$	$26$ ($12\%$)	$439$ ($69\%$)	$19.22$
$19.4$	$6602$ ($18\%$)	$4.98$	$18$ ($8.5\%$)	$416$ ($66\%$)	$17.11$
$21.9$	$4094$ ($11\%$)	$4.09$	$12$ ($5.7\%$)	$368$ ($58\%$)	$14.62$
$24.3$	$2046$ ($5.5\%$)	$2.92$	$6$ ($2.8\%$)	$301$ ($48\%$)	$11.68$
$26.7$	$497$ ($1.3\%$)	$1.58$	$0$ ($0.0\%$)	$149$ ($24\%$)	$7.52$

Table 2. Table 2: Properties of the snapshots from the N 𝑁 N -body model N0.3 with 33% initial BH and NS retention. We list the age in units of the initial half-mass relaxation time, the total bound mass in M ⊙ subscript M direct-product \,{\rm M}_{\odot} , the half-mass radius in pc, the number of BHs in the cluster, the number of NSs in the cluster and the Jacobi radius r J subscript 𝑟 J r_{\rm J} in pc calculated as in Section  3.2 . For the bound mass and the number of BHs and NSs, we also give the percentage relative to the initial values in brackets.

Age	$M_{\mathrm{CL}}$	Half-mass radius	Number of BHs	Number of NSs	$r_J$
(${\tau }_{\mathrm{rh},0}$)	($M_{\odot }$)	(pc)			(pc)
$0$	$34404$ ($100\%$)	$2.06$	$66$ ($100\%$)	$211$ ($100\%$)	$29.33$
$2.3$	$31131$ ($90\%$)	$3.07$	$22$ ($33\%$)	$199$ ($94\%$)	$28.14$
$4.7$	$29061$ ($84\%$)	$3.13$	$8$ ($12\%$)	$189$ ($90\%$)	$27.50$
$7.0$	$26824$ ($78\%$)	$3.12$	$1$ ($1.5\%$)	$179$ ($85\%$)	$26.78$
$9.4$	$24151$ ($70\%$)	$3.27$	$0$ ($0\%$)	$150$ ($71\%$)	$25.88$
$11.7$	$21002$ ($61\%$)	$3.57$	$0$ ($0\%$)	$119$ ($56\%$)	$24.72$
$14.1$	$17922$ ($52\%$)	$3.57$	$0$ ($0\%$)	$94$ ($45\%$)	$23.45$
$16.4$	$14861$ ($43\%$)	$3.61$	$0$ ($0\%$)	$75$ ($36\%$)	$22.06$
$18.8$	$11785$ ($34\%$)	$3.43$	$0$ ($0\%$)	$59$ ($28\%$)	$20.45$
$21.1$	$8808$ ($26\%$)	$3.28$	$0$ ($0\%$)	$46$ ($22\%$)	$18.59$
$23.5$	$6077$ ($18\%$)	$2.88$	$0$ ($0\%$)	$37$ ($18\%$)	$16.49$
$25.8$	$3651$ ($11\%$)	$2.61$	$0$ ($0\%$)	$26$ ($12\%$)	$13.99$
$28.1$	$1489$ ($4.3\%$)	$1.89$	$0$ ($0\%$)	$15$ ($7.1\%$)	$10.53$
$30.5$	$268$ ($0.8\%$)	$1.08$	$0$ ($0\%$)	$7$ ($3.3\%$)	$6.06$

Table 3. Table 3: Properties of the snapshots from the N 𝑁 N -body model N0.1 with 10% initial BH and NS retention. We list the age in units of the initial half-mass relaxation time, the total bound mass in M ⊙ subscript M direct-product \,{\rm M}_{\odot} , the half-mass radius in pc, the number of BHs in the cluster, the number of NSs in the cluster and the Jacobi radius r J subscript 𝑟 J r_{\rm J} in pc calculated as in Section  3.2 . For the bound mass and the number of BHs and NSs, we also give the percentage relative to the initial values in brackets.

Age	$M_{\mathrm{CL}}$	Half-mass radius	Number of BHs	Number of NSs	$r_J$
(${\tau }_{\mathrm{rh},0}$)	($M_{\odot }$)	(pc)			(pc)
$0$	$33476$ ($100\%$)	$2.04$	$22$ ($100\%$)	$56$ ($100\%$)	$29.06$
$2.3$	$31156$ ($93\%$)	$2.46$	$5$ ($23\%$)	$54$ ($96\%$)	$28.15$
$4.7$	$28719$ ($86\%$)	$2.76$	$0$ ($0\%$)	$50$ ($89\%$)	$27.41$
$7.0$	$25388$ ($76\%$)	$3.27$	$0$ ($0\%$)	$24$ ($43\%$)	$26.33$
$9.4$	$22112$ ($66\%$)	$3.52$	$0$ ($0\%$)	$19$ ($34\%$)	$25.16$
$11.7$	$18956$ ($57\%$)	$3.60$	$0$ ($0\%$)	$12$ ($21\%$)	$23.91$
$14.1$	$15809$ ($47\%$)	$3.65$	$0$ ($0\%$)	$7$ ($13\%$)	$22.53$
$16.4$	$12625$ ($38\%$)	$3.52$	$0$ ($0\%$)	$5$ ($8.9\%$)	$20.94$
$18.7$	$9772$ ($29\%$)	$3.42$	$0$ ($0\%$)	$3$ ($5.4\%$)	$19.26$
$21.1$	$6926$ ($21\%$)	$3.11$	$0$ ($0\%$)	$1$ ($1.8\%$)	$17.23$
$23.4$	$4373$ ($13\%$)	$2.75$	$0$ ($0\%$)	$0$ ($0\%$)	$14.85$
$25.8$	$2172$ ($6.5\%$)	$2.29$	$0$ ($0\%$)	$0$ ($0\%$)	$11.89$
$28.1$	$6780$ ($2.0\%$)	$1.48$	$0$ ($0\%$)	$0$ ($0\%$)	$8.19$

Table 4. Table 4: Properties of the snapshots from the N 𝑁 N -body model N0 with no initial BH and NS retention. We list the age in units of the initial half-mass relaxation time, the total bound mass in M ⊙ subscript M direct-product \,{\rm M}_{\odot} , the half-mass radius in pc and the Jacobi radius r J subscript 𝑟 J r_{\rm J} in pc calculated as in Section  3.2 . For the bound mass we also give the percentage relative to the initial value in brackets.

Age	$M_{\mathrm{Cl}}$	Half-mass radius	$r_J$
(${\tau }_{\mathrm{rh},0}$)	($M_{\odot }$)	(pc)	(pc)
$0$	$33042$ ($100\%$)	$2.04$	$28.93$
$2.3$	$30886$ ($93\%$)	$2.26$	$28.07$
$4.7$	$27539$ ($83\%$)	$2.93$	$27.03$
$7.0$	$24261$ ($73\%$)	$3.25$	$25.93$
$9.3$	$21223$ ($64\%$)	$3.46$	$24.80$
$11.7$	$18167$ ($55\%$)	$3.53$	$23.56$
$14.0$	$15053$ ($46\%$)	$3.56$	$22.15$
$16.3$	$12026$ ($36\%$)	$3.47$	$20.58$
$18.7$	$9183$ ($28\%$)	$3.18$	$18.84$
$21.0$	$6367$ ($19\%$)	$2.89$	$16.75$
$23.3$	$3955$ ($12\%$)	$2.55$	$14.33$
$25.7$	$1951$ ($5.9\%$)	$2.12$	$11.42$
$28.0$	$455$ ($1.4\%$)	$1.06$	$7.18$

Table 5. Table 5: Mass bins of the different snapshots of N 𝑁 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 subscript 𝑀 𝑗 M_{j} and the mean mass m j subscript 𝑚 𝑗 m_{j} in units of M ⊙ subscript M direct-product \,{\rm M}_{\odot} . 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
	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$
(${\tau }_{\mathrm{rh},0}$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)
$2.4$	$2498$	$0.13$	$3288$	$0.2$	$4792$	$0.32$	$6155$	$0.5$	$3630$	$0.72$	$185$	$0.82$	$3650$	$0.6$	$2346$	$0.78$	$1521$	$1.08$	$841$	$1.53$	$1504$	$12.6$
$4.9$	$2024$	$0.13$	$2702$	$0.2$	$4004$	$0.32$	$5287$	$0.5$	$3194$	$0.72$	$164$	$0.82$	$3152$	$0.6$	$2075$	$0.78$	$1391$	$1.08$	$789$	$1.52$	$955$	$11.2$
$7.3$	$1626$	$0.13$	$2182$	$0.2$	$3354$	$0.32$	$4572$	$0.5$	$2878$	$0.72$	$152$	$0.82$	$2775$	$0.6$	$1877$	$0.78$	$1299$	$1.09$	$766$	$1.52$	$596$	$10.1$
$9.7$	$1252$	$0.13$	$1722$	$0.2$	$2773$	$0.32$	$3943$	$0.5$	$2585$	$0.72$	$138$	$0.82$	$2475$	$0.6$	$1704$	$0.78$	$1230$	$1.09$	$748$	$1.53$	$538$	$9.79$
$12.1$	$892$	$0.13$	$1280$	$0.2$	$2149$	$0.32$	$3283$	$0.5$	$2235$	$0.72$	$118$	$0.82$	$2104$	$0.6$	$1507$	$0.78$	$1154$	$1.09$	$725$	$1.53$	$379$	$8.81$
$14.6$	$568$	$0.13$	$874$	$0.2$	$1543$	$0.32$	$2549$	$0.5$	$1894$	$0.72$	$104$	$0.83$	$1728$	$0.6$	$1317$	$0.79$	$1055$	$1.09$	$706$	$1.53$	$280$	$7.99$
$17.0$	$317$	$0.13$	$528$	$0.2$	$987$	$0.32$	$1818$	$0.51$	$1519$	$0.72$	$82$	$0.83$	$1299$	$0.6$	$1060$	$0.79$	$957$	$1.1$	$675$	$1.54$	$192$	$7.38$
$19.4$	$144$	$0.13$	$269$	$0.2$	$544$	$0.33$	$1172$	$0.51$	$1119$	$0.72$	$62$	$0.83$	$881$	$0.6$	$817$	$0.79$	$831$	$1.1$	$642$	$1.54$	$119$	$6.63$
$21.9$	$47.5$	$0.13$	$98$	$0.2$	$230$	$0.33$	$611$	$0.51$	$728$	$0.73$	$41.4$	$0.83$	$496$	$0.61$	$541$	$0.79$	$647$	$1.11$	$574$	$1.56$	$80$	$6.68$
$24.3$	$8.99$	$0.13$	$22.7$	$0.21$	$63$	$0.33$	$213$	$0.52$	$325$	$0.73$	$15.7$	$0.83$	$204$	$0.61$	$272$	$0.8$	$418$	$1.13$	$474$	$1.57$	$30.7$	$5.11$
$26.7$	$0.12$	$0.12$	$0.24$	$0.24$	$4.52$	$0.35$	$14.4$	$0.53$	$50$	$0.74$	$2.51$	$0.84$	$15.4$	$0.62$	$47.2$	$0.81$	$119$	$1.18$	$243$	$1.63$	$0$	$0$

Table 6. Table 6: Mass bins of the different snapshots of N 𝑁 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 subscript 𝑀 𝑗 M_{j} and the mean mass m j subscript 𝑚 𝑗 m_{j} in units of M ⊙ subscript M direct-product \,{\rm M}_{\odot} . 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
	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$
(${\tau }_{\mathrm{rh},0}$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)
$2.3$	$2746$	$0.13$	$3609$	$0.2$	$5268$	$0.32$	$6732$	$0.5$	$3935$	$0.72$	$200$	$0.81$	$3975$	$0.6$	$2516$	$0.78$	$1622$	$1.08$	$302$	$1.52$	$224$	$10.2$
$4.7$	$2465$	$0.13$	$3276$	$0.2$	$4845$	$0.32$	$6340$	$0.5$	$3783$	$0.72$	$196$	$0.82$	$3805$	$0.6$	$2423$	$0.78$	$1588$	$1.08$	$285$	$1.51$	$53$	$6.62$
$7.0$	$2113$	$0.13$	$2869$	$0.2$	$4365$	$0.32$	$5889$	$0.5$	$3637$	$0.72$	$184$	$0.82$	$3603$	$0.6$	$3603$	$0.78$	$1546$	$1.08$	$269$	$1.50$	$6.31$	$6.31$
$9.3$	$1731$	$0.13$	$2386$	$0.2$	$3783$	$0.32$	$5356$	$0.5$	$3457$	$0.72$	$179$	$0.82$	$3329$	$0.6$	$2235$	$0.78$	$1478$	$1.08$	$218$	$1.45$	$0$	$0$
$11.7$	$1337$	$0.13$	$1782$	$0.2$	$3118$	$0.32$	$4691$	$0.5$	$3205$	$0.72$	$164$	$0.82$	$3005$	$0.6$	$2087$	$0.78$	$1354$	$1.07$	$168$	$1.41$	$0$	$0$
$14.0$	$964$	$0.13$	$1406$	$0.2$	$2471$	$0.32$	$4003$	$0.5$	$2941$	$0.72$	$155$	$0.82$	$2675$	$0.6$	$1933$	$0.78$	$1244$	$1.06$	$130$	$1.38$	$0$	$0$
$16.3$	$635$	$0.13$	$990$	$0.2$	$1861$	$0.32$	$3308$	$0.5$	$2636$	$0.72$	$143$	$0.82$	$2274$	$0.61$	$1772$	$0.79$	$1140$	$1.05$	$102$	$1.36$	$0$	$0$
$18.7$	$385$	$0.13$	$622$	$0.2$	$1274$	$0.32$	$2595$	$0.51$	$2263$	$0.72$	$126$	$0.82$	$1859$	$0.61$	$1559$	$0.79$	$1024$	$1.05$	$80$	$1.35$	$0$	$0$
$21.0$	$189$	$0.13$	$341$	$0.2$	$779$	$0.33$	$1862$	$0.51$	$1847$	$0.73$	$104$	$0.82$	$1390$	$0.61$	$1329$	$0.79$	$904$	$1.04$	$61$	$1.34$	$0$	$0$
$23.3$	$74$	$0.13$	$144$	$0.2$	$404$	$0.33$	$1169$	$0.52$	$1394$	$0.73$	$83$	$0.82$	$949$	$0.61$	$1055$	$0.8$	$756$	$1.04$	$49.1$	$1.33$	$0$	$0$
$25.7$	$19.1$	$0.13$	$45.8$	$0.2$	$144$	$0.33$	$587$	$0.52$	$892$	$0.73$	$55$	$0.81$	$533$	$0.61$	$729$	$0.8$	$612$	$1.04$	$43.4$	$1.32$	$0$	$0$
$28.0$	$2.92$	$0.13$	$5.93$	$0.2$	$21.5$	$0.35$	$148$	$0.53$	$367$	$0.74$	$20.7$	$0.83$	$161$	$0.61$	$349$	$0.81$	$393$	$1.04$	$19.7$	$1.31$	$0$	$0$
$30.3$	$0.17$	$0.17$	$0.39$	$0.39$	$0$	$0$	$6.1$	$0.55$	$38.5$	$0.74$	$3.32$	$0.83$	$8.71$	$0.62$	$69$	$0.85$	$133$	$1.06$	$9.18$	$1.31$	$0$	$0$

Table 7. Table 7: Mass bins of the different snapshots of N 𝑁 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 subscript 𝑀 𝑗 M_{j} and the mean mass m j subscript 𝑚 𝑗 m_{j} in units of M ⊙ subscript M direct-product \,{\rm M}_{\odot} . 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
	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$
(${\tau }_{\mathrm{rh},0}$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)
$2.3$	$2779$	$0.13$	$3658$	$0.2$	$5338$	$0.32$	$6844$	$0.5$	$3976$	$0.72$	$204$	$0.81$	$4039$	$0.6$	$2538$	$0.78$	$1644$	$1.08$	$82$	$1.52$	$53$	$10.6$
$4.7$	$2394$	$0.13$	$3220$	$0.2$	$4800$	$0.32$	$6374$	$0.5$	$3801$	$0.72$	$200$	$0.82$	$3827$	$0.6$	$2451$	$0.78$	$1577$	$1.08$	$75$	$1.50$	$0$	$0$
$7.0$	$1954$	$0.13$	$2683$	$0.2$	$4103$	$0.32$	$5703$	$0.5$	$3560$	$0.72$	$186$	$0.81$	$3515$	$0.6$	$2317$	$0.78$	$1334$	$1.06$	$33$	$1.36$	$0$	$0$
$9.4$	$1516$	$0.13$	$2142$	$0.2$	$3405$	$0.32$	$5015$	$0.5$	$3290$	$0.72$	$180$	$0.82$	$3183$	$0.6$	$2171$	$0.78$	$1185$	$1.05$	$26$	$1.35$	$0$	$0$
$11.7$	$1131$	$0.13$	$1652$	$0.2$	$2757$	$0.32$	$4323$	$0.5$	$3011$	$0.72$	$168$	$0.82$	$2824$	$0.6$	$2009$	$0.78$	$1065$	$1.04$	$16.3$	$1.36$	$0$	$0$
$14.1$	$776$	$0.13$	$1207$	$0.2$	$2121$	$0.32$	$3606$	$0.5$	$2706$	$0.72$	$158$	$0.82$	$2455$	$0.6$	$1838$	$0.78$	$935$	$1.03$	$9.26$	$1.33$	$0$	$0$
$16.4$	$480$	$0.13$	$787$	$0.2$	$1504$	$0.32$	$2845$	$0.51$	$2345$	$0.72$	$147$	$0.82$	$2064$	$0.61$	$1635$	$0.79$	$811$	$1.02$	$6.51$	$1.3$	$0$	$0$
$18.7$	$262$	$0.13$	$456$	$0.2$	$971$	$0.32$	$2148$	$0.51$	$2025$	$0.73$	$130$	$0.82$	$1651$	$0.61$	$1414$	$0.79$	$711$	$1.01$	$3.91$	$1.3$	$0$	$0$
$21.1$	$107$	$0.13$	$227$	$0.2$	$528$	$0.33$	$1411$	$0.51$	$1600$	$0.73$	$113$	$0.82$	$1179$	$0.61$	$1169$	$0.79$	$592$	$1.01$	$1.3$	$1.3$	$0$	$0$
$23.4$	$33.4$	$0.13$	$81$	$0.2$	$227$	$0.33$	$745$	$0.52$	$1121$	$0.73$	$89$	$0.82$	$724$	$0.61$	$871$	$0.79$	$481$	$1.0$	$0$	$0$	$0$	$0$
$25.8$	$5.55$	$0.13$	$14.1$	$0.2$	$57$	$0.33$	$291$	$0.53$	$596$	$0.74$	$55$	$0.83$	$305$	$0.61$	$511$	$0.8$	$337$	$1.01$	$0$	$0$	$0$	$0$
$28.1$	$0.14$	$0.14$	$0.48$	$0.24$	$3.68$	$0.33$	$53$	$0.55$	$196$	$0.75$	$21.6$	$0.83$	$57$	$0.61$	$192$	$0.82$	$155$	$1.01$	$0$	$0$	$0$	$0$

Table 8. Table 8: Mass bins of the different snapshots of N 𝑁 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 subscript 𝑀 𝑗 M_{j} and the mean mass m j subscript 𝑚 𝑗 m_{j} in units of M ⊙ subscript M direct-product \,{\rm M}_{\odot} . 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
	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$	$M_j$	$m_j$
(${\tau }_{\mathrm{rh},0}$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)
$2.3$	$2743$	$0.13$	$3613$	$0.2$	$5300$	$0.32$	$6817$	$0.5$	$3987$	$0.72$	$210$	$0.82$	$4016$	$0.6$	$2541$	$0.78$	$1659$	$1.08$
$4.7$	$2285$	$0.13$	$3091$	$0.2$	$4614$	$0.32$	$6174$	$0.5$	$3742$	$0.72$	$197$	$0.81$	$3703$	$0.6$	$2378$	$0.78$	$1356$	$1.06$
$7.0$	$1847$	$0.13$	$2543$	$0.2$	$3944$	$0.32$	$5515$	$0.5$	$3489$	$0.72$	$185$	$0.81$	$3385$	$0.6$	$2218$	$0.78$	$1135$	$1.04$
$9.4$	$1450$	$0.13$	$2055$	$0.2$	$3292$	$0.32$	$4874$	$0.5$	$3236$	$0.72$	$174$	$0.81$	$3074$	$0.6$	$2081$	$0.78$	$987$	$1.03$
$11.7$	$1071$	$0.13$	$1589$	$0.2$	$2665$	$0.32$	$4210$	$0.5$	$2984$	$0.72$	$159$	$0.81$	$2708$	$0.6$	$1917$	$0.78$	$867$	$1.02$
$14.1$	$730$	$0.13$	$1147$	$0.2$	$2053$	$0.32$	$3495$	$0.51$	$2667$	$0.72$	$140$	$0.81$	$2329$	$0.6$	$1747$	$0.78$	$745$	$1.01$
$16.4$	$453$	$0.13$	$747$	$0.2$	$1466$	$0.32$	$2763$	$0.51$	$2336$	$0.72$	$123$	$0.82$	$1945$	$0.61$	$1557$	$0.78$	$636$	$1.0$
$18.8$	$236$	$0.13$	$432$	$0.2$	$935$	$0.32$	$2045$	$0.51$	$1961$	$0.73$	$107$	$0.82$	$1550$	$0.61$	$1362$	$0.79$	$555$	$0.99$
$21.1$	$98$	$0.13$	$196$	$0.2$	$506$	$0.33$	$1316$	$0.51$	$1518$	$0.73$	$87$	$0.82$	$1098$	$0.61$	$1086$	$0.79$	$463$	$0.99$
$23.5$	$26.7$	$0.13$	$68$	$0.2$	$213$	$0.33$	$707$	$0.52$	$1045$	$0.73$	$66$	$0.82$	$681$	$0.61$	$795$	$0.79$	$353$	$0.98$
$25.8$	$4.32$	$0.13$	$13$	$0.21$	$52$	$0.34$	$248$	$0.53$	$578$	$0.74$	$43$	$0.83$	$287$	$0.61$	$473$	$0.8$	$252$	$0.98$
$28.1$	$0.11$	$0.11$	$0.61$	$0.2$	$2.31$	$0.33$	$25.4$	$0.54$	$139$	$0.76$	$10.9$	$0.84$	$29$	$0.62$	$138$	$0.82$	$111$	$0.98$

Table 9. Table 9: Results from the MCMC fitting process for the snapshots from the N 𝑁 N -body model N1 with initial 100% BH and NS retention. We list here the age of each snapshot in units of the initial half-mass relaxation time, the total cluster mass in M ⊙ subscript M direct-product \,{\rm M}_{\odot} , the half-mass radius in pc, the dimensionless central concentration for the global mean mass, the central mean mass, and the ES stars, the truncation parameter g 𝑔 g , the mass segregation parameter δ 𝛿 \delta , the anisotropy radius for the global mean mass and the central mean mass in pc, the anisotropy parameter η 𝜂 \eta and the truncation radius r t subscript 𝑟 t r_{\rm t} in pc. All parameters except the dimensionless central concentration for the global mean mass, and the ES stars, the anisotropy radius for the global mean mass and the truncation radius r t subscript 𝑟 t r_{\rm 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 ​ σ 1 𝜎 1\sigma uncertainties.

Age	$M_{\mathrm{Cl}}$	$r_h$	$W_0$	$W_{0,\mathrm{CMM}}$	$W_{0,\mathrm{ES}}$	$g$	$\delta$	$r_a$	$r_{a,\mathrm{CMM}}$	$\eta$	$r_t$
(${\tau }_{\mathrm{rh},0}$)	($M_{\odot }$)	(pc)						(pc)	(pc)		(pc)
$2.4$	${30455}_{-14}^{+14}$	${5.11}_{-0.02}^{+0.02}$	${1.7}_{-0.1}^{+0.1}$	${16.9}_{-0.4}^{+0.4}$	${2.8}_{-0.1}^{+0.1}$	${1.43}_{-0.02}^{+0.02}$	${0.350}_{-0.008}^{+0.009}$	${4.3}_{-0.2}^{+0.15}$	${7.5}_{-0.6}^{+0.7}$	${0.17}_{-0.02}^{+0.02}$	${31.2}_{-0.2}^{+0.2}$
$4.9$	${25790}_{-12}^{+13}$	${6.13}_{-0.02}^{+0.02}$	${1.8}_{-0.1}^{+0.1}$	${20.4}_{-0.2}^{+0.2}$	${3.3}_{-0.2}^{+0.03}$	${1.25}_{-0.02}^{+0.02}$	${0.368}_{-0.008}^{+0.008}$	${8.4}_{-0.4}^{+0.41}$	${3.6}_{-0.4}^{+0.4}$	$-{0.26}_{-0.03}^{+0.04}$	${28.37}_{-0.09}^{+0.1}$
$7.3$	${22140}_{-11}^{+11}$	${6.59}_{-0.02}^{+0.02}$	${2.40}_{-0.07}^{+0.07}$	${23.6}_{-1.0}^{+0.7}$	${4.1}_{-0.1}^{+0.1}$	${1.05}_{-0.01}^{+0.01}$	${0.362}_{-0.005}^{+0.005}$	${20}_{-2.3}^{+2.6}$	${6.2}_{-0.9}^{+0.8}$	$-{0.39}_{-0.06}^{+0.06}$	${26.73}_{-0.06}^{+0.06}$
$9.7$	${19175}_{-9.9}^{+9.3}$	${6.61}_{-0.02}^{+0.01}$	${2.68}_{-0.08}^{+0.07}$	${18}_{-1.2}^{+1.2}$	${4.3}_{-0.10}^{+0.08}$	${0.91}_{-0.01}^{+0.01}$	${0.331}_{-0.007}^{+0.007}$	${61}_{-31}^{+\infty }$	${29}_{-17}^{+4770}$	$-{0.4}_{-3.6}^{+0.2}$	${25.56}_{-0.05}^{+0.05}$
$12.1$	${15878}_{-9.1}^{+8.5}$	${6.51}_{-0.02}^{+0.02}$	${2.74}_{-0.06}^{+0.06}$	${20}_{-1.5}^{+1.2}$	${4.4}_{-0.08}^{+0.08}$	${0.8}_{-0.01}^{+0.02}$	${0.360}_{-0.007}^{+0.006}$	${102}_{-43}^{+\infty }$	${39}_{-18}^{+208}$	$-{0.4}_{-0.3}^{+0.2}$	${23.96}_{-0.05}^{+0.09}$
$14.6$	${12659}_{-8.2}^{+7.7}$	${6.23}_{-0.02}^{+0.02}$	${2.85}_{-0.06}^{+0.06}$	${19}_{-1.2}^{+1.6}$	${4.4}_{-0.08}^{+0.08}$	${0.72}_{-0.02}^{+0.02}$	${0.377}_{-0.007}^{+0.007}$	–	${4416}_{-2895}^{+2814}$	$-4_{-4}^{+4}$	${22.39}_{-0.04}^{+0.06}$
$17.0$	${9464}_{-7.3}^{+6.8}$	${5.64}_{-0.02}^{+0.02}$	${3.09}_{-0.06}^{+0.07}$	${20}_{-1.2}^{+1.3}$	${4.6}_{-0.08}^{+0.08}$	${0.64}_{-0.02}^{+0.02}$	${0.390}_{-0.007}^{+0.007}$	–	${3740}_{-2428}^{+2365}$	$-4_{-4}^{+4}$	${20.27}_{-0.05}^{+0.07}$
$19.4$	${6623}_{-6.9}^{+6.2}$	${4.96}_{-0.03}^{+0.03}$	${3.20}_{-0.09}^{+0.08}$	${18}_{-1.7}^{+2.9}$	${4.4}_{-0.1}^{+0.1}$	${0.64}_{-0.02}^{+0.02}$	${0.42}_{-0.01}^{+0.01}$	–	${1017}_{-741}^{+2149}$	$-4_{-4}^{+4}$	${18.23}_{-0.05}^{+0.07}$
$21.9$	${4107}_{-5.7}^{+5.0}$	${4.03}_{-0.04}^{+0.04}$	${3.5}_{-0.13}^{+0.13}$	${17}_{-2.1}^{+3.0}$	${4.3}_{-0.1}^{+0.1}$	${0.71}_{-0.04}^{+0.04}$	${0.42}_{-0.02}^{+0.02}$	–	${3041}_{-1941}^{+1843}$	$-4_{-4}^{+4}$	${15.9}_{-0.2}^{+0.3}$
$24.3$	${2050}_{-5.0}^{+4.2}$	${2.92}_{-0.04}^{+0.04}$	${4.3}_{-0.23}^{+0.04}$	${23}_{-5.4}^{+4.0}$	${4.5}_{-0.2}^{+0.4}$	${0.7}_{-0.04}^{+0.04}$	${0.50}_{-0.03}^{+0.02}$	–	${874}_{-576}^{+600}$	$-4_{-4}^{+4}$	${12.9}_{-0.1}^{+0.1}$
$26.7$	${495}_{-3.9}^{+2.7}$	${1.7}_{-0.1}^{+0.08}$	$5_{-1.2}^{+1.1}$	${8.9}_{-0.7}^{+1.0}$	${2.7}_{-1.1}^{+1.2}$	${0.7}_{-0.1}^{+0.2}$	${0.9}_{-0.3}^{+0.4}$	${72}_{-54}^{+266}$	${21}_{-11}^{+12}$	$-4_{-4}^{+4}$	${8.7}_{-0.2}^{+0.4}$

Table 10. Table 10: Results from the MCMC fitting process for the snapshots from the N 𝑁 N -body model N0.3 with initial 33% BH and NS retention. We list here the age of each snapshot in units of the initial half-mass relaxation time, the total cluster mass in M ⊙ subscript M direct-product \,{\rm M}_{\odot} , the half-mass radius in pc, the dimensionless central concentration for the global mean mass, the central mean mass, and the ES stars, the truncation parameter g 𝑔 g , the mass segregation parameter δ 𝛿 \delta , the anisotropy radius for the global mean mass and the central mean mass in pc, the anisotropy parameter η 𝜂 \eta and the truncation radius r t subscript 𝑟 t r_{\rm t} in pc. All parameters except the dimensionless central concentration for the global mean mass, and the ES stars, the anisotropy radius for the global mean mass and the truncation radius r t subscript 𝑟 t r_{\rm 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 ​ σ 1 𝜎 1\sigma uncertainties.

Age	$M_{\mathrm{Cl}}$	$r_h$	$W_0$	$W_{0,\mathrm{CMM}}$	$W_{0,\mathrm{ES}}$	$g$	$\delta$	$r_a$	$r_{a,\mathrm{CMM}}$	$\eta$	$r_t$
(${\tau }_{\mathrm{rh},0}$)	($M_{\odot }$)	(pc)						(pc)	(pc)		(pc)
$2.3$	${31166}_{-14}^{+13}$	${3.14}_{-0.01}^{+0.01}$	${4.22}_{-0.06}^{+0.06}$	${47.8}_{-0.7}^{+0.2}$	${7.76}_{-0.08}^{+0.07}$	${1.52}_{-0.02}^{+0.02}$	${0.370}_{-0.003}^{+0.003}$	${5.9}_{-0.1}^{+0.2}$	${25.0}_{-2.8}^{+2.6}$	${0.44}_{-0.02}^{+0.03}$	${32.0}_{-0.4}^{+0.4}$
$4.7$	${29103}_{-13}^{+14}$	${3.18}_{-0.01}^{+0.01}$	${4.65}_{-0.03}^{+0.02}$	${13.26}_{-0.05}^{+0.05}$	${8.30}_{-0.03}^{+0.02}$	${1.22}_{-0.01}^{+0.01}$	${0.356}_{-0.001}^{+0.002}$	${7.8}_{-0.1}^{+0.1}$	${13.0}_{-0.3}^{+0.3}$	${0.35}_{-0.02}^{+0.02}$	${28.6}_{-0.2}^{+0.1}$
$7.0$	${26860}_{-13}^{+14}$	${3.13}_{-0.01}^{+0.01}$	${4.69}_{-0.04}^{+0.04}$	${11.23}_{-0.04}^{+0.04}$	${9.03}_{-0.04}^{+0.04}$	${1.11}_{-0.01}^{+0.01}$	${0.414}_{-0.003}^{+0.003}$	${6.7}_{-0.2}^{+0.2}$	${6.5}_{-0.3}^{+0.3}$	$-{0.03}_{-0.04}^{+0.03}$	${27.55}_{-0.1}^{+0.1}$
$9.3$	${24196}_{-12}^{+12}$	${3.38}_{-0.01}^{+0.01}$	${5.38}_{-0.08}^{+0.08}$	${13.5}_{-0.1}^{+0.1}$	${10.75}_{-0.07}^{+0.08}$	${1.102}_{-0.009}^{+0.009}$	${0.458}_{-0.006}^{+0.006}$	${6.3}_{-0.2}^{+0.2}$	${4.0}_{-0.2}^{+0.2}$	$-{0.45}_{-0.04}^{+0.05}$	${26.8}_{-0.1}^{+0.1}$
$11.7$	${21047}_{-11}^{+11}$	${3.64}_{-0.01}^{+0.01}$	${5.67}_{-0.09}^{+0.1}$	${13.8}_{-0.2}^{+0.2}$	${10.97}_{-0.09}^{+0.08}$	${1.03}_{-0.01}^{+0.01}$	${0.460}_{-0.009}^{+0.009}$	${7.5}_{-0.2}^{+0.3}$	${4.0}_{-0.2}^{+0.3}$	$-{0.67}_{-0.05}^{+0.05}$	${25.56}_{-0.09}^{+0.09}$
$14.0$	${17972}_{-10}^{+9}$	${3.68}_{-0.01}^{+0.009}$	${5.7}_{-0.1}^{+0.1}$	${13.9}_{-0.3}^{+0.2}$	${11.0}_{-0.1}^{+0.1}$	${0.97}_{-0.01}^{+0.01}$	${0.49}_{-0.01}^{+0.01}$	${7.9}_{-0.3}^{+0.3}$	${3.8}_{-0.3}^{+0.3}$	$-{0.81}_{-0.06}^{+0.06}$	${24.34}_{-0.08}^{+0.09}$
$16.3$	${14905}_{-9}^{+9}$	${3.60}_{-0.01}^{+0.008}$	${5.5}_{-0.09}^{+0.1}$	${12.0}_{-0.1}^{+0.1}$	${10.10}_{-0.07}^{+0.07}$	${0.89}_{-0.01}^{+0.01}$	${0.50}_{-0.01}^{+0.01}$	${9.6}_{-0.5}^{+0.6}$	${4.5}_{-0.2}^{+0.2}$	$-{0.97}_{-0.07}^{+0.08}$	${22.92}_{-0.07}^{+0.07}$
$18.7$	${11820}_{-9}^{+8}$	${3.47}_{-0.01}^{+0.01}$	${6.1}_{-0.1}^{+0.1}$	${11.9}_{-0.2}^{+0.2}$	${10.16}_{-0.09}^{+0.08}$	${0.85}_{-0.01}^{+0.01}$	${0.48}_{-0.01}^{+0.01}$	${10.4}_{-0.7}^{+0.9}$	${4.7}_{-0.3}^{+0.3}$	$-{1.1}_{-0.1}^{+0.1}$	${21.33}_{-0.07}^{+0.08}$
$21.0$	${8832}_{-8}^{+7}$	${3.32}_{-0.04}^{+0.02}$	${6.0}_{-0.1}^{+0.1}$	${11.6}_{-0.2}^{+0.2}$	${9.7}_{-0.1}^{+0.1}$	${0.82}_{-0.02}^{+0.03}$	${0.53}_{-0.02}^{+0.02}$	${12}_{-1}^{+2}$	${6.0}_{-0.5}^{+0.6}$	$-{1.2}_{-0.2}^{+0.2}$	${19.6}_{-0.1}^{+0.2}$
$23.3$	${6089}_{-8}^{+7}$	${2.94}_{-0.02}^{+0.02}$	${6.7}_{-0.2}^{+0.2}$	${11.2}_{-0.2}^{+0.2}$	${9.4}_{-0.1}^{+0.1}$	${0.8}_{-0.03}^{+0.03}$	${0.51}_{-0.02}^{+0.02}$	${14}_{-2}^{+4}$	${5.7}_{-0.7}^{+0.9}$	$-{1.8}_{-0.4}^{+0.4}$	${17.3}_{-0.2}^{+0.2}$
$25.7$	${3658}_{-6}^{+6}$	${2.57}_{-0.03}^{+0.02}$	${7.3}_{-0.2}^{+0.2}$	${10.9}_{-0.2}^{+0.2}$	${9.1}_{-0.2}^{+0.2}$	${0.75}_{-0.03}^{+0.03}$	${0.48}_{-0.03}^{+0.03}$	${14}_{-4}^{+7}$	${5.5}_{-0.8}^{+1.4}$	$-{2.3}_{-0.8}^{+0.6}$	${14.91}_{-0.1}^{+0.1}$
$28.0$	${1494}_{-4}^{+3}$	${1.92}_{-0.04}^{+0.04}$	${6.4}_{-0.3}^{+0.3}$	${9.4}_{-0.3}^{+0.4}$	${7.4}_{-0.2}^{+0.2}$	${0.76}_{-0.05}^{+0.05}$	${0.65}_{-0.07}^{+0.07}$	${2621}_{-1921}^{+8346}$	${1034}_{-689}^{+696}$	$-4{.}_{-4.1}^{+4.1}$	${11.6}_{-0.2}^{+0.2}$
$30.3$	${268}_{-3}^{+2}$	${1.10}_{-0.06}^{+0.06}$	${7.1}_{-0.6}^{+0.4}$	${8.4}_{-0.5}^{+0.7}$	${4.7}_{-1.1}^{+1.0}$	${0.6}_{-0.1}^{+0.2}$	${0.6}_{-0.2}^{+0.2}$	${1140}_{-782}^{+1444}$	${737}_{-497}^{+504}$	$-4{.}_{-4.1}^{+4.1}$	${6.5}_{-0.2}^{+0.5}$

Table 11. Table 11: Results from the MCMC fitting process for the snapshots from the N 𝑁 N -body model N0.1 with initial 10% BH and NS retention. We list here the age of each snapshot in units of the initial half-mass relaxation time, the total cluster mass in M ⊙ subscript M direct-product \,{\rm M}_{\odot} , the half-mass radius in pc, the dimensionless central concentration for the global mean mass, the central mean mass, and the ES stars, the truncation parameter g 𝑔 g , the mass segregation parameter δ 𝛿 \delta , the anisotropy radius for the global mean mass and the central mean mass in pc, the anisotropy parameter η 𝜂 \eta and the truncation radius r t subscript 𝑟 t r_{\rm t} in pc. All parameters except the dimensionless central concentration for the global mean mass, and the ES stars, the anisotropy radius for the global mean mass and the truncation radius r t subscript 𝑟 t r_{\rm 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 ​ σ 1 𝜎 1\sigma uncertainties.

Age	$M_{\mathrm{Cl}}$	$r_h$	$W_0$	$W_{0,\mathrm{CMM}}$	$W_{0,\mathrm{ES}}$	$g$	$\delta$	$r_a$	$r_{a,\mathrm{CMM}}$	$\eta$	$r_t$
(${\tau }_{\mathrm{rh},0}$)	($M_{\odot }$)	(pc)						(pc)	(pc)		(pc)
$2.3$	${31177}_{-15}^{+15}$	${2.51}_{-0.01}^{+0.01}$	${5.41}_{-0.03}^{+0.04}$	${16.2}_{-0.1}^{+0.2}$	${8.84}_{-0.03}^{+0.03}$	${1.47}_{-0.01}^{+0.01}$	${0.295}_{-0.002}^{+0.002}$	${11.7}_{-0.4}^{+0.3}$	${41}_{-1}^{+1}$	${0.67}_{-0.02}^{+0.02}$	${31.1}_{-0.2}^{+0.3}$
$4.7$	${28763}_{-14}^{+14}$	${2.82}_{-0.01}^{+0.01}$	${6.18}_{-0.06}^{+0.07}$	${12.9}_{-0.1}^{+0.1}$	${11.56}_{-0.05}^{+0.05}$	${1.30}_{-0.01}^{+0.01}$	${0.387}_{-0.006}^{+0.006}$	${4.33}_{-0.08}^{+0.07}$	${4.6}_{-0.2}^{+0.2}$	${0.05}_{-0.02}^{+0.02}$	${30.1}_{-0.2}^{+0.2}$
$7.0$	${25430}_{-12}^{+11}$	${3.36}_{-0.02}^{+0.01}$	${5.62}_{-0.08}^{+0.09}$	${12.7}_{-0.1}^{+0.1}$	${11.46}_{-0.07}^{+0.08}$	${1.12}_{-0.01}^{+0.01}$	${0.453}_{-0.007}^{+0.008}$	${5.4}_{-0.1}^{+0.1}$	${4.9}_{-0.2}^{+0.2}$	$-{0.11}_{-0.03}^{+0.03}$	${27.7}_{-0.1}^{+0.1}$
$9.4$	${22153}_{-12}^{+11}$	${3.56}_{-0.01}^{+0.01}$	${5.9}_{-0.1}^{+0.1}$	${12.9}_{-0.1}^{+0.1}$	${11.60}_{-0.08}^{+0.09}$	${1.04}_{-0.01}^{+0.01}$	${0.46}_{-0.01}^{+0.01}$	${6.6}_{-0.2}^{+0.2}$	${4.5}_{-0.2}^{+0.2}$	$-{0.45}_{-0.05}^{+0.05}$	${26.2}_{-0.1}^{+0.1}$
$11.7$	${18998}_{-10}^{+10}$	${3.67}_{-0.01}^{+0.01}$	${5.8}_{-0.1}^{+0.1}$	${12.2}_{-0.1}^{+0.1}$	${11.22}_{-0.07}^{+0.08}$	${0.94}_{-0.01}^{+0.01}$	${0.46}_{-0.01}^{+0.01}$	${8.0}_{-0.3}^{+0.3}$	${4.6}_{-0.2}^{+0.2}$	$-{0.68}_{-0.06}^{+0.05}$	${24.76}_{-0.08}^{+0.08}$
$14.1$	${15846}_{-10}^{+9}$	${3.69}_{-0.02}^{+0.02}$	${6.1}_{-0.1}^{+0.1}$	${11.7}_{-0.1}^{+0.1}$	${11.00}_{-0.07}^{+0.08}$	${0.85}_{-0.01}^{+0.01}$	${0.46}_{-0.01}^{+0.01}$	${8.7}_{-0.4}^{+0.4}$	${4.4}_{-0.2}^{+0.2}$	$-{0.97}_{-0.07}^{+0.08}$	${23.40}_{-0.1}^{+0.06}$
$16.4$	${12658}_{-8}^{+8}$	${3.58}_{-0.02}^{+0.02}$	${6.1}_{-0.2}^{+0.1}$	${11.8}_{-0.2}^{+0.2}$	${10.96}_{-0.09}^{+0.1}$	${0.84}_{-0.01}^{+0.01}$	${0.50}_{-0.02}^{+0.02}$	${9.3}_{-0.5}^{+0.6}$	${4.5}_{-0.3}^{+0.3}$	$-{1.1}_{-0.1}^{+0.1}$	${21.90}_{-0.08}^{+0.09}$
$18.7$	${9799}_{-8}^{+8}$	${3.45}_{-0.01}^{+0.02}$	${6.8}_{-0.2}^{+0.2}$	${12.1}_{-0.2}^{+0.2}$	${11.1}_{-0.1}^{+0.1}$	${0.80}_{-0.01}^{+0.01}$	${0.49}_{-0.02}^{+0.02}$	${10.1}_{-0.7}^{+0.9}$	${4.5}_{-0.3}^{+0.4}$	$-{1.4}_{-0.2}^{+0.1}$	${20.2}_{-0.1}^{+0.1}$
$21.1$	${6945}_{-8}^{+7}$	${3.06}_{-0.03}^{+0.03}$	${6.8}_{-0.2}^{+0.2}$	${11.1}_{-0.2}^{+0.2}$	${10.4}_{-0.1}^{+0.1}$	${0.84}_{-0.02}^{+0.02}$	${0.53}_{-0.02}^{+0.02}$	${14}_{-2}^{+3}$	${6.6}_{-0.7}^{+0.8}$	$-{1.7}_{-0.4}^{+0.3}$	${18.5}_{-0.1}^{+0.2}$
$23.4$	${4384}_{-7}^{+6}$	${2.77}_{-0.03}^{+0.03}$	${7.4}_{-0.1}^{+0.2}$	${10.8}_{-0.2}^{+0.2}$	${10.1}_{-0.1}^{+0.1}$	${0.80}_{-0.02}^{+0.03}$	${0.53}_{-0.03}^{+0.03}$	${1594}_{-1418}^{+8138}$	${504}_{-471}^{+491}$	$-4{.}_{-4}^{+4}$	${15.9}_{-0.2}^{+0.2}$
$25.8$	${2177}_{-5}^{+4}$	${2.29}_{-0.05}^{+0.05}$	${8.2}_{-0.2}^{+0.3}$	${10.7}_{-0.3}^{+0.3}$	${9.9}_{-0.2}^{+0.2}$	${0.79}_{-0.04}^{+0.05}$	${0.49}_{-0.05}^{+0.05}$	${1099}_{-772}^{+2822}$	${466}_{-304}^{+305}$	$-4{.}_{-4}^{+4}$	${13.2}_{-0.2}^{+0.2}$
$28.1$	${681}_{-3}^{+2}$	${1.42}_{-0.06}^{+0.04}$	${7.5}_{-0.3}^{+0.3}$	${8.6}_{-0.3}^{+0.3}$	${8.0}_{-0.2}^{+0.2}$	${0.68}_{-0.06}^{+0.09}$	${0.5}_{-0.1}^{+0.1}$	${1338}_{-874}^{+1719}$	${858}_{-584}^{+589}$	$-4{.}_{-4}^{+4}$	${9.2}_{-0.2}^{+0.3}$

Table 12. Table 12: Results from the MCMC fitting process for the snapshots from the N 𝑁 N -body model N0 with no initial BH and NS retention. We list here the age of each snapshot in units of the initial half-mass relaxation time, the total cluster mass in M ⊙ subscript M direct-product \,{\rm M}_{\odot} , the half-mass radius in pc, the dimensionless central concentration for the global mean mass, the central mean mass, and the ES stars, the truncation parameter g 𝑔 g , the mass segregation parameter δ 𝛿 \delta , the anisotropy radius for the global mean mass and the central mean mass in pc, the anisotropy parameter η 𝜂 \eta and the truncation radius r t subscript 𝑟 t r_{\rm t} in pc. All parameters except the dimensionless central concentration for the global mean mass, and the ES stars, the anisotropy radius for the global mean mass and the truncation radius r t subscript 𝑟 t r_{\rm 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 ​ σ 1 𝜎 1\sigma uncertainties.

Age	$M_{\mathrm{Cl}}$	$r_h$	$W_0$	$W_{0,\mathrm{CMM}}$	$W_{0,\mathrm{ES}}$	$g$	$\delta$	$r_a$	$r_{a,\mathrm{CMM}}$	$\eta$	$r_t$
(${\tau }_{\mathrm{rh},0}$)	($M_{\odot }$)	(pc)						(pc)	(pc)		(pc)
$2.3$	${30921}_{-15}^{+13}$	${2.30}_{-0.01}^{+0.01}$	${5.65}_{-0.07}^{+0.07}$	${11.9}_{-0.1}^{+0.1}$	${11.29}_{-0.05}^{+0.05}$	${1.57}_{-0.01}^{+0.01}$	${0.415}_{-0.007}^{+0.007}$	${4.7}_{-0.1}^{+0.1}$	${7.0}_{-0.3}^{+0.4}$	${0.45}_{-0.02}^{+0.02}$	${33.0}_{-0.2}^{+0.3}$
$4.7$	${27575}_{-13}^{+13}$	${3.01}_{-0.01}^{+0.01}$	${5.9}_{-0.1}^{+0.1}$	${13.1}_{-0.1}^{+0.1}$	${12.09}_{-0.08}^{+0.08}$	${1.28}_{-0.02}^{+0.01}$	${0.439}_{-0.008}^{+0.008}$	${4.60}_{-0.09}^{+0.08}$	${5.4}_{-0.3}^{+0.3}$	${0.17}_{-0.02}^{+0.02}$	${29.4}_{-0.3}^{+0.1}$
$7.0$	${24303}_{-13}^{+12}$	${3.32}_{-0.01}^{+0.01}$	${5.7}_{-0.1}^{+0.1}$	${12.8}_{-0.1}^{+0.1}$	${11.93}_{-0.07}^{+0.08}$	${1.16}_{-0.01}^{+0.01}$	${0.466}_{-0.009}^{+0.009}$	${5.8}_{-0.2}^{+0.2}$	${5.1}_{-0.2}^{+0.3}$	$-{0.14}_{-0.04}^{+0.04}$	${27.12}_{-0.07}^{+0.09}$
$9.4$	${21266}_{-11}^{+11}$	${3.55}_{-0.02}^{+0.01}$	${6.1}_{-0.1}^{+0.1}$	${13.1}_{-0.1}^{+0.1}$	${12.17}_{-0.09}^{+0.1}$	${1.08}_{-0.01}^{+0.01}$	${0.47}_{-0.01}^{+0.01}$	${7.4}_{-0.3}^{+0.3}$	${4.8}_{-0.2}^{+0.3}$	$-{0.51}_{-0.05}^{+0.05}$	${25.74}_{-0.1}^{+0.09}$
$11.7$	${18207}_{-10}^{+10}$	${3.59}_{-0.02}^{+0.01}$	${6.2}_{-0.1}^{+0.1}$	${12.7}_{-0.2}^{+0.2}$	${11.93}_{-0.1}^{+0.09}$	${0.99}_{-0.01}^{+0.01}$	${0.47}_{-0.01}^{+0.01}$	${8.1}_{-0.3}^{+0.4}$	${4.6}_{-0.3}^{+0.3}$	$-{0.75}_{-0.07}^{+0.06}$	${24.28}_{-0.08}^{+0.09}$
$14.1$	${15089}_{-9}^{+9}$	${3.65}_{-0.01}^{+0.01}$	${6.6}_{-0.2}^{+0.2}$	${13.3}_{-0.2}^{+0.2}$	${12.3}_{-0.1}^{+0.1}$	${0.94}_{-0.01}^{+0.01}$	${0.48}_{-0.01}^{+0.01}$	${9.0}_{-0.4}^{+0.5}$	${4.5}_{-0.3}^{+0.3}$	$-{0.95}_{-0.08}^{+0.08}$	${23.05}_{-0.1}^{+0.07}$
$16.4$	${12058}_{-9}^{+8}$	${3.50}_{-0.01}^{+0.01}$	${7.0}_{-0.1}^{+0.2}$	${12.3}_{-0.1}^{+0.1}$	${11.86}_{-0.1}^{+0.09}$	${0.86}_{-0.01}^{+0.01}$	${0.45}_{-0.02}^{+0.02}$	${9.7}_{-0.6}^{+0.8}$	${4.4}_{-0.2}^{+0.2}$	$-{1.3}_{-0.1}^{+0.1}$	${21.40}_{-0.06}^{+0.07}$
$18.8$	${9209}_{-8}^{+7}$	${3.22}_{-0.02}^{+0.02}$	${6.4}_{-0.1}^{+0.1}$	${11.0}_{-0.2}^{+0.2}$	${10.73}_{-0.09}^{+0.1}$	${0.81}_{-0.02}^{+0.02}$	${0.52}_{-0.02}^{+0.02}$	${11}_{-1}^{+1}$	${5.5}_{-0.4}^{+0.4}$	$-{1.4}_{-0.2}^{+0.2}$	${19.62}_{-0.07}^{+0.1}$
$21.1$	${6384}_{-7}^{+6}$	${2.95}_{-0.03}^{+0.03}$	${7.5}_{-0.2}^{+0.2}$	${11.7}_{-0.2}^{+0.2}$	${11.2}_{-0.1}^{+0.2}$	${0.86}_{-0.03}^{+0.02}$	${0.50}_{-0.02}^{+0.02}$	${12}_{-1}^{+2}$	${4.7}_{-0.4}^{+0.5}$	$-{2.0}_{-0.3}^{+0.3}$	${17.8}_{-0.1}^{+0.2}$
$23.5$	${3966}_{-6}^{+5}$	${2.54}_{-0.03}^{+0.03}$	${7.3}_{-0.2}^{+0.2}$	${10.6}_{-0.2}^{+0.2}$	${10.1}_{-0.1}^{+0.2}$	${0.79}_{-0.03}^{+0.03}$	${0.55}_{-0.04}^{+0.04}$	${16}_{-4}^{+11}$	$7_{-1}^{+2}$	$-{2.5}_{-1.0}^{+0.7}$	${15.2}_{-0.1}^{+0.2}$
$25.8$	${1955}_{-5}^{+4}$	${2.10}_{-0.04}^{+0.04}$	${7.4}_{-0.3}^{+0.3}$	${10.3}_{-0.3}^{+0.3}$	${9.5}_{-0.2}^{+0.1}$	${0.81}_{-0.04}^{+0.04}$	${0.66}_{-0.06}^{+0.06}$	${1369}_{-983}^{+3098}$	${592}_{-391}^{+390}$	$-4_{-4}^{+4}$	${12.6}_{-0.3}^{+0.3}$
$28.1$	${455}_{-3}^{+2}$	${1.01}_{-0.05}^{+0.04}$	${7.6}_{-0.4}^{+0.3}$	${8.7}_{-0.3}^{+0.4}$	${8.2}_{-0.3}^{+0.3}$	${1.0}_{-0.1}^{+0.1}$	${0.6}_{-0.2}^{+0.2}$	${19}_{-10}^{+17}$	${12}_{-6}^{+6}$	$-4_{-4}^{+4}$	${8.4}_{-0.5}^{+0.7}$

Equations

$$f_j\left(E,J^2\right)=A_j\exp \left(-\frac{J^2}{2r_{\mathrm{a},j}^2{s}_j^2}\right)E_{\gamma }\left(g,-\frac{E-\varphi \left(r_{\mathrm{t}}\right)}{s_j^2}\right),$$

$$E_{\gamma }\left(a,x\right)=\{\begin{array}{ll}\exp \left(x\right)& a=0\\ \exp \left(x\right)P\left(a,x\right)& a>0\end{array}$$

$$r_{\mathrm{a},j}=r_{\mathrm{a}}{\mu }_j^{\eta }$$

$${\mu }_j=\frac{m_j}{\bar{m}}$$

$$s_j=s{\mu }_j^{-\delta }$$

$${\rho }_j=\int f_j\left(E,J^2\right)d^{3}v$$

$${\nabla }^2\hat{\varphi }=-9\sum _j{\alpha }_j{\hat{\rho }}_j$$

$$W_0^{\ast }=W_0{\left(\frac{{\bar{m}}^{\ast }}{\bar{m}}\right)}^{2\delta }$$

$${\hat{r}}_{\mathrm{a}}^{\ast }={\hat{r}}_{\mathrm{a}}{\left(\frac{{\bar{m}}^{\ast }}{\bar{m}}\right)}^{(\eta +\delta )}$$

$$r_{\mathrm{J}}={\left(\frac{G{M}_{\mathrm{Cl}}}{2{\Omega }^2}\right)}^{1/3}$$

$$E_{\mathrm{crit}}=\varphi \left(r_{\mathrm{J}}\right)=-\frac{G{M}_{\mathrm{Cl}}}{r_{\mathrm{J}}}$$

$${\mathfrak{L}}_B=\frac{M_{\mathrm{Back}}}{V{M}_{\mathrm{tot}}}$$

$$V=\frac{4}{3}\pi {\left(2r_{\max }\right)}^3\times \frac{4}{3}\pi {\left(2v_{\max }\right)}^3$$

$$\mathfrak{L}=\frac{f(E,J^2)}{M_{\mathrm{tot}}}+\frac{M_{\mathrm{Back}}}{V{M}_{\mathrm{tot}}}$$

$${\sigma }_k^2=\frac{{\sum }_i^N{m}_i{\left[v_{k,i}-⟨v_k⟩\right]}^2}{{\sum }_i^N{m}_i}k=r,\theta ,\varphi$$

$$⟨v_k⟩=\frac{{\sum }_i^N{m}_i{v}_{k,i}}{{\sum }_i^N{m}_i}k=r,\theta ,\varphi$$

$$\beta \equiv 1-\frac{{\sigma }_{\mathrm{t}}^2}{2{\sigma }_{\mathrm{r}}^2}$$

$$\kappa =\frac{2K_{\mathrm{r}}}{K_{\mathrm{t}}}$$

$$m_j{s}_j^2=m_i{s}_i^2\forall j,ij\ne i$$

$$m_j{s}_j^2=m_j{\sigma }_{1d,j0}^2\forall j$$

Full text

††pagerange: Testing lowered isothermal models with direct $N$-body simulations of globular clusters - II. Multimass models–12††pubyear: 2017

1 Introduction

The amount of available data from observations of globular clusters (GCs) is steadily increasing. With the arrival of the ESA–Gaia data (Gaia Collaboration et al., 2016), we are entering the era of high-precision kinematics, allowing us to study properties of GCs with unprecedented detail. This calls for adequate methods of analysing and describing them in an equally detailed way. Despite the fact that GCs are thought to be free of dark matter (Baumgardt et al., 2010; Ibata et al., 2013; Baumgardt, 2017), and to have evolved to spherical and isotropic configurations as the result of two-body relaxation, GCs are complex systems to model. They consist of stars and stellar remnants with different masses and luminosities and primordial and dynamically processed binary stars (Heggie, 1975; Goodman & Hut, 1989; Hut et al., 1992; Heggie et al., 2006; Trenti et al., 2007). The mass and luminosity functions depend on the stellar initial mass function (IMF), age and metallicity. GC stellar populations display chemical anomalies (Gratton et al., 2004) and broadened main sequences (MS), possibly the result of variations in the helium abundance (Milone et al., 2014). Furthermore GCs evolve in a galactic tidal field that influences their evolution and present-day properties (Chernoff & Weinberg, 1990; Johnston et al., 1999; Takahashi & Portegies Zwart, 2000; Baumgardt & Makino, 2003; Küpper 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 cluster, Zonoozi et al. (2011) modelled a low-mass GC and finally Heggie (2014) and Wang et al. (2016) presented the first $N$-body simulations of GCs with $N\sim 10^{6}$. The faster Monte Carlo method allows to explore the parameter of the initial conditions to some extent (Heggie & Giersz, 2008; Giersz et al., 2013). To infer properties for a large number of GCs with models with several degrees of freedom, static models that are fast to calculate are required. By using relatively simple models, that are motivated by the underlying physical processes that drive their evolution, differences between models and observations can be used to increase our understanding (Binney & McMillan, 2011).

In the context of GCs, the King (1966) models are often compared to observations, although they cannot describe all GCs successfully. For example, McLaughlin & van der Marel (2005) find that the more extended Wilson (1975) models are better in describing the surface brightness profiles of some Galactic GCs. In addition, both King and Wilson models have isothermal cores, which are not able to describe the late stages of core collapse (Lynden-Bell & Eggleton, 1980; Cohn, 1980). The models we are going to test and discuss in the context of GCs are multimass, anisotropic and spherical models (hereafter multimass models), which describe the properties of GCs considering their stellar mass function (MF) in the form of mass bins. This formulation allows for different behaviour of the different components. These models are defined by a distribution function (DF) which is a solution of the collisionless Boltzmann equation assuming a Maxwellian velocity distribution that is ‘lowered’ to mimic the effect of a negative escape energy as the result of the galactic tides. The multimass formulation of a King model was first introduced by Da Costa & Freeman (1976, we note that a formulation of a multimass model was already presented in ). Gunn & Griffin (1979) extended these models including radial anisotropy as formulated by Eddington (1915) for isothermal models (see also Michie, 1963).

Since their introduction, multimass models have been successfully used in a multitude of studies such as in Illingworth & 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 & Lee, 2000; Sollima et al., 2015).

Davoust (1977) realized that the DFs of the Woolley (1954), King (1966) and Wilson (1975) models can be written as a single DF with one additional integer parameter. Gomez-Leyton & Velazquez (2014) further generalized this formulation, allowing to calculate models in between the three classical models. (Gieles & Zocchi, 2015, hereafter GZ15) took up these formulation and added radial anisotropy as defined in the Michie–King models (Michie, 1963) and multiple mass components as in Gunn & Griffin (1979). GZ15 introduced a power-law dependence between mass and anisotropy radius for each mass bin, while Gunn & Griffin (1979) argued that was not necessary because most events that influence the anisotropy are mass independent or not very important. GZ15 also implemented the possibility to change the degree of mass segregation with an additional parameter $\delta$ that describes the relation between the velocity scale and mass, which in most models is assumed to be equal to $1/2$.

Despite their success in describing observational data, multimass models have been criticized for several assumptions made in their construction (see McLaughlin 2003 and Meylan & Heggie 1997). One such aspect is that multimass models have more parameters than single-mass models: in the formulation by GZ15, there are $2N_{\rm MB}+5$ parameters and $2$ scales, with $N_{\rm MB}$ being the number of mass bins, compared to $3$ parameters and $2$ scales for the single-mass model. It is therefore easier to fit multimass models to the data because they have more degrees of freedom (McLaughlin, 2003). Not only the selection of the right number of mass bins, but also how they are defined is criticized as a ‘usual compromise between convenience and realism’, as Meylan & Heggie (1997) put it. Given the numerous studies successfully using multimass models, this problem does not seem to be too much of a concern, but we nevertheless explored it in our study. Another assumption for which multimass models have been criticized is the assumption of equipartition of energy (McLaughlin, 2003; Trenti & van der Marel, 2013). Indeed the velocity scale is usually assumed to scale with the mass as $m^{-1/2}$, but we note that evolutionary multimass models only achieve partial equipartition (Merritt 1981; Miocchi 2006; GZ15; Bianchini et al. 2016) as the result of the escape velocity.

Several aspects of multimass models were already analysed with the help of Fokker–Planck (Takahashi & Lee, 2000) and $N$-body simulations (Sollima et al., 2015). The goal of this study is to compare the multimass models in the formulation by GZ15 to a set of $N$-body models to assess the quality of the former and to analyse whether some of the above mentioned criticism is justified. In this comparison we do not include any source of uncertainties, such as observational biases, to see how good the models are under ideal conditions. Hence, we determine the MF of the multimass models directly from the $N$-body data.

The comparison is done by fitting the multimass models to snapshots from different $N$-body models using a Markov Chain Monte Carlo (MCMC) method. Additionally, we study the new parameters which are now available in the extended formulation of the models by GZ15. In particular, the continuous truncation parameter as introduced by Gomez-Leyton & Velazquez (2014) and the parameter that controls the mass dependence of the anisotropy for each mass bin. Furthermore, we study the behaviour of the mass segregation parameter $\delta$, which in previous studies was fixed to $\delta=1/2$. By letting this parameter free, we can test whether this assumption is justified. By varying the amount of stellar-mass black holes (BHs) and neutron stars (NSs) retained in the different $N$-body models, we also study their impact on the cluster as well as on the different parameters of the best-fitting models.

In Zocchi et al. (2016), a similar analysis was presented for single-mass models. This comparison showed that the single-mass models are successful in describing the different phases of the dynamical evolution. Zocchi et al. (2016) studied the development of radial anisotropy in GCs and found that the models can be used to put limits on the expected amount of radial anisotropy.

This paper is structured as follows: in the next section, we give a brief overview of the multimass models. Then in Section 3, we discuss how the $N$-body models were generated and we discuss their properties. In Section 4, we present the method used for the analysis and the challenges we encountered. The radial profiles of density, velocity dispersion and anisotropy are discussed in Section 5. In Section 6, we discuss the values of the best-fitting model parameters and scales, and their implications. Finally, in Section 7, we discuss our results and present our conclusions.

2 The limepy Models

The multimass models used in this study are provided by the limepy (Lowered Isothermal Model Explorer in PYthon)111https://github.com/mgieles/limepy software package (GZ15). A model has different components, each representing stars in a mass range, characterized by a mean and total mass. The DF of the $j$th mass component is given by

[TABLE]

for $E<\phi(r_{\rm t})$ and [math] otherwise. The specific energy $E={v^{2}}/{2}+\phi(r)$ is one of the two integrals of motion, where $v$ is the velocity and $\phi(r)$ the specific potential at a distance $r$ from the centre. The energy $E$ is lowered by the potential at the truncation radius $\phi(r_{\rm t})$. The function $E_{\gamma}$ is defined as

[TABLE]

with $P\left(a,x\right)\equiv\gamma\left(a,x\right)/\Gamma\left(a\right)$ the regularized lower incomplete gamma function. The other integral of motion is the specific angular momentum $J=rv_{\rm t}$, where $v_{\rm t}$ is the tangential component of the velocity vector.

The anisotropy radius $r_{\rm a}$ is a parameter that controls how anisotropic the model is. The system is isotropic in the centre, radially anisotropic in the intermediate part and near $r_{\rm t}$ it is isotropic again. For small values of $r_{\rm a}$, the models are strongly anisotropic and for values of $r_{\rm a}$ lager than $r_{\rm t}$, the models are completely isotropic. GZ15 include a power-law dependence between mass and anisotropy radius:

[TABLE]

with $\mu_{j}$ the dimensionless mean mass of stars in the $j$th mass component, defined as:

[TABLE]

where $m_{j}$ is the mean mass of stars in the $j$th component and $\bar{m}$ a reference mass which we set equal to the global mean mass. If $\eta$ is set to zero the anisotropy radius is independent of the mass as in Gunn & Griffin (1979). limepy expects $r_{\rm a}$ to be input in units of the King radius ($r_{\rm 0}$), $\hat{r}_{\rm a}=r_{\rm a}/r_{\rm 0}$, hence $\hat{r}_{\rm a}$ is the parameter we vary.

The truncation parameter $g$ was introduced by Gomez-Leyton & Velazquez (2014) and describes the polytropic part near the escape energy. The polytropic index $n$ relates to $g$ as $n=g+3/2$ and this formulation allows to calculate models in-between the classical models: for $g=0$ and $r_{\rm a}\gg r_{\rm t}$, the DF is identical to the one from the Woolley model (Woolley, 1954). A Michie–King (Michie, 1963; King, 1966) model is reproduced for $g=1$ and for $g=2$ one gets the non-rotating Wilson model (Wilson, 1975). The range of possible values for the model parameter are $0\leq g\leq 3.5$ because as discussed in Gomez-Leyton & Velazquez (2014) and GZ15, there are no finite models above $g=3.5$. The final parameter needed to define the models is the dimensionless central potential $W_{0}$ (King, 1966, $\hat{\phi}_{0}$ in GZ15) which specifies how centrally concentrated the model is. It is a boundary condition for solving Poisson’s equation.

Besides these parameters, there are also two constants which define the physical scales of the model: one is the global velocity scale $s$ and the other is the normalization constant $A$ which sets the phase space density. Instead of these scales, the code needs as input the total cluster mass $M_{\rm Cl}$ and a radial scale $r_{\rm scale}$ (which can be $r_{\rm 0}$, the half-mass radius $r_{\rm h}$, the viral radius $r_{\rm v}$ or $r_{\rm t}$), which are internally converted to $A$ and $s$.

The velocity scale $s_{j}$ is deduced from $s$ as:

[TABLE]

It is usually assumed that $\delta=1/2$, but in this study we determine 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 in each component ($M_{j}$), which the user provides together with $m_{j}$. It must be noted that $M_{\rm Cl}$ is a required input parameter, independent from the $M_{j}$ parameters, because the latter are only used to compute the relative masses in each component. Only after the model is solved, $\sum_{j}M_{j}=M_{\rm Cl}$.

Given these five parameters ($g$, $W_{0}$, $\delta$, $r_{\rm a}$ and $\eta$) and two scales ($M_{\rm Cl}$, $r_{\rm scale}$) together with the description of the mass bins ($M_{j},m_{j}$) limepy first calculates the density for each mass bin via:

[TABLE]

Then, the dimensionless Poisson equation is solved

[TABLE]

with $\alpha_{j}=\rho_{j,0}/\rho_{0}$, $\hat{\rho}_{j}=\rho_{j}/\rho_{j,0}$ and the dimensionless positive potential $\hat{\phi}=(\phi(r_{\rm t})-\phi)/s^{2}$, is iteratively solved by varying $\alpha_{j}$ until the calculated $M_{j}$ converges to the input values. After the model is solved, it is scaled to $M_{\rm Cl}$ and $r_{\rm scale}$. We can then find the likelihood for any phase space coordinate using the DF (equation 1).

In equation (4), we set $\bar{m}$ to the mean mass of the cluster. In the formulation by Da Costa & Freeman (1976); Gunn & Griffin (1979) and GZ15, $\bar{m}$ is the central density weighted mean mass. After performing several comparisons, we found that models calculated by the two different formulations give the same results within the numerical uncertainties. Furthermore we found that using the global mean mass instead speeds up the calculation, especially for models with BHs. When using the global $\bar{m}$ the meaning of two model parameters is modified compared to Da Costa & Freeman (1976); Gunn & Griffin (1979) and GZ15: $W_{0}$ and $r_{\rm a}$ both represent their value for a hypothetical mass group with a mass of $\bar{m}$. Besides computational improvement, this change from the original formulation also allows us to compare the multimass $W_{0}$ value with the single-mass $W_{0}$ value, as both represent the $W_{0}$ value for the mean mass group.

One can easily translate the values given in one $\bar{m}$ definition ($W_{0}$, $\hat{r}_{\rm a}$) to another $\bar{m}^{*}$ definition ($W_{0}^{*}$, $\hat{r}_{{\rm a}}^{*}$) by applying the following two equations:

[TABLE]

[TABLE]

The $\delta$ term in equation (9) comes from the $r_{0}$ dependence of $\hat{r}_{\rm a}$.

As further improvement to the original formulation of the limepy models we found that radially anisotropic models can be constructed faster if one first calculates the $M_{j}$ array of the corresponding isotropic model and then uses this model as starting point to solve the anisotropic model. As with the previous improvement the differences are only of numerical nature. This procedure is now implemented in the current distribution of limepy.

3 Description of the $N$-body models

For the computation of the $N$-body data, we use the approach presented in Trenti et al. (2010): the stellar evolution is done first and separately from the dynamical evolution. We do this because Galactic GCs have different dynamical ages but have all roughly the same physical age of around $12{\,\rm Gyr}$. We consider models with different retention fractions of NSs and BHs and analyse them at various dynamical ages. Temporal units are always expressed in units of the initial half-mass relaxation time ($\tau_{\rm{rh,0}}$) of the $N$-body model.

3.1 Set-up of the $N$-body models

For this analysis, we run four $N$-body models, with different amounts of NSs and BHs. Each $N$-body model was set up as a cluster with $N=10^{5}$ stars initially following the Hénon isochrone model (Hénon, 1959) with $r_{\rm h}=2.25{\,\rm pc}$. As IMF we adopted a Kroupa IMF (Kroupa, 2001) in the mass range between $0.1\,{\rm M}_{\odot}$ and $100\,{\rm M}_{\odot}$ without any primordial binaries. Then, by using the fitting formula by Hurley et al. (2000) and assuming a metallicity of ${\rm[Fe/H]}=-2$, the stars were evolved to an age of about $12{\,\rm Gyr}$.

We mimic the effect of supernova kick velocity by removing a certain fraction of NSs and BHs from the initial conditions described above. The retention fraction of NSs and BHs after supernova kicks is highly uncertain (Repetto et al., 2012; Mandel, 2016). To bracket all possible cases, we consider four different values for the fraction of remnants that we retain in the cluster: $100\%$ (all the remnants are retained, simulation N1), $33\%$ (simulation N0.3), $10\%$ (simulation N0.1) and $0\%$ (all the remnants are removed, simulation N0). The initial half-mass relaxation time for all four clusters was $\tau_{\rm{rh,0}}=350\,{\rm Myr}$ before the stellar evolution and the removal of the dark remnants, after these steps the $\tau_{\rm{rh,0}}$ values are $412\,{\rm Myr}$ for N1, $426\,{\rm Myr}$ for N0.3, $427\,{\rm Myr}$ for N0.1 and $428\,{\rm Myr}$ for N0.

The clusters are evolved on a circular orbit with a circular velocity of $V_{\rm c}=220\,\rm km/s$ at a distance of $R_{\rm G}=4{\,\rm kpc}$, in a singular isothermal galactic potential to mimic a galaxy. The equation of motion is solved in an inertial reference frame centred on the cluster.

These four stellar systems were then dynamically evolved 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 object reaching a distance greater than twice the Jacobi radius ($r_{\rm 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 the simulations were also only evolved dynamically.

A snapshot of each cluster is taken every ${\,\rm Gyr}$, resulting in 48 snapshots222The snapshots can be retrieved from http://astrowiki.ph.surrey.ac.uk/dokuwiki/doku.php?id=tests:collision:mock_data:challenge_2 (11 for model N1, 13 for model N0.3, 12 for model N0.1 and 12 for model N0) which we fit the multimass models to (Section 4).

3.2 Selecting bound objects

Because multimass models describe bound objects in a cluster, we removed any unbound object from the $N$-body models. We discuss here how we selected the unbound objects for each $N$-body snapshot.

First, we determine the Jacobi radius

[TABLE]

in which $M_{\rm Cl}$ is the total mass within $r_{\rm J}$ and $\Omega={V_{\rm c}}/{R_{\rm G}}$ is the orbital angular velocity. As a first guess, we set $M_{\rm Cl}$ equal to the total mass of all stars in the snapshot and then determine $r_{\rm J}$ through an iterative approach.

With $r_{\rm J}$ determined we are now able to calculate the specific critical energy which is equal to the potential at $r_{\rm J}$

[TABLE]

The true critical energy is different as equation (11) neglects the tides. We adopted this definition nevertheless to be consistent with the multimass models, which also do not account for the changed potential due to tidal effects. We considered an object bound if it is within $r_{\rm J}$ and for its energy it holds $E_{i}<E_{\rm crit}$, and we only used these bound objects in the rest of this analysis.

3.3 Properties of the $N$-body models

Fig. 1 shows how $M_{\rm Cl}$ for the four different $N$-body models evolves over the course of the simulation. It is apparent that the cluster with 100% initial BH and NS retention (simulation N1) has the highest initial mass loss, but there seems to be no direct correlation between the number of BHs and NSs and initial mass loss as can be seen from the other three $N$-body models. Over the course of evolution the four models seem to have aligned their mass-loss rate which is in accordance with the findings of Lee & Ostriker (1987) and Gieles et al. (2011) that the escape rate of clusters with the same mass mainly depends on the tidal field, which is the same for all four models.

In Fig. 2, we have plotted the evolution of $r_{\rm h}$ for the four different $N$-body models. As can be seen in the figure, increasing the retention fraction of the BHs leads to an expansion of the cluster: simulation N1, with $100\%$ NS and BH retention has an $r_{\rm h}$ which is on average twice as large as the $r_{\rm h}$ from simulation N0 with no NS and BH retention. The cluster in simulation N0.3 loses all of its BHs at around $7\tau_{\rm{rh,0}}$ and for the rest of the simulation its $r_{\rm h}$ resembles that of N0. The effect of stellar-mass BHs on the radius evolution was also described by Mackey et al. (2008). The global evolution of $r_{\rm h}$ however is essentially the same, independent of BHs and NSs retained, and follows the description in Gieles et al. (2011). In the first half of their lifetime, the clusters are in the expansion-dominated phase while in the second half the clusters are in the evaporation-dominated phase during which $r_{\rm h}$ decreases again until total dissolution.

Fig. 3 shows the relative number evolution for the BHs and NSs in the simulation N1, N0.3 and N0.1. As can be seen in the bottom figure, the cluster from $N$-body simulation N1 with $100\%$ initial BHs and NSs retention is the only cluster that retains its BHs almost until to the end of its lifetime, while the cluster from $N$-body simulation N0.3 loses all its BHs at around $7\tau_{\rm{rh,0}}$ and the one from simulation N0.1 already at $2\tau_{\rm{rh,0}}$. Looking at the NSs in the top of Fig. 3, we see their initial loss is not as strong as for the BHs. But as soon as all BHs have left the cluster, the NSs escape rate increases such that the cluster N0.1 loses all its NS at around $\sim 21\tau_{\rm{rh,0}}$. Only the cluster from $N$-body simulation N1 has a population of NS left at the end of its lifetime. After all BHs are lost, NSs, then being the most massive objects in the cluster, segregate to the centre and are then ejected from the cluster due to interactions they experience with each other.

Tables 1–4 list various properties of the different snapshots, such as the dynamical age, the bound mass and the number of NSs and BHs.

3.4 Mean mass at different radii

As in Peuten et al. (2016), we find that the mean mass profile is independent of the remnant retention fraction. In Fig. 4, we plot this for all four $N$-body models at four different times in their evolution: $2.3\tau_{\rm{rh,0}}$, $7.1\tau_{\rm{rh,0}}$, $16.5\tau_{\rm{rh,0}}$ and $26\tau_{\rm{rh,0}}$. Looking at the different times we see that the overall behaviour is the same for all models independent of their dark remnant population. Some divergences between the different $N$-body models can be seen in the first snapshot at $2.3\tau_{\rm{rh,0}}$ but over the course of evolution these differences diminish. The profiles get flatter over time. This is comparable to the behaviour found for a set of $N$-body models where the dynamical and stellar evolution were done concurrently in Peuten et al. (2016). Here, the evolution over time is less strong because the stellar evolution was done before the dynamical evolution. We are not aware of a theory providing an explanation for this attractor solution of $\bar{m}(r)$, but for single-mass system it is known that after several relaxation times the evolution becomes self-similar (Hénon, 1961, 1965). Also it had been shown that the evolution of radii and mass of the multimass systems is comparable to those of single-mass systems (Lee & Goodman, 1995; Gieles et al., 2010) but faster. Furthermore Giersz & Heggie (1996, 1997) showed in multimass $N$-body models that after some time the mean mass in Lagrangian shells stops evolving and, to a first approximation, stays constant. Although we do not explore this here, this result could be used as a theoretical prior when comparing multimass models to data.

4 Method

To determine the best-fitting multimass models for each snapshot we use the MCMC software package emcee (Foreman-Mackey et al., 2013), which is a pure-python implementation of the Goodman & Weare’s Affine Invariant MCMC Ensemble sampler (Goodman & Weare, 2010). The python implementation makes it straightforward to couple it with limepy. Furthermore, we benefit from the fact that we created a distributed grid computing version of python’s map() function, thereby dynamically distributing efficiently the workload from several MCMC runs over all available CPU cores at the University of Surrey Astrophysics computing facilities.

The fitting process consists of computing a multimass model based on the input parameters provided by the MCMC walker position in parameter space and the current mass bin description (see the next section). Then the likelihood for each star in all mass bins is calculated using the DF (see equation 14) and the phase space position of the star from the $N$-body snapshot. By randomly varying the walker positions in parameter space, the MCMC algorithm tries to find those parameters which maximize the product of these individual likelihoods. The best-fitting value of each parameter is estimated as the median of the marginalized posterior distribution using all walker positions from all chains after removing the initial burn-in phase. This generally coincides with the value of the parameter providing the largest likelihood. For the $1\sigma$ errors, we use the values from the 16th and 84th percentiles.

4.1 Determining the mass bins

As mentioned in Section 1, the mass bin selection for multimass models is in most cases a choice of convenience (Meylan & Heggie, 1997) as throughout the literature there is no general rule on how to select the best. This can be partially explained by the fact that most publications consider different data for their analyses, and have different research targets, leading to different approaches on how to set up the MF. However, we do know everything about our $N$-body models, and this allows us to test the mass bin selection for multimass models. In particular, we want to understand what the minimum number of mass bins is to get a stable result and how to choose the bins.

For this analysis we use the $N$-body snapshot of simulation N1 at the time of $2.9\,\tau_{\rm{rh,0}}$. As we wanted to trace the overall evolution of the different star types, we opted against mixing them and therefore we give every star type at least one bin. This means that we have at least five mass bins, one each for the MS stars, the evolved stars333In this work, every post MS star which is not a remnant is called an ES. (ES), the white dwarfs (WD), NSs and BHs. Looking at the BHs, NSs and ESs we decided to not further split them into several bins given the fact that they have either a rather small range of possible masses and/or are to low in number to justify the split. This leaves us with MS stars and WDs which are both numerous and do have a large mass range: $0.1-0.83\,{\rm M}_{\odot}$ for the MS stars and $0.55-1.44\,{\rm M}_{\odot}$ for the WDs in our $N$-body snapshots.

First, we determine how the bin selection influences the results of the analysis. For this, we choose four different binning methods: a logarithmic binning, a linear binning, a binning where in each mass bin there is an equal number of stars and a binning where there is the same amount of mass in each mass bin. We fixed the number of WD mass bins to one and then for each bin type we calculated the multimass model repeatedly with increasing number of MS star bins. The general idea here is that with increasing number of mass bins, the overall parameters like $M_{\rm Cl}$, $r_{\rm h}$, etc., should converge to the value one would get for the ideal case, where each star has its own mass bin. We find that the results are almost independent of the way one chooses the binning: the number of MS mass bins needed to converge is the same and the difference between the different models are for all properties generally less than $5\%$. We therefore choose for the further analysis the logarithmic binning.

Then we determine the minimum number of bins needed to get stable results, as increasing the number of bins also increases the computation time of the models. We varied the number of mass bins of the MS stars and the WDs independently from each other. Here again, we see that with increasing number of mass bins the different quantities converge. We find that we need at least two WDs mass bins and at least four MS stars mass bins for the different quantities to converge. Increasing the number of bins any further does not improve the results (values are comparable within $5\%$). For our further analysis, we opt to use five MS stars mass bins and three WDs mass bins. Therefore, in total we consider eleven mass bins (MS: 5; ES: 1; WD: 3; NS: 1; BH: 1) to set up our multimass models. Tables 5–8 list the mass bins for all $N$-body snapshots used.

4.2 Artificial background population

Before we present the results, we discuss a particularity which we encountered in our analysis. The potential of limepy models is spherical, however, the true potential of the cluster is triaxial because of the effect of tides. Also the Lagrange points of the cluster, through which stars can escape (Fukushige & Heggie, 2000; Baumgardt, 2001; Küpper et al., 2010; Claydon et al., 2017), are not accounted for in the multimass model. Therefore, the models are not able to describe the objects near the critical energy correctly. From this, it follows that some of the objects which are unbound in the true potential are still found bound in our definition and cannot be described by the model correctly. These objects pose a problem because they drive the fit to unrealistic parameter values.

To cope with this problem, we introduced an artificial background population with a constant likelihood (i.e. uniform distribution) in phase space. We added the artificial background population with a total mass of around $1\%$ of the original cluster mass to the $N$-body snapshot. This background population has the same MF as the cluster. The upper limit for the maximal distance and velocity are chosen to be twice the maximal values from the original snapshot ($r_{\rm max}$ and $v_{\rm max}$).

We describe the likelihood function of the background model as

[TABLE]

where $\,M_{\rm tot}$ is the total mass of the snapshot including the artificial background and $M_{\rm Back}$ is the mass of the background only. The phase space volume $V$ is defined as

[TABLE]

The total likelihood of an object for a given model is calculated as

[TABLE]

When integrated over the whole phase space volume within $2v_{\rm max}$ and $2r_{\rm max}$ the first term equals to $M_{\rm Cl}/\,M_{\rm tot}$ and the second to $M_{\rm Back}/\,M_{\rm tot}$, giving a total likelihood of unity, as required.

4.3 MCMC results

We initiate the MCMC walkers in a randomly chosen sphere in parameter space. For some snapshots, we run several fits with different initial conditions to test for any divergence. We chose flat priors restricted mainly by currently observed values for the parameters and/or by the range in which they are considered physically valid. For the MCMC fitting, we started out with around $500$ walkers and found good fits for the $N$-body model N0 without BHs and NSs. For the other models, prominently those with BHs, converging fits were only achieved with at least $2000$ walkers. On average, each MCMC chain was run for $1000$ iterations and convergence was reached after around $300$ iterations, which we trimmed from the MCMC chains for the calculation of the best-fitting parameters. The MCMC chain took on average longer to converge in snapshots with BHs than in snapshots without. In some cases, we also had to adjust the emcee scale parameter $a$, which is generally set to $2$, to increase the acceptance rates (for details on how this affects the MCMC algorithm see the discussion in Foreman-Mackey et al. 2013, their Eq. 2).

In Figs. 5 and 6, we show the marginalized posterior probability distribution for each parameter as well as the 2D projections of the posterior probability distribution representing the covariance between the different fitting parameters for two MCMC runs. Figs. 5 and 6 show the results of the fitting to the $N$-body model N0 at $2.3\tau_{\rm{rh,0}}$, and to the $N$-body model N1 at $17.0\tau_{\rm{rh,0}}$, respectively. The obvious difference between the two models is that for model N1 there are two parameters, namely $r_{\rm a}$ and $\eta$ that do not converge to a single value. The stellar system in this particular $N$-body snapshot is isotropic: values of $r_{\rm a}$ larger than $r_{\rm t}$ generate isotropic models, equally likely to reproduce the data, and for this reason, the values of $r_{\rm a}$, and consequently of $\eta$, cannot be constrained.

Looking at the 2D projections of the posterior probability distribution in Figs. 5 and 6, one can see that they are nearly circular for most of the parameter pairs. This shows that when using the full phase space information of each star, degeneracies between the different parameters can be alleviated.

Tables 9–12 list the best-fitting parameters for all the $N$-body snapshots we considered.

5 Comparison of multimass models and $N$-body models

In the first part of our analysis, we compare the best-fitting multimass models with the results directly computed from the $N$-body snapshots for each model at all times.

5.1 Mass density profile

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 in each radial bin there are at least $30$ objects and each radial bin has a minimal radial width of $0.15{\,\rm pc}$. We assumed Poisson errors for the uncertainties of the binned data and define the position uncertainty by the 16th and 84th percentiles of the distribution of the positions of the objects in each bin. In Fig. 7, we compare the mass density profiles for the three models N1, N0.3 and N0 at four different times: $2.3\tau_{\rm{rh,0}}$ which is the first snapshot for each $N$-body model, $7.1\tau_{\rm{rh,0}}$, $16.5\tau_{\rm{rh,0}}$ and $26\tau_{\rm{rh,0}}$ the snapshot at the end of the clusters lifetime. For clarity we only show one mass bin per stellar type. We did not include the results of model N0.1, since they are similar to the results of model N0. Together with the best-fitting result from the multimass models, we also plot the results from the walker positions at the last iteration of the MCMC routine, reflecting the uncertainties of the results.

As can be seen from Fig. 7, the best-fitting multimass models reproduce the mass density profiles of the different mass components. Differences are only found in the outermost regions and innermost regions as well as for cases where the number of objects in a mass bin is low. For the outer regions, a difference is expected: as already discussed in Section 4.2, the models assume that the cluster is spherically symmetric, which is not the case in the $N$-body models as the clusters are slightly elongated due to tidal forces.

For the differences in the most central parts one sees that the mass density is underestimated for the heavier mass bins while for the lighter mass bins it is overestimated. This could be explained by the fact that these models are post core collapse, therefore their density profiles are slightly different from isothermal models (Lynden-Bell & Eggleton, 1980). Given that the differences in the centre are small, one can see that multimass models are able to describe post-collapse models.

The overall agreement between multimass and $N$-body models for the density profiles is consistent with the findings of Sollima et al. (2017) who fitted Michie–King models to observational data of NGC 5466, NGC 6218 and NGC 6981: for all three GCs, the multimass models reproduce the observed mass density profiles (see their fig. 6). Differences are only found in the outer regions, most likely for the same reasons as discussed above.

When one compares the different models in Fig. 7 at the same dynamical ages it can be seen that models without BHs are denser and the stars are found far more concentrated than in models with BHs. The BHs in the centre ‘push’ the lower mass stars out of the core, which results in a large core radius ($r_{\rm c}$) as well as a larger $r_{\rm h}$, an effect already studied by Peuten et al. (2016) for the cluster NGC 6101. In the evolution of model N0.3, one can see how the cluster changes when all BHs have been lost: the central regions get efficiently populated by the next lighter objects and the resulting mass density profile of the cluster looks as concentrated as the one from model N0 at that same dynamical age, leaving no clue about its diminished BH population.

5.2 Velocity dispersion

For the comparison of the velocity dispersion profiles in Fig. 8, we used the same snapshots as for the mass density profiles. For the calculations of the velocity dispersions, we are using a mass-weighted approach to make the values comparable to the values from the multimass model as they are calculated for the mean mass of each mass bin. The velocity dispersion is therefore calculated as:

[TABLE]

with

[TABLE]

the mass-weighted mean velocity for each component. The calculation 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 data from the $N$-body models. As with the mass density profiles, small difference can be seen in the outermost regions. In the plot of model N0.3 at $7.1\tau_{\rm{rh,0}}$, there is no value from the $N$-body snapshot for the BHs as there is only one BH left, in which case $\sigma$ is undefined.

When comparing the different models at the same dynamical age we find that in clusters with BHs, the velocity dispersions for the different mass bins are smaller than in clusters without BHs (see discussion in Section 6.5). As the cluster is losing its BH population (see for example the evolution of model N0.3), the velocity dispersions of the different mass bins increase to the values seen in model N0 which had all its BHs removed before the actual $N$-body evolution, again leaving no hint of the lost BH population. In Section 6.5, we will look again at this relation and discuss an explanation for this behaviour.

5.3 Radial anisotropy

In this section, we consider the anisotropy of the velocities. In Section 5.3.1, we consider the anisotropy profile within the cluster and in Section 5.3.2 we consider the global anisotropy of the cluster as a whole.

5.3.1 Anisotropy profiles

The anisotropy parameter $\beta$ is defined as (Binney & Tremaine, 1987):

[TABLE]

with $\sigma_{\rm r}$ the radial velocity dispersion and $\sigma_{\rm t}$ the tangential velocity dispersion. For $\beta<0$, the orbits are tangentially biased, for $\beta=0$ they are isotropic, for $0<\beta<1$ they are radially biased and for $\beta=1$ they are radial.

In Fig. 9, we compare the anisotropy profiles from the best-fitting multimass models for a selection of mass bins to the anisotropy profiles from the $N$-body snapshots. As the $\beta$ parameter is more affected by random scatter, we had to bin the data from the $N$-body snapshots differently than in the previous two plots. We varied the number of objects per bin, such that the average uncertainty in $\beta$ is $\leq 0.1$ and there are no more than $10$ radial bins per mass bin to not overcrowd the plots. For the ESs and the BHs, the average $\beta$ uncertainty was always well above $0.1$ which is why we only show one radial bin for each in all snapshots.

Comparing the predictions from the best-fitting multimass model with the results from the $N$-body data, we find that when the snapshot has some degree of radial anisotropy the multimass models qualitatively reproduce them. This can be seen best with the mass bins from the MS stars. Also differences between the best-fitting multimass prediction and the binned data can be seen at the outer regions of the cluster. When some of the mass bins are tangentially anisotropic the best-fitting model is isotropic as our multimass models cannot describe any other kind of anisotropy.

Looking at the data from the snapshots itself, we see that the heaviest mass bins become more radially anisotropic, while the low-mass bins become first isotropic and then tangential anisotropic. We will discuss the evolution of the anisotropy further in Section 6.7 when we analyse the best-fitting $\eta$ parameter.

5.3.2 Global anisotropy

To quantify the global anisotropy we use the parameter $\kappa$, that was introduced by Polyachenko & Shukhman (1981) and is defined as

[TABLE]

where $K_{\rm r}=0.5\sum_{i}m_{i}v_{{\rm r},i}^{2}$ is the radial component of the kinetic energy and $K_{\rm t}=0.5\sum_{i}m_{i}v_{{\rm t},i}^{2}$ the tangential component. For $\kappa=1$, the models are isotropic, for $\kappa>1$ they are radially anisotropic and for $\kappa<1$ they are tangentially anisotropic.

In Fig. 10, we compare the values of $\kappa$ obtained for the best-fitting multimass models and for the $N$-body snapshots. For the uncertainties of the $N$-body data, we used Poisson statistics. The best-fitting multimass models are able to qualitatively reproduce the overall behaviour of the $\kappa$ parameter. It can be seen that at later times, the low-mass stars are tangentially anisotropic, which cannot be reproduced by limepy. This also explains why the best-fitting value of $\kappa$ (and $\beta$) from the multimass models does not converge to unity immediately (or $\beta\simeq 0$): there are still radial orbits left and the tangential orbits are treated as isotropic, so the best-fitting results for $\kappa$ (and $\beta$) still indicate some radial anisotropy, resulting in a smoother transition from radial anisotropy to isotropy with respect to what is observed in the $N$-body data. Clusters that are dominated by tangential orbits are therefore, by construction, not well reproduced by the multimass models used in our study (see also Sollima et al. 2015).

Looking at the $N$-body data, we see that $\kappa$ of the lowest mass bins typically goes down with time, while $\kappa$ of the heaviest mass bins goes up. For model N1 with BHs, this change is faster than for model N0 with no BHs and NSs. We refer the reader to Section 6.7 for a further analysis.

6 Analysis of model parameters

We focus here on the best-fit ting parameters resulting from our fitting procedure. The two model scale parameters can be computed directly from the $N$-body data, therefore we use them to assess the quality of the multimass models. For the other five model parameters, we are analysing their evolution to see whether they can give us some further insights into the clusters. Furthermore, we also discuss the evolution of two additional quantities ($r_{\rm t}$ and $\kappa$) to assess the quality of the models.

6.1 Total cluster mass

In Fig. 11, we plot the $M_{\rm CL}$ from the best-fitting multimass model divided by the true cluster mass as measured in the $N$-body model for all four $N$-body models throughout their cluster lifetime. As can be seen in this figure, the best-fitting value is always within $1\%$ of the $N$-body value but almost none of these are consistent within the $1\sigma$ uncertainties: there is some systematic error in the multimass models which is not accounted for yet. However, given that the difference is always smaller than $1\%$ this effect is negligible. The results are comparable to the single-mass results from Zocchi et al. (2016) where an accordance within $5\%$ of the true values was found.

Sollima et al. (2015) found that single-mass models underestimated the mass of an $N$-body system by $50\%$, depending on its dynamical state. There are several differences between their study and ours which could lead to such different results. They are using simulated observations as input 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 an underestimation of the mass by a factor of two or three, especially for metal-rich GCs. Also the models used in Sollima et al. (2015) have less parameters than the ones used here, as for example they do not incorporate the variable truncation parameter $g$. Zocchi et al. (2016) showed for single-mass models that the total mass is better recovered when allowing $g$ to be free. We discuss the effect of $g$ in Section 6.4.

6.2 Half-mass radius

The second scale parameter that can be computed from the $N$-body models is $r_{\rm h}$: In Fig. 12, we plot $r_{\rm h}$ from the best-fitting multimass models divided by the true value as computed from the $N$-body data, for all four clusters throughout their lifetime. Again we find good agreement, within a few percent. As with $M_{\rm Cl}$, only a few of the data points show agreement within $1\sigma$ uncertainties. These results are comparable to the result for the single-mass case by Zocchi et al. (2016) who found an agreement within $7\%$.

6.3 Dimensionless central potential

Now that we have shown that multimass models can reproduce the most important cluster properties, we focus on analysing the other fitting parameters. First, we look at the dimensionless central potential $W_{0}$ for which the evolution over the whole lifetime for the four $N$-body models is plotted in Fig. 13. As discussed in Section 2, the $W_{0}$ value in multimass models represents the dimensionless central potential of a hypothetical mass group with a mass equal to the global mean mass. As the global mean mass increases from $(0.36\pm 0.01)\,{\rm M}_{\odot}$ to $(0.8\pm 0.2)\,{\rm M}_{\odot}$ during the evolution of the four $N$-body models, the $W_{0}$ values do not refer to the same stellar population and this needs to be kept in mind when comparing $W_{0}$ values of different $N$-body models and/or at different times. Using the central density weighted mean mass instead of the global mean mass has other issues, as can be seen for example in the snapshots of model N0.3 at $4.7\tau_{\rm{rh,0}}$ or in model N1 at $26.7\tau_{\rm{rh,0}}$ in Fig. 14. In both cases, the number of BHs decreases to (almost) zero reducing the central density weighted mean mass more than the global mean mass, explaining the more significant change of the $W_{0}$ values in this figure.

It is also possible to use the values of $W_{0}$ obtained for different mass bins for a comparison. As an example we consider the ESs mass bin, because not only does the average mass of the ESs not vary in our models [$\bar{m}_{\rm ES}=(0.821\pm 0.006)\,{\rm M}_{\odot}$], but it also represents the objects that are easiest to observe. In Fig. 15, we plot the evolution of the $W_{0}$ values of the ESs for all four $N$-body models over their entire lifetime. The uncertainties are estimated by calculating the $W_{0}$ values of the ESs for the last 10 iterations of all the walkers and then using the values from the 16th and 84th percentiles as $1\sigma$ uncertainties.

Looking at Figs 13 and 15, one can see the effect of BHs as discussed in Section 5.1: in clusters with BHs, the other stars are pushed outwards and the cluster appears less concentrated with a low value of $W_{0}$ for the mean mass stars and ESs. While clusters without BHs instead appear much more concentrated as the normal stars can occupy the centre and therefore have a higher $W_{0}$ value for the mean mass stars and for the ESs. Therefore, clusters with low $W_{0}$ value for the observable stars are much more likely to be hosting a BH population than clusters with a high $W_{0}$ value (Merritt et al., 2004; Mackey et al., 2008; Peuten et al., 2016).

6.4 Truncation parameter

The truncation parameter $g$ provides an indication of the effect of external tides on the stellar system. In Fig. 16, we plot the evolution of $g$ for the four clusters over their lifetime. The evolution is similar for all four $N$-body models: at the beginning $g$ is around $1.5$ which represents a model in between a Wilson (1975) model ($g=2$) and a King (1966) model ($g=1$). As the clusters fill their Roche volume, the tides interact with the clusters, stripping their outermost stars and thereby making the truncation in energy space steeper. This evolution is reflected in the truncation parameter $g$ decreasing as the cluster evolves, converging at the end of the lifetime to a value of around $g\approx 0.73$ which represent a model between a King (1966) and a Woolley (1954) model ($g=0$).

The results are comparable to the single-mass model findings of Zocchi et al. (2016), though they start with a higher truncation parameter, due to the fact that they start with a smaller initial $r_{\rm h}/r_{\rm J}$ ratio ($=0.01$) than we do ($r_{\rm h}/r_{\rm J}=0.7$). We must note that we only use bound objects in our analysis, which might lead to smaller values of the truncation parameter as in the outer regions of the clusters ($0.8r_{\rm J}-r_{\rm J}$) most objects are energetically unbound (Claydon et al., 2017).

As before, we see a difference between the model with BHs (N1) and the models which are mostly BH free (N0, N0.1 and N0.3): at the beginning the value of $g$ decreases faster for model N1 than for the other three models and therefore also converges quicker to its final value. This behaviour in the first half of evolution is not too surprising given that $r_{\rm h}$ of that model is roughly twice as large as the others, thereby the impact on the steepness of the truncation is stronger (see Section 3.3).

McLaughlin & van der Marel (2005) found that Wilson models are equally good or better in describing a sample of Galactic and extragalactic clusters than King models. With our results of the evolution of $g$, one can interpret McLaughlin & van der Marel (2005) findings that clusters with large $g$ are still dynamically expanding towards filling the Roche volume (see also Carballo-Bello et al. 2012).

6.5 Mass segregation

In Fig. 17, we plot the evolution of the best-fitting value for $\delta$ during the whole lifetime of the four $N$-body models. In the initial stages of their evolution, all four $N$-body models are still in the process of segregation, as they were set up without any primordial mass segregation. Over the course of evolution of the clusters, the value converges to around $\delta=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 $\delta\gtrsim 0.5$, however the results are compatible with $0.5$ within $3\sigma$. Sollima et al. (2015) found that for some of their late $N$-body snapshots, the multimass models underestimate the amount mass segregation. The same is also found for the best-fitting multimass model to the observations of NGC 6218 in Sollima et al. (2017). However, we do not find this in these models.

To further analyse the behaviour of $\delta$ over the course of the cluster evolution, we additionally plot in Fig. 18 the central velocity dispersion for the different mass bins from the $N$-body data together with the predicted central velocity dispersion from the best-fitting multimass model (see also Fig. 8). Given the results from Section 5.2, it is no surprise that the best-fitting multimass models are able to reproduce the true values.

Despite the fact that the best-fitting models have $\delta\approx 0.5$, this does not mean that the multimass models are in a state of energy equipartition, as can be seen in Fig. 18. This was already pointed out by Merritt (1981) and Miocchi (2006) as well as by GZ15. In a mass-segregated multimass model the following relation

[TABLE]

holds true with $s_{j}$ and $s_{i}$ being the velocity scale of two different mass bins. For the 1D velocity dispersion at the centre the relation:

[TABLE]

only holds true for $W_{0}\rightarrow\infty$. In the $N$-body models, we are studying here $W_{0}\ll\infty$ and therefore $\sigma_{1d,j0}<s_{j}$, which means that our multimass models are never in a state of energy equipartition despite having $\delta=0.5$.

Furthermore Bianchini et al. (2016) showed, using Monte Carlo simulations, that in clusters only objects above a certain critical mass $m_{\rm eq}$ can be in energy equipartition. The value of $m_{\rm eq}$ depends on the mass spectrum of the cluster, and is larger for cluster with a wider mass spectrum, i.e. when BHs are retained. This trend can also be seen in Fig. 18, where the number of mass bins following the $\sigma_{1d}(0)\propto m_{j}^{-1/2}$ relation, which are therefore in energy equipartition, is only greater than one for the models without BHs. For models with BHs, only the BHs can be in energy equipartition as they are the only ones which have a mass greater than $m_{\rm eq}$. This leads to the largest part of the other objects in these clusters having a smaller spread in the velocity dispersions, which is the reason for the reduced mass segregation in the observable stars in clusters with BHs.

Looking at the results in Fig. 17, we can conclude that setting the mass segregation parameter to a fixed value of $\delta=0.5$ as in its initial formulation of the multimass models by Da Costa & Freeman (1976) is indeed justified to model all but the youngest clusters. As those are not yet fully mass segregated, the value of $\delta$ must therefore be smaller, something also found by Sollima et al. (2015); Sollima et al. (2017) for young clusters.

6.6 Anisotropy radius

The last two fitting parameters are coupled together as they both determine $r_{{\rm a},j}$ for each mass bin as can be seen in equation (3).

First we focus on $r_{\rm a}$. In Fig. 19 we plot $r_{\rm a}/r_{\rm t}$ of the best-fitting model, for all four $N$-body models during their lifetime. If $r_{\rm a}/r_{\rm t}\gtrsim 1$, the cluster is isotropic (see Section 2). Considering the definition of $r_{{\rm a},j}$ (equation 3), it follows that even if $r_{\rm a}$ is well above $r_{\rm t}$ for some mass bins $r_{{\rm a},j}$ can still be below $r_{\rm t}$ and therefore these mass bins still show some degree of radial anisotropy.

Looking at Fig. 19, one can see that the model which retained all its BHs (N1) behaves differently from the other models. Model N1 is only radially anisotropic at the beginning of the lifetime and quickly becomes isotropic. The model without BHs (N0) loses its radial anisotropy more slowly and only at the end of its lifetime it becomes isotropic/tangentially anisotropic. The two models in between (N0.3 and N0.1) behave at the beginning differently: as long as they still have some BHs left their $r_{\rm a}$ value evolves independently, but after all BHs are lost their $r_{\rm a}$ value drops to the $r_{\rm a}$ value of N0 at the same dynamical age and from then on essentially follows the $r_{\rm a}$ evolution of model N0.

These results obtained for the BH-free clusters are comparable, albeit with smaller magnitude, to the results found for the single-mass case by Zocchi et al. (2016). They showed that the anisotropy radius is monotonically decreasing till the cluster reaches core collapse after which it is monotonically increasing, and it eventually becomes so large that the corresponding model is isotropic.

Before we discuss the possible physical reasons behind the evolution of $r_{\rm a}$, we first have a look at the fitting parameter $\eta$, which is needed to include a dependence of the anisotropy radius on the mass.

6.7 Mass-dependent anisotropy

The anisotropy parameter $\eta$ is a novel fitting parameter in multimass models. In Fig. 20, we plot the values of $\eta$ for the four clusters over time. We only plot this for the snapshots showing some degree of anisotropy. The most important feature is that $\eta$ evolves for all clusters from a value of $\eta\approx 0.5$ at the beginning of the clusters lifetime to a value of $\approx-2.5$ at the end of their lifetime. The model which retains BHs throughout its lifetime (N1) stops the evolution earlier with a value of $\eta\approx-0.5$.

The evolution of $\eta$ is not surprising given what we already saw in Section 5.3, where we showed that the amount of radial anisotropy is decreasing in the low-mass bins and is increasing in the high-mass bins. This is reflected in the development of $\eta$ changing from a positive value to a negative one over time. This trend is comparable to what Sollima et al. (2015) found in their analysis of the W5rh1R8.5 $N$-body model from Baumgardt & Makino (2003): they found that the low-mass stars, which are preferentially located in the cluster outer regions due to mass segregation, become tangentially anisotropic. The reason for this behaviour is that interactions occurring in the cluster centre kick stars into the cluster halo on to radial orbits (Lynden-Bell & Wood, 1968; Spitzer & Shull, 1975). As stars on radial orbits reach the cluster boundary with positive velocity, they can escape the cluster more efficiently, thereby depleting the low-mass population from stars with radial orbits, leaving only the stars with tangential orbits in the cluster.

To test the relevance of $\eta$, we rerun fits to model N0 but this time fixing $\eta=0$ comparable to the original formulation by Gunn & Griffin (1979). In these fits, the most obvious difference is that with increasing time, and therefore also with a higher absolute value of $\eta$, the uncertainties of $r_{\rm a}$ increase up to five times the value recovered in the fits with a non-fixed $\eta$ value. Therefore, the introduction of $\eta$ improves the ability to describe the data.

6.8 Truncation radius

At the end of our comparison, we look at two quantities on which we do not fit but which get computed by the multimass models and which can be calculated for the $N$-body models. In Fig. 21, we plot $r_{\rm t}$ divided by $r_{\rm J}$ as determined in Section 3.2 for the four $N$-body models over their whole lifetime. The values of $r_{\rm t}$ and their uncertainties were computed using all the walker positions of the last ten iterations of the MCMC runs. This figure shows that $r_{\rm t}$ stays within $3\%$ of the computed $r_{\rm J}$ and in all but two cases the results are consistent within their $1\sigma$ uncertainties. The largest discrepancies 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 multimass models are able to reproduce $r_{\rm t}$ accurately.

6.9 Global anisotropy parameter

In Fig. 22, we look at the evolution of the global value of $\kappa$. Here, we have plotted the comparison between the best-fitting value inferred using the walker positions of the last 10 iterations of the MCMC runs from each snapshots to the one calculated from the $N$-body snapshots directly. These are calculated applying equation (18) to all objects in the $N$-body model. For the uncertainties of the $N$-body data, we used Poisson statistics. As before the best-fitting multimass models qualitatively reproduce the overall trend as long as the $N$-body model is radially anisotropic. When the $N$-body snapshots becomes tangentially anisotropic, the best-fitting multimass models are isotropic. Hence, the multimass $\kappa$ values still shows some radial anisotropy where the true cluster is already dominated by tangentially anisotropic orbits. For all models $\kappa<1.7\pm 0.25$, from which it follows that all models are stable against radial orbit instability as discussed in Polyachenko & Shukhman (1981).

7 Discussion and conclusion

In this study, we assessed the validity of the multimass anisotropic models provided by the limepy software (GZ15) fitting them to $N$-body models. We find that the $N$-body models are well described by multimass models, a result which is fortunate given the long list of observational studies using multimass models to analyse GCs (see Section 1). Zocchi et al. (2016) showed for the single-mass case that the limepy models are able to describe clusters at all evolutionary phases. Although the agreement is not perfect, the systematic differences are negligible for most applications and parameters of interest (see also the discussion in Section 6.1).

Our comparison shows that the best-fitting total cluster masses are off by no more than $1\%$ from the true value as computed from the $N$-body snapshots. The best-fitting cluster half-mass radius is reproduced within $5\%$ and the truncation radius is reproduced within $3\%$.

We find that the mass density and velocity dispersion profiles of the different mass bins are well reproduced by the multimass models. If the $N$-body snapshot is radially anisotropic then the multimass models are generally able to reproduce it.

We show that in the $N$-body models, regardless of initial BHs and NSs retention, the truncation parameter $g$ evolves from roughly $1.5$ to about $0.7$. The general trend can be explained by the tidal effects stripping the loosely bound stars. We find that the best-fitting mass segregation parameter $\delta$ converges to a value close to $0.5$ for our $N$-body models, which is the value used in the original formulation by Da Costa & Freeman (1976). Only for young clusters which are not yet mass segregated is the best-fitting value smaller and for models with BHs it is $\sim 0.4$.

The newly introduced $\eta$ parameter shows that the anisotropy radius is mass dependent and that this mass dependence changes in our $N$-body models over time from $\eta=0.5$ where the lighter stars are more radially anisotropic to $\eta=-2.5$ ($\eta=-0.5$ for the model which initially retained all its BHs) where the heavy objects are more radially anisotropic. We find in this study that the effects which influence the anisotropy radius are more mass dependent than initially thought and therefore $\eta$ is another relevant parameter when analysing radial anisotropy with multimass models. Furthermore, we find that clusters with a BH population can be tangentially anisotropic for most of their lifetime.

The $W_{0}$ parameter for the observable ESs is lower for the clusters with BHs than for the clusters without BHs. Therefore, clusters which still harbour a stellar-mass BH population should appear less dense when looking at the observable stars. $N$-body simulation N0.1, which loses its BHs within its first $\tau_{\rm{rh,0}}$, does not show any strong differences 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 an adequate tool to study the global properties of GCs, as the results from the comparison with $N$-body models show a good agreement with their properties inferred from multimass models.

Acknowledgements

This research was done as part of the Gaia Challenge444http://astrowiki.ph.surrey.ac.uk/dokuwiki/doku.php?id=start, whose goal is to compare and improve mass modelling techniques in preparation of data releases of the ESA Gaia mission. We thank the organizers of the Third Gaia Challenge Workshop (Barcelona, 2015) for a very productive meeting. We thank the anonymous referee for constructive suggestions. We are grateful to Anna Lisa Varri and Antonio Sollima for interesting discussions. We are grateful to Sverre Aarseth and Keigo Nitadori for making nbody6 publicly available, and to Dan Foreman-Mackey for providing the emcee software and for maintaining the online documentation; we also thank Mr. Dave Munro of the University of Surrey for hardware and software support. MG acknowledges financial support from the Royal Society (University Research Fellowship) and AZ acknowledges financial support from the Royal Society (Newton International Fellowship). MP, MG and AZ acknowledge the European Research Council (ERC-StG-335936, CLUSTERS). VHB acknowledges support from the Radboud Excellence Initiative Fellowship.

Appendix A $N$-body data and MCMC results