MOCCA-Survey Database: Extra Galactic Globular Clusters. III. The population of black holes in Milky Way and Andromeda - like galaxies

Numerous observational studies have found black holes (BHs) and accreting BHs candidates in Galactic and extragalactic globular clusters (GCs) (Maccarone et al., 2007; Barnard & Kolb, 2009; Roberts et al., 2012; Miller-Jones et al., 2015; Minniti et al., 2015; Bahramian et al., 2017; Dage et al., 2018). Radial velocity measurements of a binary system in NGC 3201 provide the strongest proof that a BH exists in a Galactic GC (Giesers et al., 2018, 2019). Additionally, both electromagnetic emission and dynamical mass measurements from kinematic observations of extragalactic GCs indicate the presence of a significant fraction of unseen mass, possibly stellar-mass BHs in binary systems and intermediate-mass BH (IMBH) (Taylor et al., 2015; Dumont et al., 2022). Similarly, recent studies have also looked for indicators of the presence of BHs in GCs using numerical simulations of GC models containing sizable populations of BHs (Morscher et al., 2015; Arca-Sedda et al., 2016; Askar et al., 2018; Arca Sedda et al., 2019; Weatherford et al., 2019).

Generally, the most massive GCs should form up to several thousands of BHs in the first few Myr of cluster evolution. The natal kicks that these BHs experience upon birth and the cluster’s escape velocity have a substantial impact on how many of these BHs can be kept in the GCs, and it is crucially affected by the uncertain physics of stellar collapse (Belczynski et al., 2002, 2010; Fryer et al., 2012; Repetto et al., 2017; O’Shaughnessy et al., 2017). Retained BHs would segregate rapidly, populating the central regions of their host GC (Portegies Zwart & McMillan, 2000, 2002a; Fregeau et al., 2004; Freitag et al., 2006; Arca-Sedda et al., 2016). The dynamics in the cluster central region is thus dominated by stellar-mass BHs, and can lead to different outcomes. One possibility is that dynamical interactions among the most massive BHs lead to their ejection, freeing the cluster centre from their dark content. Another possibility is that, despite dynamics, BHs form a subsystem dominating the innermost cluster region, or they merge among themselves or with other stars to build-up an IMBH. IMBHs have masses in the range $10^{2}-10^{5}\,M_{\odot}$ and they are considered as the link between the stellar-mass and supermassive BHs (Barack et al., 2019). Recent numerical studies (Breen & Heggie, 2013a, b; Heggie & Giersz, 2014; Webb et al., 2018; Banerjee, 2018; Kremer et al., 2018, 2019) showed that the BHS is not entirely decoupled from the rest of the GC, and that the energy demands of the host GC would control the evolution of its BHS (Breen & Heggie, 2013a, b). The presence of a massive BHS can be responsible for the cluster dissolution, due to the interplay of the strong energy produced from the BHS and tidal stripping (Giersz et al., 2019).

High stellar densities are required to form massive BHs such as IMBH, with a possible scenario of IMBH formation being repeated collisions in the central regions of GCs (Portegies Zwart & McMillan, 2002b; Portegies Zwart et al., 2004; Portegies Zwart & McMillan, 2007; Giersz et al., 2015; Mapelli, 2016). Indeed, IMBH can be form through multiple stellar mergers in binary system (Di Carlo et al., 2021; González et al., 2021; Arca Sedda et al., 2021; Rizzuto et al., 2021; Rizzuto et al., 2022, see also Maliszewski et al. 2022). So far, there is still no conclusive evidence of IMBH presence in Galactic GCs, even though they are considered to potentially host an IMBH (Bash et al., 2008; Maccarone & Servillat, 2008; Lützgendorf et al., 2013; Lanzoni et al., 2013; Kamann et al., 2014; Askar et al., 2017; Arca Sedda et al., 2019; Hong et al., 2020; Arca Sedda et al., 2020).

Multiple three-body interactions can cause the formation of BH-BH binaries (BBHs), which serve as a power supply for the cluster core. The continuous interactions between the BBHs and the other objects in the GC, would harden the binaries, until they would be ejected from the cluster core or be merged, releasing gravitational waves (Portegies Zwart & McMillan, 2000; Banerjee et al., 2010; Downing et al., 2010; Wang et al., 2016; Askar et al., 2017). Similarly, stars that interact with retained BHs are forced into wider orbits, causing the GC to expand and this can postpone core-collapse (Merritt et al., 2004; Mackey et al., 2008; Gieles et al., 2010; Wang et al., 2016; Kremer et al., 2019). The presence of a BHS or of a IMBH in the central region of a GC would importantly shape the structure of the host GC (Mackey et al., 2007; Zocchi, 2015; Arca Sedda et al., 2018; Baumgardt et al., 2020).

In this work we extend the study of the GC populations for the MW and M31 exploited in our companion paper (Leveque et al., 2022b, herefter Paper II). In previous papers in this series, we set up the machinery that would be used to populate external galaxies with their GC populations by combining the results from the MOCCA-Survey Database I with the MASinGa semi-analytic tool. In this work we would like to test for the first time our machinery against the GCs properties and their BH content simulated in our models for both MW and M31 populations. In particular, we compare the orbital properties of the MW GC population with the observed properties from the Bajkova catalogue (Bajkova & Bobylev, 2021). Also, we aim to constrain the spatial distribution in the galactic halo of different GCs properties for different GC dynamical states comparing our results with previous studies (Lützgendorf et al., 2013; Askar et al., 2018; Arca Sedda et al., 2019; Weatherford et al., 2019). Then, we intend to determine the properties of the BBH mergers reported in our simulations and the inferred BBH merger rate (Banerjee, 2022; Mapelli et al., 2022). Finally, we show the properties of BBHs present at 12 Gyr that could potentially be observed, and the number of BH and BH binaries that have been transported to the nuclear star cluster (NSC) by infalling star clusters. In Appendix A we present the statistical tests of the studied populations.

In this section, we will summarize the most important ingredients of our machinery. More details about all the physical assumptions are properly described in section 2 of Paper II (Leveque et al., 2022b).

The MASinGa (Modelling Astrophysical Systems In GAlaxies) program has been used to model the GC populations (Arca-Sedda & Capuzzo-Dolcetta, 2014a; Belczynski et al., 2018; Leveque et al., 2022b; Arca-Sedda, 2022). For each GC in the population, MASinGa simulates the orbital evolution while taking into consideration the galactic tidal field and shocks, which contribute to the cluster disintegration, dynamical friction, which pulls the cluster toward the galactic center, and internal relaxation, which controls the cluster mass loss and expansion/contraction.

The number of GCs and distributions of GC masses and galactocentric positions generated by MASinGa have been used to choose appropriate models from the MOCCA-Survey Database I (Askar et al., 2017) to reproduce the initial distributions (the details are provided in Sec. 2.3 in Paper II). However, few steps have been taken in order to determine the galactocentric position of the MOCCA models in the gravitational potential of the studied galaxy. Indeed, the simple point-mass approximation for the Galactic potential was used to evolve the MOCCA models, with the central galaxy mass being contained inside the GC’s orbital radius. The initial Galactocentric distance of a model in MOCCA-Survey Database I is defined by its initial mass and its tidal radius. Only a finite set of such values were considered (see Table 1 in Askar et al., 2017), with no specific initial density profile being taken into account while modeling the MOCCA-Survey Database I models. A circular orbit at Galactocentric distances between 1 and 50 kpc were assumed, and the GC’s rotation velocity was set to to $220$ $\rm{km\,\,s^{-1}}$ for the whole range of galactocentric distances. The correct galactocentric distance for a circular orbit in an external galaxy, for a given tidal radius $r_{tidal}$ and GCs mass $M_{GC}$ can be determined knowing the tidal radius for the MOCCA model and the density distribution of the simulated galaxy. From the results presented in Cai et al. (2016), it is possible to determine the galactocentric radius of a circular orbit on which a GC will experience an equivalent mass loss if it were on an eccentric orbit. In particular, the apocenter distance $R_{apo}$ for the eccentric orbit can be determined as $R_{apo}=R_{C}\cdot(1-0.71\cdot E_{GC})^{-5/3}$ (Cai et al., 2016), with $R_{C}$ the galactocentric distance for a circular orbit, and $E_{GC}$ the initial orbital eccentricity, chosen from a thermal distribution. The initial galactocentric distance for a circular orbit $R_{C}$ is known for the MOCCA models. Finally, the initial galactocentric position has been selected within the orbit apsis. For each MOCCA model, a total of 900 representations of each MOCCA-Survey Database I model were generated, with different initial orbital eccentricity and galactocentric distances - the MOCCA-Library. In this way, the same model could be populated with different orbital parameters and in different galactocentric distance regions. Each simulated cluster from the MOCCA-Survey Databse I has a different internal dynamical evolution (such as mass loss, half-mass radius, etc.). Consequently, the dynamical evolution of MOCCA-Library models that represent different MOCCA models are diverse. Hence, models from the MOCCA-Library are defined as unique when they represent different MOCCA models (more details are provided in Sec. 2.3.1 in Paper II (Leveque et al., 2022b)).

Only MOCCA models that survived their internal dynamical evolution up to 12 Gyr were selected to represent the MASinGa GC models. Indeed, the initial galactocentric position for each MOCCA model was chosen within the same initial galactocentric radius bin of the representative MASinGa model. Instead, the initial mass was chosen randomly from the GC initial mass function (GCIMF) cumulative distribution (a power law $dN/dm=b\cdot m^{-\alpha}$ with a slope of $\alpha=2$ function has been used as GCIMF), with lower and upper limits of $2\times 10^{5}\,\,M_{\odot}$ and $1.1\times 10^{6}\,\,M_{\odot}$ respectively. This mass range was set in order to reproduce the expected initial mass range for MW and M31, being of $\sim 10^{5}-10^{7}$. Indeed, the observed masses for the survived GCs located within 17 kpc from the galactic center range between $10^{4}-2\times 10^{6}\,\,M_{\odot}$ and $5\times 10^{4}-3\times 10^{6}\,\,M_{\odot}$, for MW and M31 respectively¹¹1In this calculation, we did not take in consideration GCs with masses at 12 Gyr smaller than $10^{4}\,\,M_{\odot}$, as they would be closed to the cluster dissolution, and they would be hard to be observed in external galaxies.. According to Webb & Leigh (2015), the initial GC mass was $\sim 5$ times greater than the actual observed values (more details are provided in section 3.3 of Paper II). On the other hand, the maximum mass in the GCIMF was chosen according to the maximum initial mass in the MOCCA models, being of $1.1\times 10^{6}\,\,M_{\odot}$.²²2For an initial power-law GCIMF between $M_{low}=10^{3}-10^{4}\,M_{\odot}$ and $M_{up}=10^{7}\,M_{\odot}$, we would expect a total number of low-mass clusters ($M<10^{5}\,M_{\odot}$) being around 1500-200 for $M_{low}=10^{3}\,M_{\odot}$ and $10^{4}\,M_{\odot}$ respectively, and a total mass of the GCs population around $2-3\times 10^{8}\,M_{\odot}$. The GCIMF cumulative distribution has been normalized to the initial MASinGa models mass distribution. Each model representation was successively removed from the MOCCA-Library, in order to avoid multiple instances of the same model to populate each mass and galactocentric position bin. Our procedure would guarantee that the total number of GCs in the population, the galaxy density distribution, and the GC IMF distribution for the MOCCA population would reproduce the initial conditions in the MASinGa population. Also, this procedure guaranteed an 85% minimum coverage of not repeated unique MOCCA models for both MW and M31 populations. Finally, most of the selected MOCCA models has an initial pericenter distance within 4 kpc from the galactocentric centre, and a initial thermal distribution for the eccentricity. The models that has been delivered to the NSC have small pericenter distances ($<0.5$ kpc) and really eccentric orbits ($>0.8$). The galaxy density profile have been modelled with a Dehnen (1993) density profile family, of the form:

where $M_{g}$ is the galaxy total mass in $M_{\odot}$, $r_{g}$ the galaxy length scale in kpc and $\gamma$ the density profile slope. The Dehnen models parameters’ values have been determined fitting the observational rotational curve (data from (Eilers et al., 2019) and (Chemin et al., 2009), for MW and M31 respectively) to the Dehnen rotational curve. The initial number of GCs in our models has been set by the total GC mass (defined as a fraction of the total galaxy mass) and the GCIMF. Even though the initial total number of GCs found in our simulations is much smaller compared to the observed total number of GC in both MW and M31, we found a similar number of initial GCs in our simulations and the observed one when we limit to GCs located within 17 kpc from the galactic center. Also, we found that the number of surviving GCs in our models is similar to the observed number of GCs within the selected mass and distance range for both MW and M31. The initial conditions used in our simulations to reproduce the MW and M31 population are summarized in Table 1.

Successively, the MOCCA models selected with the prescription described above have been used to follow the internal and external dynamics evolution of the GC population. This allowed not only to follow the actual evolution for GC’s mass, and half-mass radius but also to follow the evolution of the compact objects present in the system, together with the compact object binary evolution and their survival. On the other hand, the formulae adopted in MASinGa have been used to determine the galactocentric distance and eccentricity evolution.

It is important to underline that the stellar evolution prescription adopted in MOCCA models is outdated. The BH masses prescription used in MOCCA models follows the Belczynski et al. (2002) mass fallback formulae. In particular, proposed rapid and delayed supernova mechanisms (Fryer et al., 2012) were not implemented in stellar evolution prescriptions used in MOCCA-Survey Database I models. For this reason, the final BH masses are smaller compared to the observed values from GW detections (Abbott et al., 2021, 2022) and updated stellar and binary evolution prescriptions for BH progenitors (Kamlah et al., 2022). Also, the sub-sample of MOCCA models used to populate the studied galaxies have metallcities $Z$ of 0.02, 0.006, 0.005, 0.001, and 0.0002. ($Z_{\odot}=0.02$).

As mentioned in Paper II and in Madrid et al. (2017), the MOCCA results were able to recreate the N-body simulations for galactocentric distances down to a few kpc. For this reason, the region between $2$ and $17$ kpc has been the focus of the post-processing investigation and statistical analysis of GC populations’ properties in Paper II. In a similar way, we also restricted our analysis in this work to GCs that were located in the same galactocentric zone.

For both MW and M31, 100 galaxy models were created, all formed 12 Gyr ago and evolved until the present day as in the prescription described in Paper II. To reduce statistical fluctuations and have a more robust statistical representation of the models, the average values obtained from the galaxy models and their GC population have been considered. The repeating of the same unique model might distort the structural GC parameter distribution, biassing the simulated distribution toward the attributes of the unique models that were randomly picked the most. To prevent such bias, when estimating the radial distribution of each property, only one unique model inside each radial bin was examined. The average value of each population’s measurements, as well as the standard deviation, have been calculated for each attribute.

While in Paper II we discussed and studied the global properties of the GC population reproduced by our machinery (such as mass spatial distribution, half-light radius distributions, etc.), in this work we focus on studying the properties of the simulated GC population and their BH content.

The kinematic comparison with the MW GC population in the Bajkova catalogue (Bajkova & Bobylev, 2021) shows that the GC population simulated in our models represent decently the observed kinematic properties. To model the Galactic potential, the MOCCA models were simulated using a point mass approximation and GCs were assumed to move on a circular orbit. From Fig. 3 it is possible to note that the orbital approximations deployed in the MOCCA models is in relatively good agreement with the observational data in the MW GC population, within the simulated mass range. This is a crucial and important result for our models: indeed, the GCs in our machinery were initially placed in a circular orbit around the external galaxy and then they were modelled in elliptical orbits, given the findings in Cai et al. (2016). Nonetheless, despite the limitations that this assumption would imply, the distribution of the observed GC mass and their circular orbits are in agreement with our simulations.

The distribution of the mean GC mass for different dynamical models is similar between the MW and M31 population. The BHS models tend to be more massive than the Standard and IMBH models, and the IMBH models being more massive than the Standard ones. These results would suggest that the most massive GCs might contain a BHS in their center. Similarly, some correlation between the BHS mass and the galactocentric distance might exist, with more massive BHS GCs seen at small galactocentric distances.

As already shown in Paper I, the BHS models are expected to have a larger half light radius since the central energy generation (controlled by the BHs) is much stronger compared to the other systems, implying a more expanded system. On the other hand, it is expected that the influence of the IMBH would change the central properties of the GC: due to the deeper central potential, the system is expected to be more concentrated, implying a smaller half light radius. This is indeed seen in our simulations for different galactocentric distances in the both MW and M31 populations: the mean half light radius of the BHS models is larger than the Standard and IMBH models, with IMBH models being more compact at larger galactocentric distances. In comparison with the results shown in Paper I, the results shown in this work also take into consideration the interaction between the GC and the host galaxy in the survival of the GC itself. This might be a further support to our machinery results and assumptions.

The mean total BH mass in the system is significantly larger for the IMBH model compared to the BHS and Standard ones. For these models, the total BH mass in the GC is defined by the mass of the IMBH. Indeed, the presence of the IMBH would imply a high density and short dynamical interaction time-scale that would drive out all the massive BHs from the system. Also, the mean total BH mass in the GC for the IMBH models is larger at a smaller galactocentric distance. This might imply a correlation between the formation of an IMBH in a GC and the galactocentric distance and the local galactic density. Instead, the number of BHs in the Standard models are expected to be small (if any), meaning that it is expected to have an almost null BH mass in the system.

Lützgendorf et al. (2013) examined for the existence of a possible IMBH at the center of the 14 GCs in their sample using observed surface brightness profiles and velocity dispersion profiles. Six of them have been proposed by Lützgendorf et al. (2013) to host an IMBH in their center (NGC 1904, NGC 5139, NGC 5286, NGC 6266, NGC 6388, and NGC 6715), and among these, 4 GCs have a galactocentric distance greater than 5 kpc. The number density of GCs decreases for galacocentric distances greater than about 5 kpc, with a fewer number of GCs observed at larger distances. For this reason, we have set a minimum galactocentric distance of 5 kpc for this comparison. Arca Sedda et al. (2019) classified Galactic GCs IMBH (or BHS) based on how many MOCCA models had a BHS or a IMBH, finding a total of 35 models harbouring an IMBH. 16 of the 35 IMBH reported Galactic GCs in Arca Sedda et al. (2019) are found a distances $>5$ kpc.

Askar et al. (2018) and Arca Sedda et al. (2019) reported the number of GCs candidates that would harbour a BHS in their center, both using the MOCCA-Survey Database I results to identify BHS models. In Askar et al. (2018), the authors chose models based on their central surface brightness and the observed current half-mass relaxation time, with a total of 28 Galactic BHS GCs. Arca Sedda et al. (2019) reported a total of 23 models harbouring a BHS. Instead, the authors of Weatherford et al. (2019) found that the mass segregation parameter $\Delta$, which was derived from the 2D-projected snapshots of the models published in the CMC Cluster Catalogue (Kremer et al., 2020) correlates with the number of BHs in the system. A total of 29 Galactic BHS GCs were found in their work (considering only GCs that retain a number of BH $N_{BH}>50$). More than half of the reported GCs in Askar et al. (2018), Arca Sedda et al. (2018) ($\sim 70\%$) and in Weatherford et al. (2019) are found at distance $>5$ kpc. From our models, it might be expected that the probability to discover an IMBH or BHS model is higher at larger galactocentric distances. Indeed, the total number of Standard models seems to be comparable with the IMBH and BHS ones in the outskirts of the galactic halo: the percentage over the whole population of GCs with galactocentric distance $>5$ kpc is $\sim 20\%$ for both IMBH and BHS models. The mean IMBH mass decreases with the galactocentric distances, suggesting a link for the IMBH formation to the galactic field and galactocentric position. However, this is not seen for the BHS models.

The simulated merger rate within a volume $V=1\,\rm{Gpc^{3}}$ obtained by our simulations depends on the galaxy density number assumed. Considering the galaxy density number from the Illustris simulations (Vogelsberger et al., 2014), we have a merger rate of $\sim 13-23\,\,\rm{yr^{-1}\,\,Gpc^{-3}}$, meanwhile considering the interpolated density of MWEGs from Abadie et al. (2010), we obtain a merger rate of $\sim 1.0-2.0\,\,\rm{yr^{-1}\,\,Gpc^{-3}}$. The results from the Illustris simulation could be interpreted as maximum expected merger rate. Indeed, the galaxy density number used in this work assumes that all galaxies (particularly dwarf ones) in the cosmological cube considered in our study are MW- and M31-like galaxies. On the other hand, the results from the interpolated density of MWEGs should be considered as a minimum expected merger rate, due to the set of assumptions that have been made. Indeed, the merger rates for both MW and M31 in our simulations have been determined for a sub-sample of the whole GC population. As reported in Paper II, we constrained our study for GCs with initial masses between $2\times 10^{5}$ and $1.1\times 10^{6}M_{\odot}$, and within 17 kpc from the galactic center. The most massive GCs are actually excluded in our study. These might be expected to be host to a large number of BBHs and in particular to BBHs mergers. Also, one additional explanation for the smaller reported value in our simulation is the prescription used in the MOCCA-Survey Database I. As said before, the BH masses obtained in the old prescription are smaller compared to the new updated ones, leading to a larger GW decay time in our simulations. Finally, in our study we did not consider all models that initially populated the studied galaxies. Indeed, we have considered only GCs that would survive up to 12 Gyr, with the dissolved GCs not taken into account. The latter would contribute importantly to the merger rate.

Nonetheless the assumption made in the determination of our merger rate, our values are comparable to the one reported in Mapelli et al. (2022), where a local merger rate $R\sim 4-8\,\,\rm{yr^{-1}\,\,Gpc^{-3}}$ have been reported. In their analysis, Mapelli et al. (2022) explored the cosmic evolution of isolated BBHs and dynamically formed BBHs in nuclear star clusters, GC and young star clusters. In particular, they studied the BBH merger rate considering two main different supernova channels (Fryer et al., 2012) (that are rapid and delayed models) in the BH formation, finding different results according to the chosen prescription: the merger rate in the delayed models are roughly $40-60\%$ smaller compared to the merger rate in the rapid one. The different merger rate can be explained by the different minimum BH mass in the two models ($5M_{\odot}$ for rapid and $3M_{\odot}$ for the delayed), leading to a larger GW decay time in the delayed model compared to the rapid ones.

Banerjee (2022) investigated the importance of binary evolution and cluster dynamics in producing merging BBHs over cosmic time. The author performed a population synthesis for the modelled universe, deploying direct N-body simulations to model the evolution of young massive clusters, and stellar evolutionary models for isolated binaries. The author’s estimates of the intrinsic BBH merger rate density and the cosmic evolution are consistent with the findings in the GWTC-2. In particular, the intrinsic merger rate for BBHs considering only young massive clusters determined in Banerjee (2022) is of the order $\sim 1\,\,\rm{yr^{-1}\,\,Gpc^{-3}}$ for redshift $z=0$. This value is comparable to our results, with most of the models computed in Banerjee (2022) used the rapid supernova prescription.

Finally, Askar et al. (2017) determined the local merger rate density using models from the MOCCA-Database I. Assuming a GC star formation rate from Katz & Ricotti (2013), the author determined local merger rate density as function of the star formation rate and the probability of forming a BBH per unit delay time. The authors find a merger rate $R=5.4\,\,\rm{yr^{-1}\,\,Gpc^{-3}}$. Our results are in relatively good agreement with the value reported in Askar et al. (2017). Also, our estimates place in the ballpark of Ligo-Virgo Collaboration predictions for BBH merger rate of $R=23.9\,\,\rm{yr^{-1}\,\,Gpc^{-3}}$ (Abbott et al., 2021), though being on the lower side. Our minimum values are smaller compared to the value reported in Askar et al. (2017). This supports the interpretation that the small value reported by our simulations is connected with the limitations imposed to the GC sub-population selected in our study.

The number of observable BBH merger rates in a 10 years span is quite small: our results would suggest that it would be unlikely to observe a BBH merger rate in both MW and M31 galaxies. In the future, results from the MOCCA-Survey Database II could be used to better constrain the BBH mass ratio and properties, taking advantages of the new stellar evolution features (as described in Kamlah et al. 2022) and the GC gas expulsion evolution included in the Monte Carlo method (as described in Leveque et al. 2022a).

The orbital properties of the BBHs that merged within 10 and 13 Gyr strongly depend on the formation channel of the binaries. Dynamical BBHs form due to close few-body encounters, and for this reason their orbits can have higher eccentricities and larger semi-major axis values compared to those BBHs that form through evolution of primordial binaries.

As shown in Fig. 11, the majority of the BBHs observable at 12 Gyr are located in the central galactocentric regions for both MW and M31. On the other hand, a constant mean number of BBHs per GC has been observed. This implies that the larger number of BBHs in the central region is correlated to the total number of GCs being more numerous in this region of the galactic halo, being the mean number of BBHs constant.

The radial distribution and the histograms at different galactocentric distances of semi-major axis, orbital eccentricity and mass ratio of observable BBHs shown in Fig. 12 and Fig. 13 could suggest stronger dynamical interactions involving BBHs in the central region of the galaxy halo. Indeed, the orbital eccentricity seems to be more thermal and the semi-major axis seems to be larger at smaller galactocentric distances (as shown in Fig. 14). Instead, because of larger and less dense clusters in the outskirt of the galactic halo, the number of interactions that would thermalize the binaries could be expected to be smaller. This would imply a more circular orbital and smaller semi-major axis for BBHs hosted in GCs at larger galactocentric distances. This is supported by the galactocentric distance distribution of dynamical models in the galactic halo. As shown in Fig. 4, the Standard models dominate the central region of the galactocentric halo. Given that the Standard models are more numerous compared to the BHS ones, it might be expected that the number of interactions that a BBHs undergo is larger in a Standard model than a BHS one.

We find around $1,000-3,000$ BHs are transported into the galactic NSC over a 12 Gyr time span. This might have interesting consequences on the overall population of BHs in a galactic nucleus. If NSCs form exclusively in-situ, the fraction of stars turning into a BH can be retrieved by the initial mass function, and correspond to $\sim 8\times 10^{-4}$ for a Kroupa (2001) mass function, thus a number of $N_{BH,in-situ}\sim 20,000$ for a MW-like NSC. If a fraction $f$ of NSC mass is contributed by star cluster dispersal, the actual amount becomes $N_{BH,real}\sim(1-f)\,\,N_{BH,in-situ}+N_{BH,delivered}$. In our models, we assumed $f\sim 0.1$, leading to a negligible difference between full in-situ formation and a “mixed” formation process. Combining our previous work and the present analysis, we can define a BH-transport efficiency as the ratio between the total number of delivered BHs and the total mass accreted into the galactic nucleus, this being around $\eta=2,000/(3.5\times 10^{6})\sim 5.7\times 10^{-4}\,\,M_{\odot}^{-1}$, and $N_{BH,delivered}=\eta\times M_{accreted}$. Thus, if a NSC is totally contributed by infalling clusters, we would expect a number of BHs $N_{BH,infall}=\eta\times M_{NSC}\sim 14,300$. This simplistic analysis provides us with a range of the total number of BHs that might be inhabiting the nuclear regions of MW and Andromeda, being this in the range $(1.4-2.2)\times 10^{4}$.

In this paper we expanded the study of the MW and M31 population simulated with the machinery introduced in Paper II, investigating the BH content of their GC population.

Summarizing our main results:

• •

  The kinematic properties of the MW population are in agreement with the observed ones - see Fig. 1, Fig. 2 and Fig. 3. The observed orbital properties distributions (that are, orbital eccentricity, pericenter distance, and the galactocentric distance for a circular orbit determined using the prescription in Cai et al. (2016)) in the MW population has been reproduced by our models, further confirming the reliability of our semi analytic procedure. This is a significant outcome for our machinery, since the GCs in external galaxy have been initially populated on a circular orbit, and subsequently they were modelled in elliptical orbits using the prescriptions in Cai et al. (2016).

• •

  The mean GCs mass and the mean half-light radius for models that would harbour a BHS is larger than Standard and IMBH models (see Fig. 4, Fig. 5 and Fig. 7). On the other hand, Standard models are more numerous in the central region of the galaxy, with the number of IMBH and BHS comparable to the Standard ones at larger galactocentric distances.

• •

  A maximum and a minimum value observable BBH merger rate has been determined, with value $\sim 20-35\,\,\rm{yr^{-1}\,\,Gpc^{-3}}$ using the galaxy number density from the Illustris-1 simulation (Vogelsberger et al., 2014), and of $\sim 1.0-2.0\,\,\rm{yr^{-1}\,\,Gpc^{-3}}$ using the extrapolated density of MWEGs from Abadie et al. (2010) respectively. The reported value are comparable to the value reported in previous works (Mapelli et al., 2022; Banerjee, 2022; Askar et al., 2017), although our results are on the lower-side with respect to other models (e.g. Rodriguez & Loeb (2018)). This could be due to the simplistic approach followed in calculating the merger rate. These differences can be explained by a smaller sub-sample of the whole GC population that have been considered and simulated in our machinery.

• •

  The signature of primordial or dynamically formed BBHs is imprinted in the orbital parameters for the merged binaries. Indeed, the dynamically formed binaries have greater mass ratios and more eccentric orbits than the primordial ones. Furthermore, it appears that the semi-major axis of the dynamically formed binaries is larger than that of the primordial binaries. Conversely, primordial mergers are characterised by a nearly flat eccentricity distribution and a mass-ratio clearly peaked around 0.5. The primordial merger eccentricity distribution subtly implies that dynamics might have aided the merging process, owing to the fact that isolated stellar evolution generally predicts nearly circular BBH mergers.

• •

  The observable BBHs at 12 Gyr that have been simulated in our models show different orbital properties for different galactocentric distances. The BBHs do show a larger (thermal) eccentricity, larger mass ratio ($>0.8$) and smaller semi-major axes ($<10^{2}R_{\odot}$) at smaller galactocentric distances. These spatial evolution can be explained by denser GCs in the central region, enhancing the number of strong interactions between the BBHs and the other stars in the GCs.

• •

  Most of the BH and BBH that are delivered to the NSC happens in the first $1-2$ Gyr of evolution, with a slow increase observed at later times. Also, the total number of BH binaries delivered to the NSC is $\sim 5\%$ of the total BH population delivered. A total of 1000-3000 BHs and 100-200 BBHs have transported into the nucleus over a time span of 12 Gyr. This implies a total number of BHs and BBHs lurking in NSCs being of $N_{BHs}=(1.4-2.2)\times 10^{4}$ and $N_{BBHs}=700-1,100$.

• •

  The efficiency of BH binary formation, or the total number per unit of mass, can be significantly enhanced due to dynamics in GCs. In fact, compared to isolated stellar/binary evolution in the Galactic field, GCs are about twice as efficient at producing binaries with at least one BH and even more effective at generating BBHs.

In the future we would like to extend the study of the MW and M31 population to further investigate the super massive BH and nuclear star cluster masses build up (Askar et al., 2022). Also, we intend to simulate with our machinery other galaxies and the galaxies in the local Universe. We hope to restrict and identify the observational properties, evolutionary paths, and compact object content of GCs (such as IMBH, BHS, BBHs, and X-ray binaries), and also study the gravitational microlensing phenomena in GCs. Our simulation results might be utilized to calculate the BBH merger rate in the local Universe, as well as the event rates of TDEs between the SMBH and infalling GCs.

We thank the anonymous reviewer for insightful comments that helped us clarify the presentation of the results in this paper. MG and AL were partially supported by the Polish National Science Center (NCN) through the grant UMO-2016/23/B/ST9/02732. AA acknowledges support from the Swedish Research Council through the grant 2017-04217. MAS acknowledges financial support from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 101025436 (project GRACE-BH, PI Manuel Arca Sedda). AO acknowledge support from the Polish National Science Center (NCN) grant Maestro (2018/30/A/ST9/00050). AO is also supported by the Foundation for Polish Science (FNP) and a scholarship of the Minister of Education and Science (Poland).

The data underlying this article will be shared on reasonable request to the corresponding author.