© The Authors 2026

1 Introduction

Black holes (BHs) can be broadly classified into two main categories based on their mass. Stellar-mass BHs lie at the low end of the mass scale, with masses ranging from about 3–5 M_⊙ up to roughly 100 M_⊙. These BHs form from the collapse of massive stars with initial masses greater than approximately 20 M_⊙ (see, e.g., Vink et al. 2001; Woosley et al. 2002; Fryer et al. 2012; Spera et al. 2015; Vink et al. 2021; Ugolini et al. 2025). At the opposite extreme are supermassive BHs (SMBHs), which have typical masses exceeding ∼10⁶ M_⊙. SMBHs are commonly found at the centers of galaxies and are believed to grow through complex evolutionary processes, starting from lighter BH seeds that increase in mass via various mechanisms (see, e.g., Volonteri et al. 2021; Izquierdo-Villalba et al. 2023). These seeds, which often represent intermediate states between the two BH mass categories, fall in the mass range between ∼10² and 10⁵ M_⊙ and belong to the elusive family of intermediate-mass black holes (IMBHs). IMBHs have been extensively studied through a variety of observational channels, including searches for potential radiative accretion signatures (Maccarone & Servillat 2008; Tremou et al. 2018; Lin et al. 2018; Tiengo et al. 2022), integral field spectroscopy and astrometric measurements of stellar motion (Gebhardt et al. 2005; Noyola et al. 2010; van der Marel & Anderson 2010; Lützgendorf et al. 2013a; Lanzoni et al. 2013; Kamann et al. 2014, 2016; Pechetti et al. 2022; Häberle et al. 2024), millisecond pulsar timing (D’Amico et al. 2002; Colpi et al. 2003; Ferraro et al. 2003; Kızıltan et al. 2017; Perera et al. 2017; Gieles et al. 2018; Bañares-Hernández et al. 2025), and gravitational-wave (GW) detections (Abbott et al. 2020; Abac et al. 2025b). Despite the increasing observational evidence, their existence is often difficult to confirm due to the absence of a definitive observational signature and the possibility that their effects can be mimicked by other astrophysical systems or processes. Investigating IMBH formation in different stellar environments, and its potential smoking guns, is therefore fundamental to understanding the nature of these objects and their possible links to stellar and SMBHs.

One promising route for IMBH formation is through dynamical assembly in a dense star cluster. Within this formation scenario, several mechanisms can play a role in the production and growth of the IMBH. Among them are the formation of a very massive star (VMS) from repeated stellar collisions that either collapses into an IMBH or accretes most of its mass on a stellar BH, repeated mergers among BHs and/or stars, or multiple tidal disruption events on stellar BHs (Miller & Hamilton 2002; Portegies Zwart & McMillan 2002; Giersz et al. 2015; Arca Sedda et al. 2023a; Rizzuto et al. 2023). Repeated stellar collisions may occur in light (≲10⁶ M_⊙) and compact (≲0.1–1 pc) star clusters, where core collapse sets in before massive stars evolve in compact objects. In these clusters, the density peak reached at collapse can trigger the onset of collisions among massive stars which can proceed in a runaway fashion, hence the term runaway stellar collisions (Portegies Zwart & McMillan 2002). The resulting VMS can reach a mass of ∼0.1-1% of the cluster’s initial mass, eventually collapsing into an IMBH. Conversely, hierarchical binary BH (BBH) mergers (Miller & Hamilton 2002) can trigger the build-up of an IMBH starting from a stellar BH progenitor or can yield the further growth of a seed formed from stellar collisions. This process is expected to take place preferentially in massive, compact clusters with high escape velocities (≳ a few hundred km/s).

Historically, these two fundamental formation channels have been extensively studied through direct N-body and Monte Carlo simulations (see, e.g., Giersz et al. 2015; Mapelli 2016; Kremer et al. 2020; Di Carlo et al. 2021; Rizzuto et al. 2021; Arca Sedda et al. 2023a; Vergara et al. 2024; Barber & Antonini 2024; González Prieto et al. 2024; Rantala et al. 2024; Vergara et al. 2026). However, the computational costs of the two techniques restrict the cluster parameter space that can be investigated (e.g., cluster structure, metallicity, formation history). An appealing alternative that became quite popular in the last few years is represented by semianalytic codes, which encode simplified recipes to model BH dynamics and significantly reduce the required simulation time (Antonini & Rasio 2016; Arca Sedda & Benacquista 2019; Antonini & Gieles 2020; Arca Sedda et al. 2020b; Mapelli et al. 2021; Arca Sedda et al. 2023b; Kritos et al. 2024; Torniamenti et al. 2024).

For this work, we employed the latest version of the semi-analytic population synthesis code B-POP (Arca Sedda & Benacquista 2019; Arca Sedda et al. 2020b, 2023b) to explore possible relations among different IMBH formation scenarios and the properties of their nursing environments. In Sect. 2 we present the main features of B-POP and the different models we investigated for the seeding and the clusters formation history. In Sect. 3, we briefly overview the efficiency of different IMBH formation and growth mechanisms based on the initial properties of the host clusters. Moreover, we present the results of our simulations, focusing on the emergence of IMBHs at low redshifts and their similarities to potential observed candidates. In Sect. 4 we discuss some interesting implications of our findings for different families of clusters. Finally, in Sect. 5 we summarize our main conclusions.

2 Methods

We studied the dynamical formation of IMBHs in star clusters under different seeding prescriptions and cluster formation histories using the semianalytic population synthesis code B-POP (Arca Sedda & Benacquista 2019; Arca Sedda et al. 2020b, 2023b). The code constructs a synthetic Universe filled with BBHs and characterized by different cosmic star formation rates for each cluster type, a redshift-dependent metallicity distribution, a reliable time evolution for the structural properties of star clusters, and a flexible environment where to test different assumptions for the BHs initial properties. A detailed description of the code can be found in the companion work (Arca Sedda et al., subm. to Phys. Rev. D), our previous works (Arca Sedda et al. 2023b) and in Appendix A.

2.1 The B-POP code

We used B-POP to simulate 10⁸ BBHs evenly split among young (YCs), globular (GCs), and nuclear star clusters (NSCs). In Table 1 we summarize the main parameters of our simulations, which we describe in more detail in the following along with the code’s main features.

2.1.1 Clusters

The initial mass and half-mass radius distributions of each cluster family were tuned on observations at low redshifts. YCs were modeled to match observations in the Milky Way and the Local Volume (Portegies Zwart et al. 2006; Gatto et al. 2021). While the half-mass radius of these clusters is well constrained, dynamical mass estimates are available for only a few systems. Therefore, B-POP samples the half-mass radii of YCs directly from the observed distribution, while their masses are drawn from the GCs mass distribution, reduced by 2.5 dex, in line with observations of YCs in the Milky Way and its satellites (Portegies Zwart et al. 2010). The masses and radii of GCs were extracted from the observed distributions of galactic GCs (Harris 1996; Baumgardt & Hilker 2018), while the properties of NSCs were tailored to match observations of these clusters in local galaxies (Georgiev et al. 2016).

In Fig. 1, we show the different cluster formation redshift distributions used in this work (see also Fig. 4 in our previous paper, Arca Sedda et al. 2023b). In all our simulations, we assumed that YCs followed the star formation history of their host galaxies, sampling their formation redshifts from the star formation rate (SFR) of Madau & Fragos (2017, hereafter MF17). For GCs and NSCs we considered two options. In the first, we adopted a flat redshift distribution in the range z ∈ [2, 8], as suggested by observations of extragalactic GCs by Katz & Ricotti (2013, hereafter KR13). In the second, we sampled the NSCs formation redshifts from the SFR of MF17 and the GCs formation redshifts from the distribution obtained in the cosmological simulations of El-Badry et al. (2019, hereafter EB18).

Regardless of the formation scenario, we assumed that the mean metallicity of stellar progenitors evolves with redshift as log⟨Z/Z_⊙⟩ ≃ 0.153 − 0.074 z^1.34 (MF17). For each value of Z, we assumed the metallicity distribution to be spread around the mean value following a log-normal distribution peaked at log⟨Z/Z_⊙⟩ and with dispersion σ_Z = 0.2 (Bavera et al. 2020; Santoliquido et al. 2022; Mapelli et al. 2022).

Fig. 1 Initial distribution of cluster formation redshifts. The different colors refer to different assumptions on the initial formation histories of GCs (solid lines) and NSCs (dashed lines). The initial redshift distribution of YCs is fixed in all models (gray solid line). The distributions are scaled such that the area subtended by the curve is 1.

2.1.2 Black holes

In B-POP, the masses of BHs formed from single stars or in stellar binaries are drawn from catalogs generated using external stellar evolution tools. To ensure continuity with our previous works (Arca Sedda et al. 2023b), we adopted catalogs generated with MOBSE (Giacobbo & Mapelli 2018, 2020) covering metallicities in the range Z ∈ [0.0002, 0.02]. In particular, BH masses were extracted from either (i) a simple single stellar BH mass (SSBM) spectrum, corresponding to BHs formed from single stars, or (ii) a mixed single BH mass (MSBM) spectrum, which represents the population of BHs formed from binary stellar systems (Arca Sedda et al. 2023b). In all our runs BHs had a 50% probability to have a mass sampled from the MSBM or SSBM spectrum. We note that this subtly implied the inclusion in our dynamical models of effects related to primordial binary evolution, for example star mergers or accretion episodes (Arca Sedda et al. 2023b), a feature that may have a significant impact on the emergence of single and double IMBHs with a mass <300 M_⊙ (see, e.g., Paiella et al. 2025). In both mass spectra the initial mass of a stellar-mass BH never exceeds the minimum IMBH mass, i.e., 100 M_⊙. Nonetheless, additional features implemented in B-POP allow the user to assign BH masses in and beyond the pair instability supernova (PISN) range (Woosley 2017), for example assuming that their progenitors formed from one or multiple stellar collisions. We exploited this feature to explore the impact of different prescriptions for the early seeding of IMBHs, as detailed in Sects. 2.2.2, 2.2.3, and 2.3.

In our runs a BH (IMBH) can grow dynamically via repeated BBH mergers. After each merger the retention of the remnant is ensured only if its relativistic kick is lower than the escape velocity of the environment (see more on this in Sect. 2.2.4). This kick, which arises from asymmetries in GW emission prior to the merger, is primarily influenced by the mass ratio (q = m₁/m₂ < 1) of the BBH and the magnitude and relative alignment of the initial spins (Campanelli et al. 2007; Lousto & Zlochower 2008; Lousto et al. 2012). In our simulations we always assumed a BH (IMBH) to pair up with BHs which had not undergone any previous merger. Allowing higher-generation BH secondaries can, in principle, enhance hierarchical IMBH growth. However, such mergers are rare in relatively low-mass systems and are typically limited to the third generation by cluster expansion (Torniamenti et al. 2024). In general, multiple IMBHs are unlikely to form and grow within the same environment, except for clusters with a large fraction of primordial binaries, with extreme densities, or undergoing a hierarchical formation (see, e.g., Gürkan et al. 2006; Giersz et al. 2015; Mapelli 2016; Kovetz et al. 2018; Rodriguez et al. 2019; Anagnostou et al. 2020; González et al. 2021; Arca Sedda et al. 2023a; Rantala & Naab 2025). A prescription for the pairing up of BHs from different generation is under construction for the future public release of the code (Ugolini et al., in prep.).

Concerning the clusters’ BHs initial spins, we always assumed them to be distributed isotropically and with magnitudes following a Maxwellian distribution with dispersion σ_χ = 0.2. This distribution is consistent with the inferred spin distribution of merging BBHs observed by the LVK collaboration up to the O3 run (Abbott et al. 2023).

2.2 IMBH formation channels

In the following we briefly describe the processes regulating IMBH seeding and growth considered in this work, namely the stellar collision scenario and the hierarchical BBH merger scenario. We note that the two processes are not mutually exclusive; for example, an IMBH formed from the collision of a few massive stars can indeed grow further by merging with stellar BHs. For the sake of conciseness, we provide a more in-depth discussion about these mechanisms in Appendix A.2.

2.2.1 Cluster relaxation and core-collapse

Any gravitationally bound system undergoes dynamical relaxation driven by the continuous dynamical interactions among its constituents. In a star cluster, this process occurs over a timescale (Spitzer 1969; Binney & Tremaine 2008) $$t_{\text{rel}}=4.2\text{Gyr}\left(\frac{15}{\ln {\text{Λ}}_{\text{c}}}\right)\left(\frac{R_{\text{h},\text{cl}}}{4\text{pc}}\right)\sqrt{\frac{M_{\text{cl}}}{{10}^7{\text{M}}_{\odot }}},$$(1)

called half-mass relaxation time, which depends on the cluster mass (M_cl) and the radius within which half of the cluster mass is contained, i.e., the half-mass radius (R_h,cl). Relaxation causes light stars to diffuse outward, carrying kinetic energy away from the core. This energy loss makes the core contract, which in turn increases the diffusion, creating a loop that leads the core to collapse. This phenomenon is commonly referred to as gravothermal instability (Spitzer 1969).

The contraction is finally reversed by the formation of hard binaries, that transfer kinetic energy to other core members through dynamical interactions and leads the core to bounce back. The core collapse process is, therefore, fundamentally governed by dynamical interactions between light and massive stars. These interactions also promote the segregation of the most massive stars in the cluster center on a timescale (Spitzer & Hart 1971; Portegies Zwart et al. 2004) $$t_{\text{s}}(m)=\frac{{\overline{m}}_{\ast }}{m}\frac{0.138N_{\ast }}{\log \left({\text{Λ}}_{\text{c}}\right)}{\left(\frac{R_{\text{h},\text{cl}}^3}{G{M}_{\text{cl}}}\right)}^{1/2},$$(2)

where m is the mass of the sinking star, ${\overline{m}}_{*}$ is the average mass of a star in the cluster and $N_{\ast }\sim M_{cl}/{\overline{m}}_{\ast }$ is the average total number of stars assuming an initial cluster mass M_cl. For our simulations we fix ${\overline{m}}_{\ast }=0.7{\text{M}}_{\odot }$ and assume that the core collapse happens on a timescale equal to the segregation time of the most massive stars in the cluster population, i.e., for t_cc ∼ t_s(150 M_⊙). Moreover, we define a lower bound for the core collapse timescales of t_cc ∼ 0.2 t_rel (see also Portegies Zwart & McMillan 2002).

2.2.2 Runaway stellar collisions and IMBH seeds

If the core-collapse happens on timescales shorter than ∼3–5 Myr, which is the typical lifetime of BH massive stellar progenitors, it can trigger a chain of runaway stellar collisions and the formation of a VMS (Portegies Zwart & McMillan 2002) that may collapse into an IMBH.

The maximum mass of the VMS can be conveniently rewritten as a function of the initial cluster mass, $$m_{\text{VMS}}=m_{\text{seed}}+800{\text{M}}_{\odot }\left(\frac{f_{\text{c}}}{0.2}\right)\left(\frac{M_{\text{cl},0}}{{10}^5{\text{M}}_{\odot }}\right)\left(\frac{\ln {\text{Λ}}_{\text{c}}}{10}\right),$$(3)

with m_seed and f_c being the initial mass of the stellar seed and the effective fraction of colliding binaries in the cluster.

In general, Eq. (3) neglects the impact of stellar evolution processes which can affect both the amount of mass effectively accreted onto the growing VMS and the effect of metallicity on mass loss. For instance, strong stellar winds, typical of stellar populations with a metallicity Z ≳ 0.001, can effectively quench the VMS growth (e.g., Spera & Mapelli 2017; Iorio et al. 2023). Nonetheless, mass loss is expected to be nearly negligible for metal-poor stellar populations (see for instance Fig. 2 of Mapelli 2018). The final mass of the IMBH seeded by the VMS is also highly uncertain, although recent stellar evolution tracks for stars with masses up to 500–2000 M_⊙ suggest that a relatively small fraction (<0.2–0.3) of the VMS mass is lost during the collapse to an IMBH if the metallicity is Z < 0.001 (see Fig. 4 in Costa et al. 2025). We also note that the ultimate fate of a star resulting from multiple stellar collisions may differ from that of a similar star born in isolation due to the different chemical stratification of the collided star and the additional rotation it may acquire during the collisions (Glebbeek et al. 2009; Costa et al. 2022; Ballone et al. 2023; Vergara et al. 2026).

2.2.3 Mild seeding scenario: Stellar collisions in dense star clusters

While the runaway scenario requires extreme conditions, in clusters with longer core collapse timescales numerical simulations have shown that several stellar collisions, often triggered by strong dynamical encounters among the population of primordial binaries, can lead to the rapid build-up of IMBH seeds with masses in the range 10² up to a few 10³ M_⊙(e.g., Mapelli 2016; Kremer et al. 2020; Di Carlo et al. 2021; Arca Sedda et al. 2023b; Rantala et al. 2024; Rantala & Naab 2025). Hence, we consider an additional scenario in which we assume a 20 % probability that a cluster denser than ρ_cl,0 > 10⁵ M_⊙ pc⁻³ and with metallicity Z < 0.001 forms a VMS, as recently suggested by Arca Sedda et al. (2023a) and Rantala et al. (2024); Rantala & Naab (2025). Under this scenario, we assign the mass of the VMS randomly between m_VMS,min and m_VMS,max, defined as $$m_{\text{VMS,min}}=\min \left(150,{10}^{-0.23}\times {\left(M_{\text{cl},0}/{\text{M}}_{\odot }\right)}^{0.53}\right){\text{M}}_{\odot },$$(4) $$m_{\text{VMS,max}}=\min \left(2\times {10}^4,0.02\times \left(M_{\text{cl},0}/{\text{M}}_{\odot }\right)\right){\text{M}}_{\odot }.$$(5)

The lower bound on the minimum VMS mass matches the values in Table 1 of Rantala et al. (2024) which are derived from the cluster mass – maximum stellar mass relation in Yan et al. (2023). The upper bound on the maximum VMS mass at 2 × 10⁴ M_⊙ serves as a conservative upper bound, given the uncertainties in mass loss and partial disruption of the star during each collision.

Although, in principle, this seeding mechanism can act in concert with the runaway one (Sect. 2.2.2), their development intrinsically depends on the cluster environment. In clusters that are prone to runaway collisions, the timescale for VMS build-up and collapse (∼2–5 Myr) is short enough to prevent long-term stellar collisions and to suggest that the VMS final product can grow by accreting stellar material in a tidal disruption event or by merging with other BHs. Clusters that may sustain a mild seeding, instead, are unlikely to undergo a runaway phase because stellar BHs, formed before the core-collapse, will likely settle in the centre and heat-up the core, preventing star-star mergers. For these reasons, we considered the two processes separately, as described in Sect. 2.3.

2.2.4 IMBH assembly via hierarchical binary black hole mergers

Aside from stellar collisions, an IMBH can also grow in a cluster via hierarchical BBH mergers, as summarized below.

Massive objects in a stellar population drift to the center of the parent cluster over a segregation timescale (Eq. (2)), with subsequent three-body (Lee 1995) and binary – single interactions (Goodman & Hut 1993), promoting the formation of hard binaries (Heggie 1975). Though these interactions get progressively rarer, they also become more violent and can impart a recoil to the BBH that may exceed the cluster escape velocity, leading to the premature ejection of the binary. However, if the hardening rate due to GW emission becomes dominant over the interaction rate before ejection, the BBH will merge inside the cluster. In such a case, the remnant is either expelled owing to the post-merger relativistic kick, or remains in the cluster possibly participating in another merger event. In either case, the BBH will undergo the merger only if the sum of cluster formation time, the time needed for the binary to form and harden, and the final merging time does not exceed a Hubble time. A BH typically needs to undergo a handful of mergers to reach the IMBH regime, an outcome possible only, in general, if the cluster’s escape velocity is sufficiently high to prevent its ejection.

2.3 Models descriptions

We explored three models based on different prescriptions for the IMBH seeding or the formation history of star clusters. The models are summarized in the following and in Table 2:

• Model A: all clusters with a core collapse time shorter than 5 Myr and metallicity Z < 0.001 are assumed to trigger the formation of a VMS via runaway stellar collisions, followed by the seeding of an IMBH, as described in Sect. 2.2.2. We also invoke for the seeding a minimum cluster mass of 2 × 10⁴ M_⊙ in order to have, on average, at least one seed massive star in the cluster (see Appendix B). The VMS mass (m_VMS) is computed according to Eq. (3) assuming an initial star mass m_seed = 150 M_⊙ and an effective fraction of colliding binaries f_c = 0.2. The resulting IMBH is assumed to retain 80% of the VMS mass, consistently with the findings of Costa et al. (2025) for these metallicities. We then follow the dynamical processes that may determine the pairing and merger of the IMBH with other stellar-mass BHs. If core collapse happens after 5 Myr, the runaway process is suppressed, and the IMBH formation can occur only through hierarchical BBH mergers (Sect. 2.2.4).

  In this model, we assume YCs to follow the SFR from MF17, while for both GCs and NSCs we adopt the redshift distribution of KR13. The latter assumption supports the scenario in which NSCs share their formation history with GCs, as might occur if NSCs form through the (possibly rapid) merger of GCs migrated towards the galactic center via dynamical friction, a scenario commonly referred to as the migration scenario (Tremaine et al. 1975; Capuzzo-Dolcetta 1993; Neumayer et al. 2020, and references therein).

• Model A′: this model adopts the same assumptions of Model A, except for the formation histories of GCs and NSCs. As with YCs, we assume NSCs to follow the SFR from MF17. This choice reflects a scenario in which NSCs form directly through star formation and gas infall occurring near the galactic center, a process commonly referred to as the in situ formation scenario (Loose et al. 1982). In addition, we sample the GCs formation redshift according to a Gaussian distribution with mean redshift µ_z = 3.2 and σ_z = 0.2 (Mapelli et al. 2022). This choice reflects theoretical models aimed at tracking the formation of Galactic GCs (see, e.g., EB18).

• Model B: this model is based on the same set of assumptions of Model A regarding the cluster formation history, but assumes a mild seeding scenario in which stellar collisions can lead to the formation of a VMS with a 20% probability in clusters that are both metal-poor, Z < 0.001, and very dense, ρ_cl,0 > 10⁵ M_⊙ pc⁻³ (see also Sect. 2.2.3). As for the runaway seeding, we fix a minimum cluster mass (5 × 10³ M_⊙) to ensure the presence of massive stellar seeds for the collisions (see Appendix B). The mass of the VMS is sampled as in Eq. (5) while we assume the IMBH to be 80% the mass of its progenitor. As in Model A simulation, if the cluster does not meet the seeding conditions we simply follow the dynamical evolution of a stellar-mass BH which may grow into an IMBH via hierarchical BBH mergers.

3 Results

Hereafter, we refer to IMBHs formed from a stellar BH that underwent a chain of BBH mergers as hierarchical IMBHs, while we refer to IMBHs formed through direct collapse of a stellar collision product as IMBH seeds or seeded IMBHs. We note that IMBH seeds can also grow by further merging with stellar BHs.

3.1 Theoretical expectations

In Fig. 2 (left panel) we present the input distribution of initial masses and half-mass radii of YCs, GCs, and NSCs used in our simulations.

The runaway stellar collisions seeding, which we assume in Model A and A′, comprises all the environments falling on the left side of the dashed-dotted contour at t_cc = 5 Myr. As shown in the plot, this seeding mechanism typically favors YCs as well as very light GCs and NSCs. Approximately 11% of YCs, 39% of GCs, and 20% of NSCs in B-POP fall within this runaway seeding region.

Conversely, the mild collision scenario, which we assume in Model B, selects the environments which fall above the gray dotted density contour at (half-mass) density ρ_cl,0 = 10⁵ M_⊙ pc⁻³. This seeding condition generally favors NSCs, which present the densest environments among the different cluster families. In this scenario we find that ∼19%, 9%, and 32% of YCs, GCs, and NSCs mass-radius distributions fall in the portion of the parameter space where ρ_cl,0 > 10⁵ M_⊙ pc⁻³. Assuming an additional 20 % seeding probability for these clusters (see more in Sect. 2.2.3) leads to ∼4%, 2%, and 6% of YCs, GCs, and NSCs being potential IMBH seeds hosts in Model B.

It should be noted, however, that the actual fraction of clusters producing IMBH seeds in both models depends critically on the metallicity, as we assume that the VMS progenitors are assembled only in clusters with Z < 0.001. The impact of metallicity is strictly linked to the clusters formation history and can significantly decrease the effective fraction of seeds hosts as we discuss in more details in Sect. 3.2.

The efficiency of hierarchical BBH mergers in assembling an IMBH is more complicated to define since many processes can actually lead to the ejection of the original BH during its growth. As a reference, in Fig. 2 we present, in different colors, the regions of the parameter space where a BBH (i) is ejected prior to merger due to a strong Newtonian recoil, (ii) merges and produces a BH which is then ejected by its relativistic kick, or (iii) the BBH is retained in the cluster but fails to merge within a Hubble time (see Appendix A.2 for further details). The first two regions are drawn comparing the magnitude of the recoils and kicks to the initial escape velocity of the host clusters. To evaluate the GW kicks we average over all possible spin configurations assuming them to be isotropically distributed in space, and sampling the spins magnitudes from our fiducial Maxwellian distribution (see more in Sect. 2). The last region is evaluated computing the average delay time between the BBH formation and merger in the cluster, and comparing it to a Hubble time. We adopt a 30 M_⊙ − 30 M_⊙ BBH as our reference system in all clusters, except when the cluster lies in the regime of runaway stellar collisions (i.e., t_cc < 5 Myr). In that situation, we instead consider an IMBH-BH binary, with the primary mass set to m₁ = 0.8 m_VMS (Eq. (3)) and the secondary mass to m₂ = 30 M_⊙. For the stellar-mass BBH case, only the most massive NSCs lie outside the regions where the binary is either ejected or prevented from merging. In contrast, when an IMBH-BH binary forms from a seeded IMBH, the larger inertia of the system significantly reduces the probability of ejection, making retention within the cluster more likely even for YCs and GCs. This picture becomes, however, far more complicated when accounting for the temporal evolution of the host cluster which can strongly impact the IMBH formation and growth through the hierarchical channel. For instance, the cluster mass loss and expansion will drive clusters to wander downward and leftward in the M_cl,0–R_h,0 plane shown in Fig. 2, reducing the BH interaction rate and increasing the probability of BH ejection. This is especially the case of clusters with relaxation times comparable to the dynamical timescales over which the BBHs evolve. We explore all these aspects in more detail in the next section. Finally, in Fig. 2 we also provide a simple sketch of the processes contributing to the seeding and growth of an IMBH in a star cluster as well as the ejection mechanisms mentioned above.

Fig. 2 Left panel: initial cluster mass and half-mass radius for YCs (bottom left area), GCs (central area), and NSCs (upper right area). The cluster distributions are cut at the 99% contour. Clusters with tcc < 5 Myr lie below the dot-dashed black line; those with central density ρcl,0 > 105 M⊙ pc−3 are above the shaded dotted gray line. We also identify regions of the parameter space in which (i) a BBH binary is ejected before merging due to strong Newtonian recoils (in pink), (ii) the merger remnant is typically kicked out of the cluster due to its relativistic kick (in purple), (iii) the BBH merger time exceeds a Hubble time (in yellow). Right panel: schematic overview of a BH dynamical growth in a star cluster. The colored boxes refer to the regions highlighted in the left plot.

3.2 Demographics

We find that only a small fraction of the 10⁸ BH realizations produces an IMBH: about ∼0.6% in Model A, ∼0.4% in Model A′, and ∼0.1% in Model B. In all cases, the overwhelming majority of IMBHs originates from stellar collision seeds. Specifically, these IMBHs represent ∼97–98% of the whole IMBH population in Models A and A′, and ∼88% in Model B. Hence, in all our simulations stellar collisions appear to be the primary drivers of IMBH formation in star clusters. The overall IMBH fractions remain modest mainly due to the imposed constraint on the metallicity threshold that we require for the formation of a VMS via stellar collisions. It is also important to note that, in our simulations, YCs, GCs, and NSCs are assumed to be equally abundant. Therefore, these fractions do not reflect the true relative proportions of different cluster families in the Universe, which remain highly uncertain. For reference, the Milky Way hosts one NSC, ∼200 GCs (Harris 1996; Minniti et al. 2017), and possibly up to ∼10³ massive YCs clusters (Larsen 2009).

To quantify the IMBH formation efficiency within each cluster family, we report in Table 3 the final IMBH populations, distinguishing between their formation channel (hierarchical or seeded) and the host cluster type (YCs, GCs, NSCs). The percentages shown for each mechanism are computed relative to the total number of clusters belonging to the corresponding family. From Table 3 it is clear that seeding via stellar collisions is consistently more efficient than hierarchical growth, by at least one order of magnitude, regardless of the cluster type. The environments in which this mechanism is more frequent, however, depend strongly on both the seeding conditions (Model A vs. Model B) and the assumed cluster formation histories (Model A vs. Model A′). For instance, in Model A, GCs are more likely to be hosting a seed thanks to their shorter core collapse times, whereas in Model B the seeding is restricted to very dense clusters, hence IMBH seeds are formed more frequently in NSCs. A direct comparison between Models A and A′ further illustrates the impact of formation history. In Model A′, only about 0.05% of GCs produce an IMBH seed, compared to ∼0.8% of NSCs. This trend is reversed in Model A, where GCs are nearly twice as likely as NSCs to host a seed. The difference arises because, in Model A′, GCs statistically form at lower redshifts than in Model A, and thus tend to have higher metallicities, often exceeding the threshold for runaway or mild stellar collisions. While NSCs also tend to form later in Model A′, their redshift distribution extends up to z ≳ 10, broader than that of GCs. This wider redshift range partly offsets the metallicity constraint, allowing a nonnegligible fraction of NSCs to still meet the seeding conditions.

Across all three models, hierarchical IMBH formation is efficient only in NSCs, while YCs and GCs produce virtually no IMBHs through this channel (see Table 3). This is because only NSCs reach escape velocities high enough to retain stellar-mass BHs against the typical dynamical recoils and relativistic kicks they experience during their growth (see Sect. 3.3.1).

The differences observed for the IMBH formation efficiency across formation mechanisms and cluster environments also impact the final population of mergers involving IMBHs observed in our simulations. We discuss the properties these mergers in Appendix C.

3.3 IMBHs at low redshifts: Models and observations

3.3.1 Retention efficiency and ejection mechanisms

At redshift z = 0, only around 33% (Model A), 39% (Model A′), and 61% (Model B) of all IMBHs forming in the simulations are retained in their clusters. In general, the IMBHs in our simulations can be ejected by means of three different processes: (i) Newtonian recoils, which expel an IMBH-BH binary as a result of strong binary-single interaction with the other BHs in the cluster, (ii) relativistic kicks, which eject an IMBH remnant after the binary merged, and (iii) the complete evaporation of the IMBH host cluster. The latter process mainly affects light clusters with short relaxation times which completely disrupt within a Hubble time, releasing their IMBHs in the field (see also Sect. 3.3.2).

In Table 4, we summarize the final fate of IMBHs in our simulations, categorized by formation mechanism, cluster family, and ejection process. We also report the fraction of IMBHs retained at z = 0, with all percentages normalized to the total number of IMBHs formed in each cluster type and formation scenario. We focus on the comparison between Model A and B, as Model A′ exhibits trends similar to Model A. While Models A and A′ differ in the total number of IMBH seeds due to the different choice of initial metallicity (redshift) distributions (Sect. 3.2), the relative fractions of IMBHs that are ejected or retained remain essentially unchanged. This is because these fractions are primarily set by the initial BH and IMBH seed properties and by the escape velocities and dynamical evolution of their host clusters, which are identical in both models. In addition, for hierarchical IMBHs we only refer to NSCs, since the fractions of IMBHs formed in YCs and GCs through this channel are always smaller than 0.01% (see Table 3).

Model B exhibits a substantially higher fraction of IMBH seeds retained at z = 0 compared to Model A. This arises from the higher-density threshold imposed for IMBH seeding in Model B, which preferentially selects host clusters with larger escape velocities than those in Model A. In both models, IMBHs are primarily ejected through post-merger relativistic kicks. This is particularly evident for hierarchical IMBHs, which are retained at z = 0 in less than 15% of the times. In addition to dynamical ejection mechanisms, such as kicks or recoils, we find that the premature evaporation of YCs quenches the growth of IMBH seeds in 11–17% of cases, depending on the model. In Sect. 4.1 we discuss possible observational implications for this population of orphan IMBHs.

3.3.2 Comparison to low-redshift IMBH candidates

In Fig. 3, we compare the IMBH population retained at z = 0 in our simulations with several observations of IMBH candidates in local GCs (Lützgendorf et al. 2013b; Kızıltan et al. 2017; Häberle et al. 2024), in extragalactic clusters from potential tidal distruption events (Lin et al. 2018), in galactic nuclei (Graham & Spitler 2009; Neumayer & Walcher 2012; Kormendy & Ho 2013), in dwarf galaxies, ultra-compact dwarf galaxies (UCDs) and low-mass active galactic nuclei (AGNs; Reines & Volonteri 2015; Chilingarian et al. 2018; Nguyen et al. 2018). In the following we restrict the analysis to Model A and Model B, as Model A and Model A′ produce results broadly consistent with each other. Model A is shown in the left panel of Fig. 3, while Model B is displayed in the right panel.

In Model A, IMBH seeds can form only in clusters with a core collapse time t_cc ≲ 5 Myr, a condition that translates into an upper limit on the clusters relaxation times and, consequently, their lifetime. This effect is evident in Fig. 3 for YCs and low-mass GCs, which in our simulations either evaporate before z = 0 (see also Table 4) or undergo substantial mass loss, retaining IMBHs with typical final mass ratios M_IMBH/M_cl ≳ 0.1. Conversely, hierarchical IMBHs are retained only in clusters with final masses ≳10⁷ M_⊙ and exhibit typical mass ratios M_IMBH/M_cl ≲ 10⁻², filling a region of the cluster mass–IMBH mass plane partly populated by observed NSCs. The growth proceeds uninterrupted in the densest clusters with final masses > 10⁸ M_⊙, generally leading to the exhaustion of the BH reservoir. These clusters show notable similarities with more complex stellar systems, such as UCDs and light low-mass AGNs. In Model A we observe a distinct gap around M_cl ∼ 10⁷, M_⊙ between the seeded and hierarchical populations, which lie to the left and right of the gap, respectively. The gap emerges from the combined effect of the seeding conditions and the inefficiency of hierarchical BBH mergers in clusters falling in that mass range. For completeness, in Appendix D we also discuss how the gap may be affected assuming GCs to be more massive and compact at birth. The IMBH mass – host mass relation changes significantly in Model B. In this model the seeded and hierarchical retained populations fill homogeneously the cluster mass range ∼10³–10⁹ M_⊙ and clearly overlap with the masses of putative IMBHs observed in galactic and extragalactic GCs, as well as NSCs in the local Universe. In both models we observe a tail of light IMBHs (∼100 M_⊙) in clusters heavier than ∼10⁵ M_⊙. These IMBHs are retained as part of IMRIs with merging times exceeding a Hubble time.

The clear difference between the two seeding models at z = 0 raises several interesting points. First, it highlights the limited role of runaway stellar collisions in favoring the growth of IMBHs heavier than 10⁴–10⁵ M_⊙. Second, it clearly shows the low efficiency of the hierarchical BBH merger scenario in explaining low-redshift IMBHs observations, except in very massive and dense hosts (≳10⁷ M_⊙ pc⁻³). Third, it suggests that a hybrid scenario in which an IMBH seed forms from a few collisions among massive stars (e.g., triggered by binary-single and binary-binary interactions) and possibly further grows by merging with lighter BHs, as in Model B, may be a promising dynamical formation channel for IMBHs in GC-like environments.

Fig. 3 Clusters hosting an IMBH at z = 0 in Model A (left panel) and Model B (right panel) compared to different classes of IMBH host candidate observations. We note that the masses reported for low-mass AGNs refer to the host galaxy mass (up to ∼ 100 times higher than the mass of the galactic nuclei). The gray dashed line shows the theoretical prediction for the initial mass of an IMBH seeded via stellar collisions (Eq. (3) in the left panel and Eq. (5) in the right panel).

4 Discussion

In the following we present some practical applications of our models to the study of IMBHs in YCs and GCs. In particular, we analyze the implications of premature IMBH ejection from YCs on wandering IMBHs in the Milky Way and other galaxy halos, and investigate the origin of Galactic GCs through potential IMBH signatures. These discussions will form the basis for more detailed studies in forthcoming papers.

4.1 Young clusters and wandering IMBHs

As discussed in Sect. 3.3.2, metal-poor YCs with short relaxation time can form IMBHs via stellar collisions. However, many of these IMBHs are likely to be ejected shortly after formation due to Newtonian recoils, post-merger relativistic kicks, or the evaporation of their host clusters. In this section we briefly explore whether a population of IMBHs formed in Milky Way young and open clusters could remain trapped in the Galactic potential and now wander across the Milky Way.

To explore this scenario and filter out statistical fluctuations, we resimulated 5 × 10⁷ BBHs in YCs, focusing on the runaway collision scenario (Model A) and the mild collision scenario (Model B). In all models, IMBHs are either ejected during their build-up or are retained until the cluster evaporates via internal relaxation and tidal disruption operated by the galactic potential. Upon cluster dispersal, the IMBH is gently released in the galactic environment, where it will keep following the original orbit of its former host up to a Hubble time. Clusters formed sufficiently close to the center of the galaxy can also drift inward and bring their IMBHs to the galactic nucleus, a process that possibly impacts both the growth of NSCs and SMBHs (Tremaine et al. 1975; Capuzzo-Dolcetta 1993; Ebisuzaki et al. 2001; Portegies Zwart et al. 2006; Antonini et al. 2012; Arca Sedda et al. 2015; Arca Sedda & Gualandris 2018).

Depending on the ejection mechanism, the IMBH is released in the field with a different velocity. IMBHs ejected via GW relativistic kicks receive a median kick of 20 km/s, with tails up to 10³ km/s, those ejected through dynamical interactions attain velocities of at most a few km/s, while for those released by evaporating clusters we assume a velocity ∼0 km/s. The ejection mechanism also determines whether the IMBH is expelled as a single object (e.g., in the aftermath of a BBH merger) or as a binary (in case of dynamical ejections and evaporated clusters). Regardless of the ejection mechanism, IMBH ejection velocities are generally much lower than the typical escape velocity estimated for the present-time Milky Way disk (∼500 km/s; see for reference Hunt & Vasiliev 2025), thus suggesting that expelled IMBHs should remain anchored to the Milky Way gravitational field.

To constrain the number of IMBHs produced in YCs and now wandering the Milky Way, we assume that the Galactic potential at the time of ejection resembles its present-day configuration. As a reference, we adopt the MWPotential2014 model implemented in galpy (Bovy 2015), which combines the contributions of a bulge, a disk, and a dark matter halo.

Observations of Milky Way disk clusters suggest that ∼19 000 YCs with masses ≳10⁴ M_⊙ may have formed over ∼10 Gyr (Larsen 2009). Using our catalogs, we generate 5000 realizations of this population, drawing 19 000 YCs per realization with initial masses above 10⁴ M_⊙. Each YC (IMBH) is assigned a random position in the galactic plane, and the number of IMBHs retained is determined by comparing their final recoil or kick to the Galactic escape velocity. Our analysis suggests that 273 ± 8 (26 ± 2) wandering IMBHs form in Model A (Model B), among which roughly ∼20–25% are in IMBH-BH systems (48 ± 4 for Model A and 8 ± 1 for Model B). In both models, the majority of these IMBHs have masses around ∼350 M_⊙ and ejection velocities v_ej < 10² km/s.

The potential existence of these IMBHs in the Galaxy could have interesting observational implications. For instance, single or binary IMBHs might be detectable through microlensing (Paczynski 1986). IMBHs ejected from the cluster could capture luminous companions, allowing astrometric measurements to reveal their presence in the binary, as demonstrated by Gaia’s discovery of BH–star systems (El-Badry et al. 2023a,b; Chakrabarti et al. 2023; Tanikawa et al. 2024; Gaia Collaboration 2024). Furthermore, ejected IMBH–BH binaries may emit low-frequency gravitational waves detectable by interferometers such as LISA (see, e.g., Arca Sedda et al. 2021).

Fig. 4 Comparison of simulated IMBH candidates in Model B with observational data. (a) Corner plot showing the GCs (purple) and NSCs (orange) clusters in our simulations hosting an IMBH at z = 0, compared to observations of potential IMBHs in Milky Way globular clusters and G1. (b) One-dimensional IMBH mass distribution, compared to estimated masses for G1 and 47 Tuc.

4.2 Galactic globular clusters hosting IMBHs and their (mis)classification

As discussed in Sect. 3.3.2, stellar collisions can explain potential IMBH observations at z = 0 in clusters with masses below ∼10⁷ M_⊙. However, the agreement between our simulations and candidate GCs hosts (Lützgendorf et al. 2013b; Kızıltan et al. 2017; Häberle et al. 2024) depends strongly on the assumed IMBH seeding model. Figure 3 shows that a runaway stellar collision scenario (Model A) struggles to reproduce these observations, whereas a milder stellar collision scenario (Model B) aligns significantly better with the data.

Focusing on Model B, we find that both GCs and NSCs can partially recover the same host–IMBH masses inferred for these observations (Fig. 4). To quantify this, we introduce a simple Bayesian framework to evaluate whether local IMBH host candidates are more likely GCs or NSCs, offering a systematic approach for studying the Milky Way’s GCs. This is particularly relevant since the Galactic population of GCs is expected to be mixed. Theoretical models and observations indeed suggest that a fraction of the Milky Way’s GCs formed in situ, while ∼30– 35% were accreted during galaxy mergers (see, e.g., Massari et al. 2019), with some of the latter possibly being stripped galactic nuclei. Since our simulations reveal clear links between the properties of a cluster and its IMBH, they provide a powerful tool to differentiate between genuine GCs and former NSCs among local IMBH hosts.

To summarize, our goal is to use the Bayes factor to compare two competing hypotheses for the nature of the observed IMBH hosts:

1. H_GC: the observed IMBH host is a GC;

2. H_NSC: the observed IMBH host is a NSC.

As already mentioned, we restrict our analysis to Model B simulation (right panel of Fig. 3), which shows the most significant overlap with the observed candidates. From this simulation we obtain discrete samples of the logarithm of the intrinsic cluster mass and IMBH mass, i.e., {log(M_cl), log(M_IMBH)}. To convert these samples into continuous probability densities we smooth them using a gaussian kernel density estimator (KDE). Figure 4 shows a zoom-in of these simulated distributions together with the inferred cluster and IMBH masses of a set of known Milky Way and local GCs host candidates (black squares). These masses are also listed in Table 5. We gather these estimates from a collection of different studies (Lützgendorf et al. 2013b; Kızıltan et al. 2017; Häberle et al. 2024; Bañares-Hernández et al. 2025) as follows. For the clusters with only upper limits on IMBH masses in Lützgendorf et al. (2013b), we conservatively treat these limits as representative values. For ω Centauri (ω Cen), we rely instead on two recent and somewhat conflicting analyses. Häberle et al. (2024) infer a lower limit of ∼8200 M_⊙ based on stellar kinematics, whereas Bañares-Hernández et al. (2025) favor an extended central mass of ∼2–3 × 10⁵ M_⊙, setting a lower limit of ∼6000 M_⊙ for an IMBH. To bracket these uncertainties, we model the IMBH mass in ω Cen as uniformly distributed between 6000 and 8200 M_⊙. Finally, for 47 Tucanae (47 Tuc) we adopt the cluster and IMBH mass estimates reported by Kızıltan et al. (2017). We emphasize, however, that the existence of an IMBH in 47 Tuc is still controversial, as multiple studies have demonstrated that a highly concentrated population of stellar-mass BHs in the cluster core can account for its observed dynamical properties without invoking an IMBH (see, e.g., Mann et al. 2019; Hénault-Brunet et al. 2020). For each set of inferred masses, we assign a posterior probability distribution function, which we approximate as a two-dimensional Gaussian with means and standard deviations given by the reference values in Table 5.

In Appendix E, we derive an approximated expression of the Bayes factor relative to a single observation assuming that the likelihood is effectively the same for H_GC and H_NSC, and exploiting the fact that the likelihood can be derived directly from the posterior of {log(M_cl), log(M_IMBH)} if the prior is sufficiently broad and uniform. The values of the Bayes factors’ logarithms for each host candidate are listed in the last column of Table 5 as log ${ℬ}_{\text{G}}^{\text{N}}$. We note that the rows of Table 5 are ordered according to log ${ℬ}_{\text{G}}^{\text{N}}$ in descending order, with positive (negative) values supporting a NSC (GC) hypothesis. In particular, we can consider log $B_{\text{G}}^{\text{N}}>2(<-2)$ as strongly pointing towards a NSC (GC) nature of the cluster (Jeffreys 1939). Two notable cases which emerge from our simple analysis are the G1 cluster, which orbits the Andromeda galaxy, and ω Cen in the Milky Way, both showing substantial preference for a NSC origin. G1 is hypothesized to be the disrupted nucleus of a dwarf galaxy that merged with Andromeda, a scenario supported by its broad stellar metallicity distribution (see, e.g., Ma et al. 2007). Similarly, ω Cen, the most massive cluster in the Milky Way, is also one possible example of an accreted nuclei (Hilker & Richtler 2000; Ibata et al. 2019), though the presence of an IMBH in this cluster has long been debated. Most clusters, like 47 Tuc, show a slight or negligible preference to the NSC or GC hypothesis. For reference, in Fig. 4 we also present a zoom-in of the IMBH mass distributions in GCs and NSCs and the putative masses of the IMBHs contained in G1, ω Cen and 47 Tuc. We note that, while G1 falls on the peak of the NSCs distribution found in our simulations, the ω Cen IMBH mass is more compatible with the IMBHs found in GCs. Nevertheless the large mass of this cluster ultimately determines its preference towards a NSC origin.

5 Conclusions

We investigated the formation and growth mechanisms of IMBHs in young (YCs), globular (GCs), and nuclear (NSCs) clusters by means of the new semianalytic code B-POP (Arca Sedda & Benacquista 2019; Arca Sedda et al. 2023b; Arca Sedda et al., in prep.). We focused on the potential formation of IMBHs through two mutually nonexclusive channels, namely (i) the collapse of very massive stars (VMSs) resulting from repeated stellar collisions and (ii) hierarchical BBH mergers. For the first channel, we adopted two models: Model A, assuming VMS formation via runaway stellar collisions in clusters with short core-collapse times (t_cc < 5 Myr; Portegies Zwart & McMillan 2002), and Model B, invoking a mild collision scenario in which multiple stellar collisions in extremely dense clusters drive the VMSs formation (ρ_cl,0 > 10⁵ M_⊙ pc⁻³) and motivated by recent N-body simulations (Arca Sedda et al. 2023a; Rantala & Naab 2025). In both cases we restricted the seeding to metal-poor environments (Z < 0.001) to avoid strong stellar-wind mass loss. Furthermore, we tested a separate scenario using the same seeding as Model A, but with alternative redshift distributions for NSCs and GCs (Model A′), to examine the effects of clusters formation histories on IMBH assembly.

Our main results can be summarized as follows:

• We show that seeding conditions and cluster formation histories strongly shape the IMBH formation efficiency across cluster families. The runaway stellar collision scenario (Model A), which favors systems with short core-collapse times, maximizes the efficiency in GCs. Conversely, the mild collision scenario (Model B) is the most efficient in NSCs, which comprise the densest clusters. In Model A′ NSCs prove to be much more efficient at producing IMBH seeds compared to GCs, owing to their broader distributions in redshift and metallicity, which extends down to Z < 0.001. Hierarchical IMBHs form efficiently only in NSCs;

• We observe a marked gap in the cluster mass–IMBH mass plane at M_cl ∼ 10⁷, M_⊙ for the population of IMBHs retained at z = 0 in Model A. This scarcity of IMBH hosts results from the combined effects of the seeding conditions and the inefficiency of the hierarchical channel IMBH formation channel. In contrast, Model B presents a homogeneous population of retained IMBHs which smoothly transitions from IMBHs seeded by stellar collisions to those assembled hierarchically via BBH mergers. These IMBHs have a clear overlap with potential observations of IMBHs in galactic GCs.

Additionally, we discussed possible implications for IMBHs in YCs and GCs:

• We showed that if stellar collisions are an efficient mechanism for IMBH formation, we likely expect a subpopulation of ejected IMBH seeds wandering in the potential well of Milky Way-like galaxies. The presence of such IMBHs (as single objects or in binaries) could be probed through different observational channels, such as astrometric measurements of a potential stellar companion (similarly to Gaia-BH1, Gaia-BH2, and Gaia-BH3 El-Badry et al. 2023a,b; Chakrabarti et al. 2023; Tanikawa et al. 2024; Gaia Collaboration 2024), microlensing (Paczynski 1986), or GW emission signatures (Arca Sedda et al. 2021);

• We discuss a simple Bayesian framework that exploits the correlations found in our simulations between the clusters and IMBHs masses to probe the nature of Milky Way IMBH host candidates. This is particularly interesting as some of the Milky Way GCs (e.g., ω Cen) are believed to actually be stripped nuclei of galaxies that interacted with our galaxy in the past. We apply this approach comparing the systems found in Model B to the inferred masses in Lützgendorf et al. (2013b); Kızıltan et al. (2017); Häberle et al. (2024); Bañares-Hernández et al. (2025). We find that most clusters lie in a region where the GCs and NSCs hypotheses are equally probable with the exception of the G1 cluster in the Andromeda galaxy and ω Cen, which strongly support a NSC origin.

Overall, our work highlights the critical role of stellar collisions in producing IMBHs over a broad range of clusters. Moreover, we propose possible ways of probing the presence of such IMBHs and of disentangling their formation mechanisms and the formation histories of their hosts. We also emphasize the fundamental importance of simultaneously considering the codependences among different formation mechanisms, cluster families, and their evolutionary histories, as well as environmental factors, to effectively investigate IMBH production and growth. In this context, the B-POP code represents a valuable tool for these kinds of studies as it integrates all these aspects within a flexible and efficient framework.

Acknowledgements

The authors acknowledge Gor Oganesyan, Filippo Santoliquido, and Konstantinos Kritos for useful discussion and comments. MAS acknowledges funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 101025436 (project GRACE-BH). MAS and CU acknowledge financial support from the MERAC Foundation under the 2023 MERAC prize and support to research. This research made use of NUMPY (Harris et al. 2020), SCIPY (Virtanen et al. 2020), PANDAS (Pandas Development Team 2020). For the plots we used MATPLOTLIB (Hunter 2007).

References

Appendix A B-POP overview

We briefly outline the typical B-POP workflow for a purely dynamical population of BBHs, i.e., BBHs formed and evolving in star clusters, and review the main theoretical prescriptions implemented in the code.

Appendix A.1 Workflow

At the beginning of each simulation, B-POP

1. determines the number of clusters from the total number of binaries and the relative abundances of YCs, GCs, and NSCs clusters;

2. assigns cluster masses and half-mass radii by sampling their initial distributions;

3. assigns each cluster a formation redshift and metallicity;

4. iterates over the cluster population;

5. checks whether each cluster can form an IMBH via stellar collisions (Sect. 2.2); if so, an IMBH is assigned, otherwise a stellar-mass BH is sampled from the MOBSE catalogues;

6. assigns a spin to the central BH or IMBH;

7. computes the total number of BHs available for dynamical interactions and mergers (see more in Appendix A.2);

8. follows the dynamical evolution of the BH (or IMBH) while self-consistently evolving the cluster mass and half-mass radius;

9. moves to the next cluster–BH pair if the BH (or IMBH) is ejected, fails to merge within a Hubble time, or if the host cluster evaporates.

At the end of the simulation, the code provides a realistic sample of the dynamical population of BBH mergers forming in the Universe from YCs, GCs, and NSCs.

Appendix A.2 BH dynamics

In the following we present a more focused overview of the modeling of BH dynamics in star clusters implemented in B-POP. We refer to our companion paper (Arca Sedda et al., subm. to Phys. Rev. D.) for a detailed description of the clusters mass and half-mass radius evolution modeling.

Given a cluster of initial mass M_cl,0, B-POP firstly sets up the number of possibly interacting BHs per cluster (n_BHs) assuming it to be comparable to the number of BHs inside the cluster scale radius, i.e., $$n_{\text{BHs}}=f_{\text{BH}}{f}_{\text{reten}}{f}_{\text{encl}}{M}_{\text{cl},0},$$(A.1)

with f_BH = 0.0008 the fraction of cluster mass in stellar BHs according to a Kroupa (2001) IMF, f_reten = 0.5 the fraction of BHs retained in the cluster, and f_encl the fraction of cluster mass enclosed within the cluster typical radius.

For each cluster the code follows the evolution of a single BH (potentially an IMBH seed) and assumes the reservoir of surrounding BHs to be n_BHs. The initial absolute time considered by the code is the formation time of the cluster which varies depending on the cluster type and the assumed formation history (see also Sect. 2.3). To this timescale B-POP adds up the BH stellar progenitor lifetime, taken from the input stellar catalogues. We note that the BH can be already in a binary with a certain probability which is specified as one of B-POP input parameters. Afterwards, B-POP considers the BH or BBH segregation timescale to the center of the cluster which was already introduced in Eq. (2). If the BH is not part of a binary, once it segregates to the core it can pair up with another BH via either three-body or binary-single scatterings. In the first scenario the BBH forms as a result of the dynamical interaction of three unbound objects over a timescale (Lee 1993) $$\begin{array}{cc}& t_{\text{3bb}}=4\text{Gyr}\left(\frac{{10}^6{\text{M}}_{\odot }{\text{pc}}^{-3}}{{\rho }_{\text{c}}}\right){\left({\zeta }^{-1}\frac{{\sigma }_{\text{c}}}{30\text{km/s}}\right)}^9\times \\ & {\left(\frac{10m_{\ast }}{m}\right)}^{9/2}{\left(\frac{10{\text{M}}_{\odot }}{m}\right)}^{-5},\end{array}$$(A.2)

where ρ_c and σ_c are the core density and velocity dispersion, while ζ ≤ 1 parametrizes the level of energy equipartition between the heavy and light populations of stars in the cluster. In the second scenario the BH is captured by a pre-existing stellar binary on a timescale (Miller & Lauburg 2009) $$\begin{array}{cc}{t}_{\text{bs}}=3\text{Gyr}\left(\frac{0.01}{f_{\text{b}}}\right)& \left(\frac{{10}^6{\text{pc}}^{-3}}{n_{\ast }}\right)\times \\ & \left(\frac{{\sigma }_{\text{c}}}{30\text{km/s}}\right)\left(\frac{10{\text{M}}_{\odot }}{\left(m+m_1+m_2\right)}\right)\left(\frac{1\text{au}}{a_{\text{h}}}\right).\end{array}$$(A.3)

In Eq. (A.3), m₁ and m₂ correspond to the masses of the two stars, f_b to the fraction of primordial stellar binaries in the cluster, n_∗ to the star number density, while the orbital separation $$a_{\text{h}}\sim 59\text{a}.\text{u}.\left(\frac{m_1+m_2}{30{\text{M}}_{\odot }}\right)\left(\frac{30\text{km/s}}{{\sigma }_{\text{c}}}\right)$$

represents the typical hard binary separation.

A hard binary is characterized by a binding energy greater than the typical kinetic energy of stars in the cluster (G m₁m₂/a > m_∗ σ²). On the contrary, a binary with a binding energy smaller than the average stellar kinetic energy is referred as to a soft binary.

A hard binary will continue to harden (to get tighter) on average as a result of close dynamical interactions with the other objects in the cluster, whereas a soft binary will continue to soften (Heggie 1975). BHs are typically in hard binaries since they are among the most massive objects in a cluster.

Though dynamical interactions favor the shrinking of the orbital separation they can also result in the premature expulsion of the BBH from the cluster. These interactions also give off some energy to the binary center of mass, producing a recoil which gets increasingly stronger as the binary gets harder. In B-POP we evaluate the critical orbital distance at which the binary recoil exceeds the escape velocity of the cluster and ejects the binary as (Antonini & Rasio 2016) $$a_{\text{ej}}=0.07\text{a}.\text{u}.\frac{\mu m_{\text{p}}^2}{\left(m_1+m_2+m_{\text{p}}\right)\left(m_1+m_2\right)}{\left(\frac{v_{\text{esc}}}{50\text{km}/\text{s}}\right)}^{-2},$$(A.4)

where m₁, m₂ are the masses of the BHs in the binary, µ = m₁m₂/(m₁ + m₂) is the binary reduced mass, and m_p is the average mass of a BH perturber. This dynamical hardening competes with the hardening driven by GW emission which becomes efficient at an orbital separation: $$a_{\text{GW}}=0.05\text{a}.\text{u}.{\left(\frac{m_1+m_2}{20{\text{M}}_{\odot }}\right)}^{3/5}\left(\frac{{\left(m_2/m_1\right)}^{1/5}}{{\left(1+m_2/m_1\right)}^{2/5}}\right)\times {\left(\frac{{\sigma }_{\text{c}}}{30\text{km/s}}\right)}^{1/5}{\left(\frac{{10}^6{\text{M}}_{\odot }{\text{pc}}^{-3}}{{\rho }_{\text{c}}}\right)}^{1/5}.$$(A.5)

The code computes the orbital separations in Eq. (A.5) and (A.4) and considers the maximum of the two as the critical separation orbit, i.e., a_crit = max(a_ej, a_GW). This orbital separation is assumed to be reached after a timescale (Miller & Hamilton 2002) $$t_{\text{crit}}\sim 5\frac{m_1+m_2}{m_{\text{p}}}{t}_{2-1},$$(A.6)

with t₂₋₁ the typical timescale over which the binary hardens due to single-binary interactions, which is evaluated in B-POP as (Gültekin et al. 2004; Antonini & Rasio 2016) $$\begin{array}{cc}& t_{2-1}=0.02G\text{yr}{\zeta }^{-1}\left(\frac{{10}^6{\text{pc}}^{-3}}{n_{\ast }}\right)\left(\frac{{\sigma }_{\text{c}}}{30\text{km/s}}\right)\times \\ & \left(\frac{0.05\text{a}.\text{u}.}{a_{\text{h}}}\right)\left(\frac{20{\text{M}}_{\odot }}{m_1+m_2}\right)\sqrt{\frac{10m_{\ast }}{m_{\text{p}}}}.\end{array}$$(A.7)

If a_crit = a_ej B-POP resamples the BBH eccentricity from a thermal distribution and checks if t_crit is smaller than the GW merger time (t_GW) evaluated as in Peters & Mathews (1963). If t_GW < t_crit the code assumes the binary to merge within the cluster on a merger time t_mer = t_GW. This corresponds to a scenario in which the BBH merger is caused by a sudden eccentricity change induced by strong dynamical interactions. If t_GW > t_crit the binary is considered to be ejected and to merge over a timescale t_crit + t_GW. Conversely, whenever a_crit = a_GW B-POP assumes the BBH to merge within the cluster on a timescale t_crit + t_GW. If the total time needed to reach the merger is longer than a Hubble time the binary is considered to not be merging. When the BBH successfully merges within the host cluster the code evaluates the relativistic kick induced on the remnant BH and compares it to the cluster escape velocity at the moment of the merger. If the kick does not exceed the cluster escape velocity, the BH can undergo subsequent cycles of evolution and experience multiple mergers, until it is (i) ejected by a Newtonian recoil in a binary interaction or by a post-merger relativistic kick as a single object, (ii) retained in a binary that does not merge within a Hubble time, or (iii) left without companions after exhausting all other BHs in the cluster.

Appendix B Depletion of massive stars due to stellar collisions and minimum cluster mass

Stellar collisions typically involve the most massive stars in clusters, both because of their large collision cross sections and because of their overabundance in the cluster core due to mass segregation. Consequently, stellar collisions can deplete the cluster of potential stellar-mass BH progenitors, potentially impacting the subsequent growth of IMBH seeds via BBH mergers. In this section we show that this depletion process does not significantly affect the population of stellar-mass BH progenitors at the low metallicities required to sustain the growth of VMSs through stellar collisions.

We consider a Kroupa-like initial mass function (IMF; Kroupa 2001) defined as $$\xi (m):=\frac{dN}{dm}=\{\begin{array}{cc}A\times m^{-\alpha }& m\in [0.08,0.5]{\text{M}}_{\odot }\\ B\times m^{-\beta }& m\in [0.5,150]{\text{M}}_{\odot }\end{array},$$(B.1)

with α = 1.35, β = 2.35, and normalization constants A, B chosen such that ∫ξ(m) × m dm = 1 and ξ(m) is continuous at m = 0.5 M_⊙. Assuming that stellar collisions preferentially consume stars from the heaviest to the lightest, we estimate the minimum mass cutoff m_cut such that the fraction of stellar mass above m_cut equals the fraction of mass accreted through stellar collisions. While this is a simplification, since collisions are stochastic and may involve lighter stars, it allows us to identify the stellar mass range that primarily contributes to the final VMS mass. The fraction of mass accreted by a seed star via stellar collisions can be simply expressed as $$f_{\ast }=\frac{m_{\text{VMS}}-m_{\text{seed}}}{M_{\text{cl}}},$$(B.2)

with m_VMS (m_seed) the final (seed) mass of the VMS. In a runaway stellar collisions scenario (Model A) this fraction can be computed using Eq. (3) and corresponds to a value f_∗,A ∼ 0.008. Equating f_∗,A to the integral of the high-mass tail of the IMF, $${\displaystyle {\int }_{m_{\text{cut}}}^{150{\text{M}}_{\odot }}\xi (m)\times mdm}$$(B.3)

we find m_cut ∼ 110 M_⊙. For metallicities Z < 10⁻³, this mass lies near the low-mass edge of the pair-instability supernovae (PISNe) gap (see Spera & Mapelli 2017; Iorio et al. 2023 for reference), although the exact gap boundaries depend on stellar evolution models. Therefore, we reasonably neglect depletion in our simulations, assuming that most of the VMS final mass originates from stars that would otherwise have undergone PISNe.

Similarly, for the mild collisions seeding scenario (Model B), we can estimate the maximum expected mass fraction from m_VMS,max in Eq. (5), f_∗,B ∼ 0.02, and find a cutoff mass from Eq. (B.3) m_cut ∼ 93 M_⊙, which indicates a mild depletion of massive stars. We note, however, that IMBH seeds typically do not merge with all available BHs in the cluster, hence we expect that neglecting the depletion in this case would still not affect significantly our final results.

Our calculations implicitly assume that the cluster stellar population closely follows the theoretical IMF. However, if the cluster contains too few stars, the high-mass end of the IMF will be poorly sampled. To estimate the minimum cluster mass required to host at least one star heavier than m_cut, we compute the average number of stars in the cluster as $N_{\ast }=\frac{M_{\text{cl}}}{\langle m_{\text{Kroupa}}\rangle }\approx \frac{M_{\text{cl}}}{0.7{\text{M}}_{\odot }}$, and solve for M_cl such that $N_{\ast }\times {\displaystyle {\int }_{m_{\text{cut}}}^{150{\text{M}}_{\odot }}\xi (m)dm=1}$. This yields a minimum cluster mass of approximately 2 × 10⁴ M_⊙ for the runaway seeding scenario (Model A) and 5 × 10³ M_⊙ for the alternative seeding scenario (Model B). Accordingly, we allow seeding only in clusters heavier than these thresholds for the respective models.

Appendix C Mergers and IMRIs

In Sect. 3.2 we have studied the IMBH formation efficiency on the overall population of simulated BHs. In the following we focus on the subset of these BHs which undergo at least one merger in our simulations and discuss the frequency of mergers containing an IMBH in the different models. In Table C.1 we present the different merger efficiencies (ϵ_BBHs) we find for YCs, GCs, and NSCs. These efficiencies represent the fraction of BHs that experience at least one merger out of the total simulated BH population. Their values are the same for Model A and B but change slightly in Model A′ for GCs and NSCs. This indicates that the choice of seeding prescription has only a minor impact on the overall number of merging BHs. In contrast, differences in the assumed formation histories for GCs and NSCs play a more significant role. In Table C.1 we also report the fraction of merging BHs which ultimately result in the merger between an IMBH and a stellar-mass BH, also known as intermediate-mass ratio inspirals (IMRIs).¹ These systems are particularly interesting since they can represent potential GW sources detectable by the next generation detectors in different frequency bands from mHz, for example LISA (Colpi et al. 2024), to deci-Hz (Arca Sedda et al. 2020a), and also the Lunar Gravitational-Wave Antenna (Ajith et al. 2025), to Hz and beyond, for example the Einstein Telescope (Abac et al. 2025a). We evaluate the fraction of IMRIs on the BBH merger population for the different cluster families, along with the median number of mergers (⟨N_mer⟩) the progenitor IMBH seed or BH undergoes in that specific environment, and the subfractions of IMRIs in which the primary is heavier than 150 M_⊙ or 10³ M_⊙, respectively.

As expected, IMRIs typically arise from IMBH seeds, with relative abundances up to two orders of magnitudes larger than the ones observed for the hierarchical channel depending on the environment and on the model considered. More specifically, the hierarchical channel is extremely inefficient compared to the seeding one in YCs while the two channels produce comparable amounts of IMRIs in NSCs. Significant differences in the IMRI fractions also emerge between the two seeding models (Model A and Model B), suggesting that future detections of these systems could serve as a valuable diagnostic for discriminating among different seeding mechanisms. The actual detectability of these sources, however, will naturally depend on the sensitivity of specific interferometers; evaluating event rates for different instruments will be our next step to translate these findings into observational forecasts.

Seeded IMBHs often undergo many more mergers (on average) than hierarchical IMBHs. Indeed, IMBH seeds are much heavier at birth than the underlying stellar-mass BH populations in the host clusters and are less prone to ejection (see also Sect. 3.3.1), hence sustaining longer merger chains. This difference across the seeded and hierarchical channel is particularly clear in Model B, in which IMBHs seeds form in very dense environments which favor the retention and fast mergers of the seed with other BHs. In NSCs half of the IMBH seeds can exceed ∼50 (∼240) mergers in Model A (B). In contrast, in YCs IMBHs seeds typically undergo only one merger before being expelled from their host. Similarly, hierarchical IMBHs forming in YCs and GCs typically emerge from a single BBH merger of relatively heavy stellar-mass BHs. Most of these IMBHs present masses at the lowest end of the IMBH mass spectrum (∼100 M_⊙), hence no IMRI is found with primary mass above 150 M_⊙ in these environments ( f_>150M_⊙ = 0). High-mass IMRIs (> 10³ M_⊙) occur at a rate of approximately one per 10²–10⁴ BBH mergers in our models. The production of such heavy binaries is especially interesting for YCs and GCs, where, as we have previously mentioned, the pure hierarchical formation channel is otherwise highly inefficient.

Appendix D Gap dependence on the initial conditions

As discussed in Sect. 3.3.2, we observe a dearth of retained IMBHs at z = 0 for clusters with masses around ∼10⁷ M_⊙ when assuming the seeding to be efficient only for clusters with very short core-collapse times, namely < 5 Myr (Model A). This gap in the cluster mass–IMBH mass plane, which we do not observe in Model B, is caused by the combined effect of the seeding condition and the initial mass–half-mass radius distributions assumed for the clusters. To assess the effect of assuming different initial conditions on the gap we run a separate simulation of 10⁸ BHs in GCs where we assume (i) the masses of the clusters to be enhanced by 0.5 dex and (ii) the half-mass radii to be decreased of 1 dex. These conditions are chosen to artificially enhance the number of clusters in the runaway seeding regime (see Eqs. (2) and (1)). In Fig. C.1 we compare the obtained population of retained IMBHs to the one obtained for IMBHs in GCs in Model A (see Fig. 3). As expected, we find an increase of IMBH seeds hosts. We also observe an increase of systems above the M_IMBH/M_cl = 1 contour, since the overall decrease in the GCs relaxation times corresponds to higher mass losses over the simulation time. Interestingly, we also observe an enhanced formation of IMBHs through hierarchical mergers of stellar-mass BHs, corresponding to the points lying below the tail of seeded IMBHs and mainly clustered at M_cl ≳ 10⁵ M_⊙. These IMBHs are more easily produced under these initial conditions thanks to the increase in the escape velocity of GCs, which is of the order of $\sim \sqrt{M_{\text{cl,}0}/R_{\text{h},0}}=0.75\text{dex}\sim \mathrm{5.6.}$. These IMBHs mildly overlap with galactic GCs host candidates though their abundance is still very limited with respect to the systems found in Model B (right panel of Fig. 3).

Fig. C.1 GCs hosting an IMBH at z = 0, compared to observed host candidates in local GCs (squares; Lützgendorf et al. 2013b) and a potential TDE (right triangle; Lin et al. 2018), for two sets of initial mass–half-mass radius conditions. The orange-edged star marks the ω Cen cluster (Häberle et al. 2024). As in Fig. 3, the dashed lines correspond to mass ratios MIMBH/Mcl = 10−1, 10−2, 10−3. The top panel shows the results from a B-POP simulation of 108 GCs, while the bottom panel shows the results from a simulation in which the same GCs have initial masses higher by 0.5 dex and half-mass radii smaller by 1 dex.

Appendix E Bayes factor for cluster classification

In Sect. 4.2 we employed the Bayes factor to assess the relative support for two competing hypotheses regarding the nature of some GCs IMBH host candidates (Lützgendorf et al. 2013b; Kızıltan et al. 2017; Häberle et al. 2024; Bañares-Hernández et al. 2025): a globular cluster origin (H_GC) or a nuclear star cluster origin (H_NSC). The Bayes factor is defined as $${ℬ}_{\text{G}}^{\text{N}}=\frac{\mathcal{Z}\left(d\mid H_{\text{NSC}}\right)}{\mathcal{Z}\left(d\mid H_{\text{GC}}\right)},$$(E.1)

where $$\mathcal{Z}\left(d\mid H_i\right)={{\displaystyle \int }}^{\text{}}d\theta ℒ\left(d\mid \theta ,H_i\right)\pi \left(\theta \mid H_i\right),H_i\in \left\{H_{\text{GC}},H_{\text{NSC}}\right\}$$(E.2)

denotes the marginal likelihood, or evidence, of the data d under hypothesis H_i. The evidence is the likelihood of the data d given the intrinsic parameters θ, averaged over their prior distribution π(θ|H_i) under hypothesis H_i. In our case, π(θ|H_i) corresponds to the two–dimensional distributions of cluster and IMBH masses obtained from Model B simulation (see Fig. 4).

Since our simulations do not provide direct access to ℒ(d, |, θ, H_i), we instead apply Bayes’ theorem to rewrite the likelihood in terms of the observed posterior probability distributions. These posteriors are approximated as Gaussians based on the values reported in Table 5 and we can denote them as p(θ | d, H₀), where H₀ represents the specific model assumptions adopted in each observational study. By Bayes’ theorem, $$ℒ\left(d\mid \theta ,H_0\right)=\frac{p\left(\theta \mid d,H_0\right)\mathcal{Z}\left(d\mid H_0\right)}{\pi \left(\theta \mid H_0\right)},$$(E.3)

with π(θ | H₀) the prior adopted in the corresponding inference, and 𝒵(d | H₀) the associated evidence.

In general, the likelihood may depend explicitly on the model assumed, i.e., ℒ (d | θ, H₀) ≠ ℒ (d | θ, H_GC) ≠ ℒ (d | θ, H_NSC). In the present context, however, the clusters observables (e.g., velocity dispersions, luminosity profile) are governed by the same dynamical and stellar–evolutionary processes, regardless of their nature, whether globular or nuclear. Consequently, we can assume a common likelihood ℒ (d | θ) = ℒ (d | θ, H₀), and attribute all model dependence to the prior π(θ|H_i). The prior π(θ|H₀) entering Eq. (E.3) varies across observational studies and it is often not explicitly specified. For tractability, we approximate it as broad and uniform in both log₁₀(M_cl) and log₁₀(M_IMBH). Hence, substituting Eq. (E.3) into Eq. (E.1), both 𝒵 (d | H₀) and π(θ|H₀) cancel out and the Bayes factor reduces to $${ℬ}_{\text{G}}^{\text{N}}\sim \frac{{{\displaystyle \int }}^{\text{}}d\theta p\left(\theta \mid d,H_0\right)\pi \left(\theta \mid H_{\text{NSC}}\right)}{{{\displaystyle \int }}^{\text{}}d\theta p\left(\theta \mid d,H_0\right)\pi \left(\theta \mid H_{\text{GC}}\right)}.$$(E.4)

Moreover, since the priors are assumed to be broad, no relevant regions of parameter space are truncated, and the integrals can be evaluated over the full domain of π(θ|H_i). Specifically, we estimate the integrals in Eq. (E.4) using Monte Carlo summation as $${\displaystyle \sum }_j^{N_{\text{s}}}\pi \left({\theta }_j\mid H_i\right),\text{where}{\theta }_j\sim p\left(\theta \mid d,H_0\right),$$(E.5)

i.e., sampling θ_j from the host candidate properties posterior (p(θ|d, H₀)) and computing the probability of such a system emerging from our simulation as a NSC or GC (π(θ|H_NSC) or π(θ|H_GC)).

To actually select one hypothesis over the other we would need to include in Eq. (E.4) an extra term, the so-called prior-odds of the two hypotheses, i.e., π(H_NSC)/π(H_GC). Including this term we would obtain $$O_G^N=\frac{p\left(H_{\text{NSC}}\mid d\right)}{p\left(H_{\text{GC}}\mid d\right)}={ℬ}_G^N\frac{\pi \left(H_{\text{NSC}}\right)}{\pi \left(H_{\text{GC}}\right)},$$(E.6)

where $O_G^N$ is called posterior-odds. In our context, the prior odds encode the initial degree of belief about the nature of the observed hosts and can vary substantially depending on the adopted criteria. For example, one could attempt to evaluate it starting from the number of mergers the Milky Way has experienced with dwarf galaxies (i.e., the number of potential accreted nuclei) and the number of in situ GCs IMBH hosts. However, these two quantities are highly nontrivial to estimate. Moreover, the prior odds may differ for each candidate depending on properties such as metallicity spread, density profile, or galactic position.

For simplicity, in this work we focus on estimating only the Bayes factor (Eq. E.4), which measures the relative support for the globular vs. nuclear hypothesis given the data. This corresponds to evaluating the posterior odds under the assumption of equal prior probabilities for H_NSC and H_GC.

All Tables

All Figures

	Fig. 1 Initial distribution of cluster formation redshifts. The different colors refer to different assumptions on the initial formation histories of GCs (solid lines) and NSCs (dashed lines). The initial redshift distribution of YCs is fixed in all models (gray solid line). The distributions are scaled such that the area subtended by the curve is 1.
In the text

Fig. 2 Left panel: initial cluster mass and half-mass radius for YCs (bottom left area), GCs (central area), and NSCs (upper right area). The cluster distributions are cut at the 99% contour. Clusters with tcc < 5 Myr lie below the dot-dashed black line; those with central density ρcl,0 > 105 M⊙ pc−3 are above the shaded dotted gray line. We also identify regions of the parameter space in which (i) a BBH binary is ejected before merging due to strong Newtonian recoils (in pink), (ii) the merger remnant is typically kicked out of the cluster due to its relativistic kick (in purple), (iii) the BBH merger time exceeds a Hubble time (in yellow). Right panel: schematic overview of a BH dynamical growth in a star cluster. The colored boxes refer to the regions highlighted in the left plot.

In the text

Fig. 3 Clusters hosting an IMBH at z = 0 in Model A (left panel) and Model B (right panel) compared to different classes of IMBH host candidate observations. We note that the masses reported for low-mass AGNs refer to the host galaxy mass (up to ∼ 100 times higher than the mass of the galactic nuclei). The gray dashed line shows the theoretical prediction for the initial mass of an IMBH seeded via stellar collisions (Eq. (3) in the left panel and Eq. (5) in the right panel).

In the text

	Fig. 4 Comparison of simulated IMBH candidates in Model B with observational data. (a) Corner plot showing the GCs (purple) and NSCs (orange) clusters in our simulations hosting an IMBH at z = 0, compared to observations of potential IMBHs in Milky Way globular clusters and G1. (b) One-dimensional IMBH mass distribution, compared to estimated masses for G1 and 47 Tuc.
In the text

Fig. C.1 GCs hosting an IMBH at z = 0, compared to observed host candidates in local GCs (squares; Lützgendorf et al. 2013b) and a potential TDE (right triangle; Lin et al. 2018), for two sets of initial mass–half-mass radius conditions. The orange-edged star marks the ω Cen cluster (Häberle et al. 2024). As in Fig. 3, the dashed lines correspond to mass ratios MIMBH/Mcl = 10−1, 10−2, 10−3. The top panel shows the results from a B-POP simulation of 108 GCs, while the bottom panel shows the results from a simulation in which the same GCs have initial masses higher by 0.5 dex and half-mass radii smaller by 1 dex.

In the text