Issue 
A&A
Volume 617, September 2018



Article Number  A69  
Number of page(s)  11  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/201832919  
Published online  19 September 2018 
The black hole retention fraction in star clusters
^{1}
Astronomical Institute of Charles University, Prague, Czech Republic
email: pavlik@sirrah.troja.mff.cuni.cz
^{2}
Observatory and Planetarium of Prague, Prague, Czech Republic
^{3}
HelmholtzInstitut für Strahlen und Kernphysik (HISKP), Universität Bonn, Bonn, Germany
^{4}
European Southern Observatory, Garching bei München, Germany
email: tjerabko@eso.org
^{5}
School of Mathematics and Physics, University of Queensland, St. Lucia 4072, Australia
Received:
28
February
2018
Accepted:
28
May
2018
Context. Recent research has been constraining the retention fraction of black holes (BHs) in globular clusters by comparing the degree of mass segregation with Nbody simulations. They are consistent with an upper limit of the retention fraction being 50% or less.
Aims. In this work, we focus on direct simulations of the dynamics of BHs in star clusters. We aim to constrain the effective distribution of natal kicks that BHs receive during supernova (SN) explosions and to estimate the BH retention fraction.
Methods. We used the collisional Nbody code nbody6 to measure the retention fraction of BHs for a given set of parameters, which are: the initial mass of a star cluster, the initial halfmass radius, and σ_{BH}, which sets the effective Maxwellian BH velocity kick distribution. We compare these direct Nbody models with our analytic estimates and newest observational constraints.
Results. The numerical simulations show that for the onedimensional velocity kick dispersion σ_{BH} < 50 km s^{−1}, clusters with radii of 2 pc and that are initially more massive than 5 × 10^{3} M_{⊙} retain more than 20% of BHs within their halfmass radii. Our simple analytic model yields a number of retained BHs that is in good agreement with the Nbody models. Furthermore, the analytic estimates show that ultracompact dwarf galaxies should have retained more than 80% of their BHs for σ_{BH} ≤ 190 km s^{−1}. Although our models do not contain primordial binaries, in the most compact clusters with 10^{3} stars, we have found evidence of delayed SN explosions producing a surplus of BHs compared to the IMF due to dynamically formed binary stars. These cases do not occur in the more populous or expanded clusters.
Key words: methods: numerical / methods: analytical / galaxies: star clusters: general / stars: black holes
© ESO 2018
1. Introduction
Black holes (BHs) and their retention fraction in star clusters play an important role in the evolution of the clusters, and are relevant for other astrophysical fields including stellar evolution, the formation of intermediatemass/supermassive BHs and predictions of gravitational wave events (e.g. Belczynski et al. 2002; Favata et al. 2004; Mackey et al. 2007; Banerjee et al. 2010; Mapelli et al. 2011, 2013; Morscher et al. 2013; Ziosi et al. 2014; Banerjee 2017, 2018; Repetto et al. 2017). The retained number of BHs also constrains the maximum mass of potentially formed intermediatemass/supermassive BHs in massive stellar systems such as in globular clusters (GCs) and ultracompact dwarf galaxies (UCDs; Jeřábková et al. 2017).
Black holes can be observed and quantified either directly via gravitational waves or indirectly via (i) mass accretion (which requires the presence of gas and/or stellar companions), (ii) gravitational lensing, or (iii) mass segregation in star clusters (Baumgardt & Sollima 2017; Weatherford et al. 2017). Method (iii) has been successfully used to establish an upper limit on the BH retention fraction in GCs to be 50% or less (Peuten et al. 2016; Baumgardt & Sollima 2017). However, the constraints on the initial retention fraction of BHs remain weak. By the initial retention, which is what this study is concerned with, we mean the fraction of BHs that are retained in the star cluster by the time of the last core collapsed supernova (SN) leaving a BH. This corresponds to a time scale of about 12 Myr after the formation of the cluster.
An asymmetric SN explosion gives the newly formed remnant an initial momentum, p_{BH}, to compensate for the excess momentum of the stellar envelope going in the opposite direction, p_{env−}, that is,
where p_{env+} is the momentum vector of the envelope in the same direction as p_{BH}. The velocity kick, v_{kick}, is then the initial remnant’s velocity derived from its momentum and mass (Lyne & Lorimer 1994).
Here we provide a systematic direct Nbody survey of the initial retention fractions of BHs for different assumptions on the kick velocities for a variety of star cluster radii and masses, that is, in the range from 10^{3} to 10^{5} stars. We also compare these numerical results with our own analytic estimates on the retention fraction of BHs. Finally, we use this study to extrapolate to the initial retention fraction in larger systems, that is, GCs and UCDs.
2. Numerical models
We performed over 1500 Nbody simulations of isolated star clusters. Clusters with lower numbers of stars (N < 25k) were evolved with a singlecentralprocessingunit (CPU) nbody6 integrator. Clusters with more stars were integrated with a parallel nbody6.sse version (Aarseth 2003; Nitadori & Aarseth 2012)^{1}. All the initial conditions were set up with a random seed by the integrator as follows.
The models presented contain 1k, 3k, 10k, 25k, 50k, and 100k stars; see Table 1 for their detailed integration parameters. We assumed the canonical initial mass function (IMF, Kroupa 2001; Kroupa et al. 2013) with masses in the range from 0.08 M_{⊙} to 100 M_{⊙} and the distribution of positions of stars according to Plummer (1911) and Aarseth et al. (1974). For each model, we used three values of the initial virial radius, r_{v} = 0.5, 1, and 2 pc; the value of the initial Plummer radius is and the halfmass radius is r_{h} ≈ 1.305 r_{P} (Kroupa 2008). Neither primordial binaries nor a gas component were included in the models.
For the stellar evolution, we invoked the algorithm by Hurley et al. (2000) for single stars and Hurley et al. (2002) for binary stars. For that we assumed a metallicity [Fe/H] = −1.30 which corresponds to the average metallicity of GCs in the Galaxy (Baumgardt & Sollima 2017). We have also performed one additional calculation of a vastly expanded star cluster (r_{v} = 20.0 pc) with the same metallicity and IMF and no primordial binaries in order to estimate the time scale needed for a single star to evolve into a BH. In an extended cluster like this, dynamical effects and binary evolution may be neglected. The last BH appeared at ≈11.7 Myr from the beginning of the integration. Therefore, all our simulations were closely followed over that time.
Finally, in each setup, we assumed a different value for the onedimensional (1D) dispersion, σ_{BH}, the Maxwellian distribution from which we drew the initial BH kick velocities, v_{kick,i}. It has been suggested, for example by Jonker & Nelemans (2004) or Repetto et al. (2012), that BHs may receive kicks as high as neutron stars (NSs), therefore the highest value of σ_{BH} comes from the best fit of the velocity distribution of the observed NSs (Hansen & Phinney 1997), that is, σ_{BH} = 190 km s^{−1}. We set the lowest value of σ_{BH} to 3 km s^{−1} and a moderate value to 50 km s^{−1}. In the case of our 1 pc models, we also included σ_{BH} = 10 and 25 kms^{−1}; see Table 1. The kick velocity is composed of three random deviates chosen from a Gaussian distribution with a 1D dispersion σ_{BH} (e.g. Kroupa 2008).
In this work, we do not treat a detailed kick velocity mechanism (we simply assume a Maxwellian distribution). There are currently different models of how the kicks are produced during the SN explosions: After a SN explosion, the initial velocity kick, v_{kick,i}, may be reduced by the mass that falls back onto the remnant. One possibility is that the kick velocity scales with the ratio between the mass of the envelope, m_{env}, and the mass of the star before the SN explosion as
where m_{BH} is the remnant’s mass (Aarseth 2003). However, according to a recent study by Belczynski et al. (2008) the remnant’s velocity should be scaled by a fraction of the envelope mass, Δm_{env}, that is,
Fryer et al. (2012) use the same prescription for the final magnitude of the velocity kick and they also argue that the fallback mass is proportional to the mass of the envelope, although they determined the mass differently. In either case, the recoil is given by the momentum conservation due to an asymmetric spatial distribution of the mass of the envelope; see Eq. (1), and the final velocity is scaled by the mass that falls back (Δm_{env}). Therefore, the total mass of the star before a SN, as in Eq. (2), may not be a valid scaling parameter for the fallback since it cannot help us establish by how much the fallback fraction of the envelope actually slows the remnant. Because in our case the masses of the BHs are similar, the effect of fallback on the final kick velocity is comparable to evaluating different values of σ_{BH}. Therefore, in our models, there is no rescaling of the kick velocity, that is, v_{kick} = v_{kick,i}.
Parameters of the starcluster models used in this paper.
In our calculations, if a star explodes as a SN while bound in a binary system, the remnant does not receive a kick. This is the case in less than 2% of stars that lead to SNe, so this does not significantly affect the results. Because the kick velocity dispersion σ_{BH} is an assumed quantity, it is to be interpreted as an effective kick velocity dispersion. This means that any physical process that changes the actual kick, v_{kick,i}, to the velocity of the BH when it is free streaming in the cluster (e.g. after dissociating from a binary) is not considered here.
We note that, according to Belczynski et al. (2010, 2016), for example, BH masses may grow larger for lower metallicities than according to Hurley et al. (2000). This, however, does not have any effect on our present calculations of the kick velocity because in our models, v_{kick} does not depend on the BH mass. The only possible effect could be through the cluster’s expansion as a result of larger masses of the lost BHs. In the models presented, about 6% of the mass of the cluster can be lost if all BHs escape. This value could increase by a maximum factor of two if we consider the approach of Belczynski et al. (2010). Although this would be useful to study in the future, it should not have a significant effect on our results; see, for example, Fig. 1 where we compare the cluster’s expansion by means of the halfmass and tidal radii. In the models with σ_{BH} = 190 km s^{−1}, almost all BHs escaped, while in the models with σ_{BH} = 3 km s^{−1} only a small fraction did (see Sect. 4). Nevertheless, no significant difference between those two plots can be seen.
Fig. 1. Evolution of the halfmass radii (r_{h} ≈ 0.769 r_{v}, solid lines) and tidal radii (dotted lines) in our models. We compare models with the lowest and highest σ_{BH}. Models with different N are colourcoded as given by the key at the top of the figure. 

Open with DEXTER 
3. Methods
3.1. Retention fraction
In each realisation of our models, we tracked the positions of all newly formed and existing BHs. The retention fraction, hereafter denoted as η_{BH}, corresponds to the fraction of BHs that are inside a certain radius, for example, the halfmass radius or the cutoff radius.
Both of these radii are derived with respect to the density centre, provided by nbody6 (according to Casertano & Hut 1985). The cutoff radius, r_{t}, is taken to be the tidal radius (e.g. Binney & Tremaine 1994)
where M_{C} is the mass of the cluster within this radius, r_{G} is the distance of the cluster from the centre of the Galaxy (assumed to be 5 kpc), and M_{G} is the mass of the Galaxy comprised in the radius r_{G}. According to Faber & Gallagher (1979) and BlandHawthorn & Gerhard (2016), we took M_{G} ≈ 5 × 10^{10} M_{⊙}. At each time step, the halfmass radius is calculated from the stars that are bound to the cluster, that is, up to r_{t} from Eq. (4).
Certain dynamical effects, for example, close encounters of two stars or single stars with binaries, are able to significantly reduce the fraction of BHs that are in the cluster on a time scale of hundreds of millions of years. Other effects, such as dynamical friction, can slow down escaping BHs sufficiently to maintain the BH population inside a cluster for longer. In order to reduce the influence of these effects as much as possible, we evaluate the BH retention soon after the last BH has formed in each model (given the metallicity, this is at approximately 12 Myr). We are aware that the latest BHs with velocities drawn from a Maxwellian distribution with σ_{BH} = 3 km s^{−1} may not be able to escape from the cluster by that time. On the other hand, their v_{kick} would barely exceed v_{esc}, so our results should hold (see Sect. 4 and Fig. 2).
Fig. 2. Average ratio between the number of BHs that received v_{kick} > v_{esc} during a SN explosion (N_{BH,kick}) and the number of BHs that are outside the tidal radius at 12 Myr (N_{BH,tide}). Each colour corresponds to a different kick velocity dispersion σ_{BH} listed above the plots. The error bars correspond to the Poisson uncertainties; see Eq. (9). 

Open with DEXTER 
3.2. Analytic estimate
To estimate the retention fraction analytically, we generated a set of star clusters with the canonical IMF in the same mass range as our Nbody models. First, we assume that, at the time of the kick, the systems are Plummer (1911) models with the virial radius r_{v} and that the kick velocities, v_{kick}, follow a Maxwellian distribution,
with a velocity dispersion σ_{BH}. The values of σ_{BH} used in these estimates are the same as we used in the initial conditions in the above Nbody models.
The BH retention fraction is the fraction of BHs that do not escape from the cluster. Here, we assume that the stellar remnants are lost if their velocity after a SN kick is larger than the escape velocity at their current radius r from the centre of the star cluster. This limitvelocity is defined as
where M_{C} is the mass of the cluster, G is the gravitational constant, and r_{P} is the Plummer radius from the initial conditions of our models. The predicted ratio of retained BHs is then given by
4. Results
4.1. Retention fraction
If no stellar dynamics is taken into account, for a given effective kick velocity dispersion (σ_{BH}), the value of η_{BH} should be stationary and depend only on the number of heavy stars in the model, that is, the initial cluster mass (given that we use the same IMF and stellar evolution in all models). In Fig. A.1 and Table A.1, we document the evolution of the mean retention fraction in our models. The results show that almost all BHs are kicked out of the clusters in our simulations with σ_{BH} = 190 km s^{−1} and also that the retention fraction increases with a decreasing σ_{BH}. In Fig. 1, we compare the temporal evolution of the halfmass and tidal (or cutoff) radii of our clusters with the lowest and highest σ_{BH}. The halfmass radii calculated from these two classes of models show not more than a marginal difference in the overall evolution. As the mass contained in the BHs is around 6% of the total mass, it is expected that the whole cluster would not feel their absence.
We have also evaluated dynamical effects on the retention of BHs. In Fig. 2 we plot the ratio
Here, N_{BH,kick} is the number of BHs that received a kick with a higher velocity than the escape velocity at their current distance from the cluster centre from Eq. (6), that is, v_{kick} > v_{esc}. The denominator, N_{BH,tide}, is the number of BHs that are outside the tidal radius of the cluster at 12 Myr, that is, the approximate time necessary for BHs to form in our models. We also include the Poisson uncertainties defined as the square root of the value. Based on the propagation of uncertainty, the plotted error bars correspond to
When confronted with the final η_{BH}, this ratio indicates one of three scenarios: all the BHs have been ejected because of a SN kick (Q = 1), some of the kicked BHs were dynamically retained inside the cluster (Q > 1), or how many BHs, not retained in the cluster, have been expelled dynamically through encounters with other stars or other BHs (Q < 1).
We deduced that natal kicks are a dominant factor influencing the resulting initial retention fraction, η_{BH}, in our models with σ_{BH} ≥ 10 km s^{−1}. The only models where dynamical processes have a major influence on η_{BH} are those with σ_{BH} = 3 kms (black lines in Fig. 2), where Q < 1. Nonetheless, even in those models, the retention fraction in the tidal radius is 82% or higher (see Table A.1). Ejecting a BH merely through dynamical processes is, therefore, rare on such a short time scale.
According to the single stellar evolution algorithm by Hurley et al. (2000) which is parametrised only by mass and metallicity of individual stars, the last SN explodes at ≈11.7 Myr. For the small N systems, we also see a systematic shift to a BH production beyond this time (in Fig. 3). The reason is that the IMF sampled for 1000 stars gives only a couple of massive stars (sometimes only one). Due to dynamical processes, those few massive stars sink to the centre of the cluster and form a binary star that has a very small chance of being disrupted by the lessmassive stars surrounding it. The binary stellar evolution algorithm (Hurley et al. 2002) introduces additional parameters, for example, mass transfer, accretion, and collisions, which lead to a rejuvenation of the star, causing a delayed SN explosion. This process is particularly pronounced, even by a factor of two (see the top panel of Fig. 3), for the r_{v} = 0.5 pc clusters in which the binaries that form are sufficiently compact for the rejuvenation to occur. Individual cases, for example, the 100k model with r_{v} = 1.0 pc and σ_{BH} = 190 km s^{−1} (orange circle in the middle panel of Fig. 3), or the 10k model with r_{v} = 0.5 pc and σ_{BH} = 3 and 190 km s^{−1} (black and orange circles in the middle panel of Fig. 3, respectively), also delayed the last BH formation to above 12 Myr. We interpret those as being due to rejuvenated stars as well.
Fig. 3. Time of BH formation in the presented models. Circles correspond to the time of the last BH formation (i.e. f_{BH} = 100%), triangles represent the time when a fraction f_{BH}> 97% of BHs formed. For reference, the dotted line marks 11.7 Myr, as deduced from the Hurley et al. (2000) algorithm. 

Open with DEXTER 
In the most compact models, that is, r_{v} = 0.5 pc with 1k (and also 3k) stars, we see an overproduction of BHs compared to other models (even with the same number of stars) – compare the top panel of Fig. 4 with the lower two. Any random effect of sampling of the initial conditions may be ruled out because we have done several realisations of each model (hundreds in the case of 1k models). We explain this increment of BHs by the stellar evolution in binary stars. The effect of delayed SN explosions of massive stars has been systematically studied, for example by Podsiadlowski et al. (1992) and De Donder & Vanbeveren (2003), or more recently by Zapartas et al. (2017). In our models, especially in those compact and not very populous clusters, highmass binaries are easy to form and difficult to destroy. We demonstrate this in Fig. 5 where we plot the mean lifetime of a BH progenitor in a binary before becoming a BH. In the most compact clusters (top panel), especially for a low number of stars, the binary evolution is clearly more significant than in more massive or larger clusters (middle and bottom panels). Due to mass transfer in the binary, even an initially less massive star (which would have ended as a NS) can accrete enough mass from its more massive companion to eventually become a BH. Therefore, we end up with an additional BH. This effect is not as profound in less dense star clusters (lower binary production rate) or in more populous star clusters (higher binary disruption rate).
Fig. 4. Ratio of BHs per total initial number of stars in our models. Each vertical bar goes from the lowest to the highest number of BHs found within the realisations. 

Open with DEXTER 
Fig. 5. Mean time that a star, which ended up as a BH, spent in a binary (or multiple) system before becoming a BH. The vertical uncertainties represent one half of the output period. 

Open with DEXTER 
4.2. Analytic estimate in comparison to Nbody models
We separate the results into several plots (Figs. A.2, A.3 and A.4) depending on the initial virial radius of the model, which defines the value of the Plummer radius in Eq. (6) by . The upper limit of the predicted η_{BH} is calculated as if all the BHs were ejected directly from the centre of the cluster, that is, r = 0 pc. For the lower limit, we assumed that the ejections take place at the tidal radius determined from Eq. (4) with M_{C} = M_{tot} (other parameters are the same as in the tidal radius of the models in Sect. 3.1). If our simple analytic estimate is correct, we expect the retention fraction determined from the models to be within the shaded area between 0 and r_{t}. For a better understanding of the scale, in the plots, we also include a dotted line that corresponds to the escape radius of r = 10 pc (which is roughly the radius of an expanded star cluster) and a dashdotted line equal to the initial virial radius of each cluster, that is, r_{v} = 0.5, 1.0, and 2.0 pc.
In all the figures, we see a general agreement of the analytic estimates and the results from the Nbody simulations. Especially for the retention fraction in the halfmass radius. The best results are achieved for higher initial masses of the star clusters or wider models (larger r_{v}). This is partially expected from the nature of the escape velocity in Eq. (6), which is proportional to the square root of the mass of the cluster and inversely proportional to the square root of the initial Plummer radius.
Our analytic prediction does not consider more complicated effects, for example, dynamical friction, which can slow the BHs and increase their retention in the cluster (as we have shown in Fig. 2); the dynamics of core collapse (especially for small clusters <10k); or BHs bound in binary systems (mainly in the more compact clusters). Those could very easily shift the limits upwards and yield an even better agreement between η_{BH} and the results of the numerical models. In the case of higher σ_{BH}, where the analytic predictions are very strictly giving us η_{BH} = 0, we see more fluctuations in the value of η_{BH} taken from our models with r_{v} = 0.5 or 1.0 pc. Those are also the cases that need more investigation from the point of view of the dynamical effects.
Our analytic estimate has no upper limit for the total mass of the cluster. In Figs. A.2, A.3 and A.4 we are showing the results up to 5 × 10^{6} M_{⊙} where GCs would be. The trend set by our clusters provides a good indication of how the retention fraction in GCs could behave in those plots. Even further out to the right, outside of these plots, is the mass domain of the UCDs which is documented separately in Fig. A.5. We note that in the latter case, we use the birth radius according to Eq. (7) from Marks & Kroupa (2012), which is a reasonable assumption since the UCDs are expected to significantly expand to the observed presentday radii (Dabringhausen et al. 2008, 2010). The retention fraction of BHs in young UCDs should therefore be very high, as shown in Fig. A.5, but selfconsistent modelling will be needed to make improved estimates of how the retention fraction evolves as the UCDs expand as a result of the stellarevolution mass loss.
5. Conclusions
The BH retention fraction grows with the increasing initial mass of a star cluster. This is shown by direct Nbody simulations which corroborate our analytic estimates. Therefore, we conclude that it is possible to estimate the BH retention fraction using a simple analytic formula, according to which the number of escaping BHs is the number of BHs in the tail of the Maxwellian distribution above the escape velocity from the cluster; see Eq. (7). This implies that UCDs should have retained more than 80% of their BHs for σ_{BH} ≤ 190 km s^{−1}.
The estimate agrees with our Nbody results, especially for the retention fraction within the halfmass radius and for star clusters with an initial mass greater than 10^{4} M_{⊙} and an initial virial radius r_{v} ≥ 1 pc. In the cases where the analytic prediction does not follow the numerical results, the former serves as a good lower estimate for the retention fraction. When applying the standard velocity dispersion of SN kicks (σ_{BH} ≈ 190 km s^{−1}), only a few BHs remained bound to the modelled star clusters. Those BHs were either in binaries or they received only a very small kick by chance. With such a high σ_{BH}, the only places where BHs could be retained is in very massive GCs, UCDs, or an environment comparable to a nuclear cluster. Other possibilities for retaining BHs with high σ_{BH} are: (a) the kick velocity is to be scaled by a fraction of the envelope that falls back onto the BH, for example, from Eq. (3) Belczynski et al. (2002), or (b) they might implode, leaving BH remnants without any kick, that is, putting p_{env+,−} = 0 in Eq. (1). Due to those reasons, if the upper limit of the retention fraction of BHs is about 50% in GCs (Peuten et al. 2016; Baumgardt & Sollima 2017), then the velocity distribution of the kicks cannot be just one Maxwellian distribution with a velocity dispersion of σ_{BH} ≳ 190 Km s^{−1}. This suggests that, if σ_{BH} ≳ 50 Km s^{−1}, implosions or some bimodal kick velocity distribution could be valid, as implied for neutron stars by Verbunt et al. (2017).
Another result we find here is that compact clusters provide an environment in which particularly massive binaries form dynamically and then evolve by binary star evolution to enhance the number of SNe (exploding later than the original SN) providing additional delayed BH formation. Instead of the last SN exploding at ≈12 Myr, which is typical for the metallicity and mass range of our models, several BHs needed twice this time to get to the SN stage. This is evident especially in the smallest and most compact clusters, that is, 1k stars with r_{v} = 0.5 pc, although the tendency to later SN explosion is visible in all 1k models, and individual cases of more populous models. The stellar evolution in dynamically formed binary stars (the models here have no primordial binaries) affects mostly the evolution of highmass stars in small and compact clusters as it is virtually impossible to disrupt them in such systems. The other effect of living in a binary for most of the star’s life is the overproduction of BHs. An initially relatively lessmassive star in a binary system can accrete enough mass to explode as a SN and leave an additional BH that cannot be predicted from the IMF. We have detected this in the case of our 1k model with r_{v} = 0.5 pc. More populous and less dense clusters also produce dynamically formed binaries, but these do not evolve to produce additional late SNe and BHs.
Acknowledgments
This study was supported by the Charles University grants SVV260441 and GAUK186216. We thank Sverre Aarseth and Douglas Heggie for the discussion about nbody6 and Ondřej Chrenko for helpful comments. We also greatly appreciate access to the computing and storage facilities owned by parties and projects contributing to the National Grid Infrastructure MetaCentrum, provided under the programme “Projects of Large Research, Development, and Innovations Infrastructures” (CESNET LM2015042). We are grateful to the referee for useful comments.
References
 Aarseth, S. J. 2003, Gravitational NBody Simulations (Cambridge, UK: Cambridge University Press) [CrossRef] [MathSciNet] [Google Scholar]
 Aarseth, S. J., Henon, M., & Wielen, R. 1974, A&A, 37, 183 [NASA ADS] [Google Scholar]
 Banerjee, S. 2017, MNRAS, 467, 524 [NASA ADS] [Google Scholar]
 Banerjee, S. 2018, MNRAS, 473, 909 [NASA ADS] [CrossRef] [Google Scholar]
 Banerjee, S., Baumgardt, H., & Kroupa, P. 2010, MNRAS, 402, 371 [NASA ADS] [CrossRef] [Google Scholar]
 Baumgardt, H., & Sollima, A. 2017, MNRAS, 472, 744 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Kalogera, V., Rasio, F. A., et al. 2008, ApJS, 174, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Bulik, T., Fryer, C. L., et al. 2010, ApJ, 714, 1217 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Repetto, S., Holz, D. E., et al. 2016, ApJ, 819, 108 [NASA ADS] [CrossRef] [Google Scholar]
 Binney, J., & Tremaine, S. 1994, Galactic Dynamics (Princeton, NJ: Princeton University Press) [Google Scholar]
 BlandHawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Casertano, S., & Hut, P. 1985, ApJ, 298, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Dabringhausen, J., Hilker, M., & Kroupa, P. 2008, MNRAS, 386, 864 [NASA ADS] [CrossRef] [Google Scholar]
 Dabringhausen, J., Fellhauer, M., & Kroupa, P. 2010, MNRAS, 403, 1054 [NASA ADS] [CrossRef] [Google Scholar]
 De Donder, E., & Vanbeveren, D. 2003, New Astron., 8, 817 [NASA ADS] [CrossRef] [Google Scholar]
 Faber, S. M., & Gallagher, J. S. 1979, ARA&A, 17, 135 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Favata, M., Hughes, S. A., & Holz, D. E. 2004, ApJ, 607, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749 [Google Scholar]
 Hansen, B. M. S., & Phinney, E. S. 1997, MNRAS, 291, 569 [NASA ADS] [Google Scholar]
 Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543 [NASA ADS] [CrossRef] [Google Scholar]
 Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [NASA ADS] [CrossRef] [Google Scholar]
 Jeřábková, T., Kroupa, P., Dabringhausen, J., Hilker, M., & Bekki, K. 2017, A&A, 608, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jonker, P. G., & Nelemans, G. 2004, MNRAS, 354, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Kroupa, P. 2008, in The Cambridge NBody Lectures (Berlin: Springer Verlag), eds. S. J. Aarseth, C. A. Tout, & R. A. Mardling, Lect. Notes Phys., 760, 181 [NASA ADS] [CrossRef] [Google Scholar]
 Kroupa, P., Weidner, C., PflammAltenburg, J., et al. 2013, The Stellar and SubStellar Initial Mass Function of Simple and Composite Populations, eds. T. D. Oswalt, & G. Gilmore, 115 [Google Scholar]
 Lyne, A. G., & Lorimer, D. R. 1994, Nature, 369, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Mackey, A. D., Wilkinson, M. I., Davies, M. B., & Gilmore, G. F. 2007, MNRAS, 379, L40 [NASA ADS] [CrossRef] [Google Scholar]
 Mapelli, M., Ripamonti, E., Zampieri, L., & Colpi, M. 2011, MNRAS, 416, 1756 [NASA ADS] [CrossRef] [Google Scholar]
 Mapelli, M., Zampieri, L., Ripamonti, E., & Bressan, A. 2013, MNRAS, 429, 2298 [NASA ADS] [CrossRef] [Google Scholar]
 Marks, M., & Kroupa, P. 2012, A&A, 543, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Morscher, M., Umbreit, S., Farr, W. M., & Rasio, F. A. 2013, ApJ, 763, L15 [NASA ADS] [CrossRef] [Google Scholar]
 Nitadori, K., & Aarseth, S. J. 2012, MNRAS, 424, 545 [NASA ADS] [CrossRef] [Google Scholar]
 Peuten, M., Zocchi, A., Gieles, M., Gualandris, A., & HénaultBrunet, V. 2016, MNRAS, 462, 2333 [NASA ADS] [CrossRef] [Google Scholar]
 Plummer, H. C. 1911, MNRAS, 71, 460 [NASA ADS] [CrossRef] [Google Scholar]
 Podsiadlowski, P., Joss, P. C., & Hsu, J. J. L. 1992, ApJ, 391, 246 [NASA ADS] [CrossRef] [Google Scholar]
 Repetto, S., Davies, M. B., & Sigurdsson, S. 2012, MNRAS, 425, 2799 [NASA ADS] [CrossRef] [Google Scholar]
 Repetto, S., Igoshev, A. P., & Nelemans, G. 2017, MNRAS, 467, 298 [NASA ADS] [Google Scholar]
 Verbunt, F., Igoshev, A., & Cator, E. 2017, A&A, 608, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weatherford, N. C., Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2017, ApJ, submitted [Google Scholar]
 Zapartas, E., de Mink, S. E., Izzard, R. G., et al. 2017, A&A, 601, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ziosi, B. M., Mapelli, M., Branchesi, M., & Tormen, G. 2014, MNRAS, 441, 3703 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A
Additional tables and figures
BH retention fraction in our models.
Fig. A.1. Evolution of the retention fraction, η_{BH}, within the halfmass radius (orange) and the tidal radius (red). The black line shows the ratio of BHs that have formed with respect to the final number of BHs in the cluster. The curves for 1k and 10k models are averaged over 100 and 10 realisations, respectively. 

Open with DEXTER 
Fig. A.2. A comparison of the retention fraction from the Nbody models (with r_{v} = 1.0 pc) with the analytic estimate. The limits from our analytic estimate, Eq. (6), are represented by a shaded area delimiting the escape radius from 0 pc (solid line) to the tidal radius (dashed line) determined from Eq. (4). We also include the curve for an escape radius of 10 pc (dotted line) and the virial radius (dashdotted line). The squares and crosses represent the retention fraction from our numerical simulations in the halfmass and tidal radius, respectively. 

Open with DEXTER 
Fig. A.3. As in Fig. A.2 but for the models with the initial virial radius r_{v} = 0.5 pc. 

Open with DEXTER 
Fig. A.4. As in Fig. A.2 but for the models with the initial virial radius r_{v} = 2.0 pc. 

Open with DEXTER 
Fig. A.5. Analytic estimate of the retention fraction for a broader range of masses. Here, we use the relation from Marks & Kroupa (2012, Eq. (7)) to compute the appropriate radius of the system with a given mass. The escape velocity, Eq. (6), is calculated from 0 pc (solid line), the tidal radius (dashed line) determined from Eq. (4), and the halfmass radius of the model (dashdotted line). 

Open with DEXTER 
All Tables
All Figures
Fig. 1. Evolution of the halfmass radii (r_{h} ≈ 0.769 r_{v}, solid lines) and tidal radii (dotted lines) in our models. We compare models with the lowest and highest σ_{BH}. Models with different N are colourcoded as given by the key at the top of the figure. 

Open with DEXTER  
In the text 
Fig. 2. Average ratio between the number of BHs that received v_{kick} > v_{esc} during a SN explosion (N_{BH,kick}) and the number of BHs that are outside the tidal radius at 12 Myr (N_{BH,tide}). Each colour corresponds to a different kick velocity dispersion σ_{BH} listed above the plots. The error bars correspond to the Poisson uncertainties; see Eq. (9). 

Open with DEXTER  
In the text 
Fig. 3. Time of BH formation in the presented models. Circles correspond to the time of the last BH formation (i.e. f_{BH} = 100%), triangles represent the time when a fraction f_{BH}> 97% of BHs formed. For reference, the dotted line marks 11.7 Myr, as deduced from the Hurley et al. (2000) algorithm. 

Open with DEXTER  
In the text 
Fig. 4. Ratio of BHs per total initial number of stars in our models. Each vertical bar goes from the lowest to the highest number of BHs found within the realisations. 

Open with DEXTER  
In the text 
Fig. 5. Mean time that a star, which ended up as a BH, spent in a binary (or multiple) system before becoming a BH. The vertical uncertainties represent one half of the output period. 

Open with DEXTER  
In the text 
Fig. A.1. Evolution of the retention fraction, η_{BH}, within the halfmass radius (orange) and the tidal radius (red). The black line shows the ratio of BHs that have formed with respect to the final number of BHs in the cluster. The curves for 1k and 10k models are averaged over 100 and 10 realisations, respectively. 

Open with DEXTER  
In the text 
Fig. A.2. A comparison of the retention fraction from the Nbody models (with r_{v} = 1.0 pc) with the analytic estimate. The limits from our analytic estimate, Eq. (6), are represented by a shaded area delimiting the escape radius from 0 pc (solid line) to the tidal radius (dashed line) determined from Eq. (4). We also include the curve for an escape radius of 10 pc (dotted line) and the virial radius (dashdotted line). The squares and crosses represent the retention fraction from our numerical simulations in the halfmass and tidal radius, respectively. 

Open with DEXTER  
In the text 
Fig. A.3. As in Fig. A.2 but for the models with the initial virial radius r_{v} = 0.5 pc. 

Open with DEXTER  
In the text 
Fig. A.4. As in Fig. A.2 but for the models with the initial virial radius r_{v} = 2.0 pc. 

Open with DEXTER  
In the text 
Fig. A.5. Analytic estimate of the retention fraction for a broader range of masses. Here, we use the relation from Marks & Kroupa (2012, Eq. (7)) to compute the appropriate radius of the system with a given mass. The escape velocity, Eq. (6), is calculated from 0 pc (solid line), the tidal radius (dashed line) determined from Eq. (4), and the halfmass radius of the model (dashdotted line). 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.