The STARfall pathway: forging intermediate-mass black holes through stellar collisions in dense star clusters

Stellar collisions as a pathway to intermediate-mass black holes

Intermediate-mass black holes (IMBHs, $10^{2} \le m_{\mathrm{BH}}/M_{\odot} < 10^{5}$) occupy the mass range between stellar-mass black holes (BHs) and supermassive black holes (SMBHs), yet they remain the least observationally constrained population [119, 72]. A few debated candidates have been reported in Galactic and extragalactic star clusters (see e.g., [125, 24, 168, 110, 185, 186, 93, 20, 3, 164, 129, 47]). More solid evidences include the detection of an IMBH with mass $\gtrsim(6-8) \times10^{3}\,\mathrm{M}_{\odot}$ in $\omega$ Centauri [75, 22], as well as gravitational-wave (GW) observations such as GW231123, which involved a merger of two black holes with a total mass of $190-265\,\mathrm{M}_{\odot}$ [2]. In the coming years, next-generation GW observatories (e.g., Einstein Telescope [33, 1]; LISA [6, 6]; DE-CIGO [91]) will be able to observe signals from IMBH binaries and intermediate mass-ratio inspirals, i.e., an IMBH paired with a stellar mass compact object [7, 11, 10, 40]. Although their typical mass range remains largely unexplored, the next decades are expected to play a crucial role in unveiling the nature and origin of this elusive population.

From a theoretical point of view, explaining the formation of IMBHs through standard single-star evolution is challenging. In particular, pair-instability supernovae (PISNe) are expected to completely disrupt massive stars with a helium core mass between $64-135\,\mathrm{M}_{\odot}$, opening a gap in the BH mass spectrum between $\sim100-130\,\mathrm{M}_{\odot}$ [78, 28, 182, 156, 44]. Bridging this mass gap and producing more massive BHs therefore require additional formation channels (see, e.g., [52, 51, 49, 42, 41, 23]). To consider IMBHs a distinct population rather than the extreme tail of stellar-mass BHs, processes that can sustain their formation and growth across their full mass range are required. These mechanisms are historically divided into two main complementary scenarios: high-redshift baby-galactic nuclei and dense stellar clusters.

In the “baby-galactic nuclei” scenario, IMBHs can form at high redshift ($z > 10$) in the centers of proto-galaxies through collisions and mergers of Population III stars or their remnants [139, 43] (Taniwaka 2024), direct collapse of massive gas clouds [154, 35, 101], or a combination of both [159]. Population III remnants are expected to produce seeds of $m_{\mathrm{BH}} \sim10^{2}-10^{3}\,\mathrm{M}_{\odot}$, while direct-collapse channels can yield more massive seeds, $m_{\mathrm{BH}} \sim10^{4}-10^{6}\,\mathrm{M}_{\odot}$. It is possible that mergers of Population III stars in star clusters produce IMBHs up to $10^{4}\,\mathrm{M}_{\odot}$ (e.g., [173, 180, 117]). These objects may subsequently grow by accretion/mergers into SMBHs by $z \sim6-8$ [152, 177, 178, 111].

In dense stellar clusters, IMBH formation may be driven by stellar collisions and mergers [131, 130, 62, 112, 9, 11, 16, 143, 142, 17, 98, 70, 135, 138, 128]. Repeated collisions among massive stars can produce a very massive star (VMS, $m_{*} > 150\,\mathrm{M}_{\odot}$), which may collapse into an IMBH. Further growth can occur via accretion from stellar companions or mergers with stellar-mass BHs (see, e.g., [17]). Although the efficiency of IMBH seeding via VMSs depends on cluster properties such as mass segregation and core-collapse timescales (Gürkan et al. 2004; Freitag et al. 2006a,b; Mestichelli et al. 2026), runaway stellar collisions may enable the formation of IMBHs well above the PISN mass gap. This formation channel has been recently invoked as a possible route to form massive BH seeds in high-redshift compact clusters [73, 127, 181, 147, 107] (Bocchi et al. 2026).

Both $N$-body and Monte Carlo simulations indicate that stellar mergers and collisional processes in dense clusters can lead to the formation of IMBHs with masses up to $\sim500\,\mathrm{M}_{\odot}$ [112, 144, 48, 97, 68, 69, 70, 14 c]. Semi-analytic studies further suggest that repeated BH mergers may produce even more massive IMBHs, reaching $10^{4}-10^{5}\,\mathrm{M}_{\odot}$ in dense nuclear star clusters [8, 113, 98, 39, 17]. However, many of these models typically do not follow the full dynamical evolution of the host cluster and may simplify or neglect processes such as stellar collisions or additional growth channels, with only a few notable exceptions (e.g. [8, 17, 98, 167, 172, 171, 138]). Interestingly, recent studies have shown that mass loss driven by optically-thick stellar winds [155, 166], as well as collisions and mergers [134, 138, 145], can significantly influence the final mass, and consequently the fate, of VMSs.

Direct $N$-body simulations remain the most accurate approaches for modeling the evolution of dense stellar systems. Their main limitation, however, lies in the computational cost required to follow the full dynamical evolution, which restricts the number of realizations and typically yields only one collision chain per cluster. Indeed, the growth of massive stars is usually dominated by a single runaway object [131], a behavior naturally expected from mass segregation: if several stars begin to grow, they rapidly sink toward the cluster core, where repeated interactions favor the emergence of one VMS. As a result, exploring the parameter space in a statistically meaningful way remains challenging.

To overcome this limitation, we developed STARFALL, a semi-analytic tool, designed to study the formation of VMSs and IMBH seeds in dense environments. STARFALL can follow the dynamical evolution of a massive binary or of a single massive star at the cluster centre, where interactions are most frequent, while incorporating state-of-the-art prescriptions for cluster evolution and dynamical encounters. Its key strength is the computational efficiency: each binary can be evolved in seconds to minutes, depending on the initial density. This enables the exploration of a much broader parameter space, including densities currently inaccessible to direct simulations.

In this work, we investigate the conditions under which runaway collisions can lead to formation of VMSs and IMBH seeds, and the implications for high-redshift dense environments, such as Little Red Dots (LRDs; [100, 115] (Kocevski et al. 2023)) and clusters in the Cosmic Gems Arc (CGA; [4, 32, 116]), as well as for present-day globular clusters (GCs; [110, 77, 93]), using STARfall.

In the following section, we describe the implementation of STARfall. Our results are presented in Sect. 3. In Sect. 4, we discuss their implications for VMS growth and IMBH seed formation in the CGA, GCs and LRDs. Finally, we summarize our work in Sect. 5.

Dynamical evolution of massive stellar binaries with STARfall

We simulated the dynamical evolution of massive binaries, which undergo repeated binary–single interactions and stellar collisions while orbiting in the host cluster core, with the new semi-analytical tool STARfall. While STARfall can model the dynamical evolution of a single star, initial conditions assume a massive binary because $N$-body simulations show that primordial binaries are typically the systems that trigger the first collision chains [15, 136, 118]. Even in the absence of primordial binaries, dense clusters form a massive binary in the core on very short timescales (0.1 – 0.3 Myr; [137]), providing the same effective starting point for the runaway sequence. The flowchart of STARfall is shown in Figure 1. In the following, we describe the adopted star cluster model (Sect. 2.1), the treatment of single-star evolution (Sect. 2.2), and the dynamical processes (Sect. 2.3) included in the code.

Vector (PDF)

Star cluster evolution

The evolution of the star cluster is modeled using the prescriptions implemented in the B-POP³ code (Arca Sedda et al. 2026), which tracks the time evolution of the cluster’s bound and core mass, as well as its half-mass and core radii. The main equations are reported in Appendix A, while additional details can be found in [11, 17] (Arca Sedda et al. 2026). In the following, we briefly summarize the treatment of core collapse, a key physical process for our analysis. Indeed, during this phase, the steep rise in central density can trigger runaway stellar collisions that lead to the formation of VMSs (Portegies Zwart & McMillan 2002; Freitag et al. 2006a,b; Gaburov et al. 2008).

In a semi-analytic model, the main uncertainty is related to how deep the core collapses before interactions become energetic enough to revert the contraction. In STARfall, core collapse is halted once the core radius reaches a minimum value, $R_{\mathrm{cc}}$. If the central massive binary remains bound, the minimum core radius is set by the binary scale:

where $a$ is the binary semi-major axis, $m_1$ and $m_2$ are the masses of the two components, and $M_{\mathrm{core}}$ is the mass of the core. This choice reflects the fact that, below this scale, the binary dominates the central dynamics and efficiently heats the surrounding cluster [30]. If the binary is instead ionized, we adopted:

where $R_{\mathrm{core}}$ is the cluster core radius. We assume that the collapse is halted if the mass is such that the influence radius, i.e., the region within which the gravitational potential of a single star dominates over that of the cluster, becomes larger than the wander radius, i.e., the characteristic amplitude of the star’s Brownian motion around the minimum of the cluster potential [50] (Arca Sedda et al. 2026).

In this work, we considered star clusters with masses $\log(M_{\mathrm{cl,0}}/M_\odot) = 5, 6, 7, 8$, and half-mass radii $R_{\mathrm{hm,0}} = 0.01, 0.05, 0.1, 0.5, 1$, and $3$ pc, with initial core densities ranging between $\rho_{\mathrm{core,0}} \sim10^{5.5} - 10^{14}\,M_\odot\,\mathrm{pc}^{-3}$. These configurations allow us both to compare our results with initial conditions consistent with GCs and systems in the CGA, and to explore regions of parameter space that remain inaccessible to direct $N$-body calculations. For each combination of $M_{\mathrm{cl,0}}$ and $R_{\mathrm{hm,0}}$, we modeled 500 binaries, as described below.

Single-star evolution

Stellar physical sizes and post–main-sequence evolution are modeled using the following prescriptions. The fitting formulae presented below are calibrated against PARSEC stellar tracks [44], interpolated with SEVN⁴ (SEVNtracks_parsec_v05_AGB; [157, 88]) up to $m_\ast= 300\,M_\odot$. Beyond this value, the following formulae are extrapolated. We considered metallicities $Z = 0.0002$ and $Z = 0.002$.

We first initialized stellar radii at zero-age main sequence (ZAMS) using the following mass–radius relation:

and assumed it to be constant until the end of the main sequence.

To integrate stellar expansion within STARfall, we assumed that after a time:

the stellar radius increases by a factor $f_{\mathrm{rad}} = R_\ast(\tau_{\mathrm{MS}})/R_{\mathrm{ZAMS}}$, obtained from SEVN–PARSEC tracks (see lower panel of Fig. E.1). We note that the timescale over which the massive star rapidly expands is smaller than typical dynamical timescales of the cluster. The resulting radius is compared against the prescription reported in [171]:

The final adopted stellar radius is chosen as the minimum of the two. In almost all cases, the radii reported by [171] are larger than those from SEVN–PARSEC tracks.

We assumed that all stars have a lifetime:

By construction, the condition $\tau_{\mathrm{MS}}/\tau_{\mathrm{life}} < 1$ is always satisfied. The values of the fitting parameters are listed in Table 1.

Parameter	$Z = 0.0002$	$Z = 0.002$
$R_{\mathrm{ZAMS}}$ [R$_\odot$]	$a_\ast$ 0.757 ± 0.002	0.763 ± 0.004
	$b_\ast$ 0.553 ± 0.001	0.587 ± 0.001
$\tau_{\mathrm{MS}}$ [Myr]	$a_{\mathrm{MS}}$ 194 ± 4	205 ± 4
	$b_{\mathrm{MS}}$ -1.138±0.005	-1.151±0.005
	$c_{\mathrm{MS}}$ 1.752±0.006	1.766±0.006
$\tau_{\mathrm{life}}$ [Myr]	$a_{\mathrm{life}}$ 205 ± 4	213 ± 4
	$b_{\mathrm{life}}$ -1.138±0.006	-1.146±0.006
	$c_{\mathrm{life}}$ 1.954±0.007	1.951±0.006

When a merger or collision occurs between main-sequence stars, the remnant star is rejuvenated. We computed the effective age $A_{\mathrm{MS},3}$ of the remnant star, i.e., the age of a single star whose remaining main-sequence lifetime matches that of the merger product, according to Appendix § 9 of Vergara et al. [171] (but see also, e.g., Leonard [104]; Lombardi et al. [109]; Hurley et al. [85]; Glebbeek et al. [65]; Mapelli [112]; Schneider et al. [150]; Schneider [149]; Rantala et al. [138]):

where $q = m_2/m_1 < 1$, $A_{\mathrm{MS},i}$ and $\tau_{\mathrm{MS},i}$ are the effective ages and main-sequence lifetimes of the two colliding stars and the remnant star (with subscripts $i = 1, 2$ referring to the primary and secondary stars, respectively, and $i = 3$ to the merger remnant). Note that this rejuvenation prescription contains explicitly the mass ratio between colliding stars. This parameter regulates the rejuvenation itself and indicates that when the star is very massive compared to the accreted stars its age is not substantially reduced anymore. We also note that the coefficient 0.1, which depends on the mass of the hydrogen core, is calibrated at $Z = 0.01$. At the lower metallicities considered here, this coefficient should be likely higher, possibly leading to a more efficient rejuvenation and extending the stellar lifetime by a few years. A more detailed assessment of the impact of rejuvenation will be explored in future works.

The mass loss during a stellar merger or collision is parametrized as:

where $D = 0.3$ (Glebbeek et al. [64]). The stellar lifetime and radius of the remnant are computed as previously described, where in this case $m_* = (m_1+m_2)\times(1-f_{\mathrm{loss}})$.

The primary star mass loss via stellar winds is treated following the prescription of Vink [174], extended here also to stars with masses above $900\,M_{\odot}$. We did not consider mass loss from the secondary star or from other stars involved in dynamical interactions. Since the most massive stars are rapidly accreted through repeated mergers, we expect mass loss from stars other than the primary to play only a secondary role in the evolution.

The final fate of each star (i.e. whether it forms a BH and with what mass, or instead undergoes a PISN) is determined by the relation between the pre-supernova (pre-SN) stellar mass and the remnant mass derived from SEVN-PARSEC tracks (see upper panel of Fig. E.1). This choice helps to ensure consistency across different mass-loss histories. We adopted the rapid core-collapse supernova formalism of Fryer et al. [56], while PISNe follow the prescriptions of Mapelli et al. [114]. We also note that the onset of PISNe is primarily determined by the core mass rather than the pre-SN mass of the star. However, stellar collisions and the associated rejuvenation can alter the core growth, potentially leading to core masses that differ from those of single stars of the same ZAMS mass. Therefore, adopting pre-SN mass thresholds represents a conservative assumption.

Notes. The first column indicates the stellar quantities (Eqs. (3), (4), and (6), respectively), the second one lists the fitting parameters, and the third and fourth ones report their values for $Z = 0.0002$ and $Z = 0.002$, respectively.

We discuss the implications of our astrophysical assumptions in more detail in Appendices E and F. Our prescriptions are broadly consistent with Rantala et al. [138], enabling comparison across models.

Binary pairing and processes

Observations indicate that massive stars mainly live in binaries or higher-order multiple systems [146, 122]. Their shorter mass-segregation timescales and larger cross sections compared to single stars imply that massive binaries play a central role in the dynamical evolution of star clusters. In addition, typically only one star per cluster can grow through runaway collision [131, 55]. To model this scenario, we draw the primary mass from a flat distribution in the range $m_1 = 100 - 150$ M$_\odot$, while the companion mass is determined by the binary mass ratio $q$. Following Moe & Di Stefano [122], binaries are grouped into three mass-ratio intervals, $(0 - 0.5)$, $(0.5 - 0.6)$, and $(0.6 - 1)$, containing approximately 20%, 40%, and 40% of systems, respectively. One of these intervals is selected randomly based on the percentile, and the mass ratio is then drawn from a uniform distribution within the chosen range. The companion mass is finally set to $m_2 = \max(5$ M$_\odot$, $q\,m_1)$, where the minimum mass reflects observational constraints [146, 92, 94].

The binary eccentricity is drawn from a thermal distribution. The initial semi-major axis is extracted from a logarithmically flat distribution limited between the sum of the stars’ radii and 50 AU. These assumptions permitted us to adopt a setup widely used in direct $N$-body simulations (see e.g., [179, 90] (Arca Sedda et al. 2024)).

Binary-single interactions are modeled via statistical treatment based on seminal works by Heggie [79]; Hut & Bahcall [86]; Heggie & Hut [80]; Sigurdsson & Phinney [153]. This formalism allows us to determine the outcome of the interactions (i.e., hardening or softening, dynamical exchange, binary ionization, or recoil and subsequent ejection from the cluster) based on energy considerations during the encounter. All possible cases are discussed in detail in Appendix B.

Here, we summarize the formalism we employed to model these interactions. If one considers a binary with mass $m_{\mathrm{bin}} = m_1 + m_2$, semi-major axis $a$, and eccentricity $e$, traveling through a population of objects characterized by mean mass $m_{\mathrm{mean}}$, number density $n$, and velocity dispersion $\sigma$, the interaction rate between the binary and a single object in such an environment can be written as:

where $\Sigma_{\mathrm{bs}}$ is the binary cross-section:

the term in the square brackets, accounting for the gravitational focusing, is largely greater than 1 in most applications. The associated timescale of this process is given by $t_{\mathrm{bs}} = (n\,\sigma\,\Sigma_{\mathrm{bs}})^{-1}$ [121, 9].

On the other hand, the collision cross-section is computed by assuming that the pericenter distance equals the stellar radius:

where $m_i$ and $R_i$ denote the mass and radius of the binary component, while $m_p$ and $R_p$ are the mass and radius of the perturber. We deepen this assumption in Appendix § 9. The velocity dispersion and number density adopted here correspond to those of the cluster core. To determine which binary component undergoes the collision, we compared their gravitational focusing terms; the component with the larger gravitational focusing is assumed to collide. The corresponding collision timescale, $t_{\mathrm{coll}}$, is then computed following the same procedure as for the binary–single interaction case. When a collision occurs, we neglect any recoil kicks. While this may be a strong assumption, hydrodynamical simulations [46, 109, 23] indicate that recoil velocities are typically $\lesssim\sigma$. Furthermore, as the star grows in mass, the recoil velocity is expected to become progressively smaller. However, a binary may be ionized following a collision if the kinetic energy of the perturber exceeds the binding energy of the binary. The velocity of the perturber is sampled from a Maxwellian distribution with dispersion $\sigma$.

After an interaction, if the new binary pericenter is smaller than the sum of the two stars’ radii, the stars are merged and the merger remnant is possibly paired with a new companion. If the stellar merger occurs during a resonant interaction, the merger remnant is paired with the third body and the newly formed binary is assumed to have a semi-major axis $a_{\mathrm{new}} = (m_3/m_2)a_{\mathrm{old}}$, with $m_2 < m_1$ [153]. Otherwise, after a merger or ionization, we assumed that the remnant star will be captured in another binary over a time $t_{\mathrm{bin}} \equiv t_{\mathrm{3bb}}$ [102]:

where $v_\ast\sim\sigma$, $R_{\mathrm{cl}}$ is the radius of the cluster, and $n$ its number density. The binary formed this way is assigned a semi-major axis from a logarithmically flat distribution between $a_{\mathrm{new}} = (10^{-2} - 10)a_h$, $a_h = Gm_{\mathrm{mean}}/4\sigma^2$ [151], which are typical values observed in self-consistent $N$-body models, and an eccentricity extracted from a thermal distribution.

In STARfall, the scheme above is implemented following a multi-stepped strategy (see Figure 1). We assume that the cross section and the corresponding interaction timescale (either $t_{\rm bs}$ for a binary–single interaction if the binary remains bound or $t_{\rm coll}$ for a possible collision) quantify the probability of such an interaction occurring in a time interval $dt$. In particular, at each new time-step $dt$, the timescales $t_{\rm bs}$ and $t_{\rm coll}$ are evaluated. Then,

• if $t_{\rm bs} < \min(t_{\rm coll},dt)$, the system undergoes a binary–single interaction, with the outcome (e.g., hardening, dynamical exchange, or ionization) determined following Sigurdsson & Phinney [153], as described in Appendix B;

• if $t_{\rm coll} < \min(t_{\rm bs},dt)$, a collision event is assumed to occur;

• if the binary is ionized and $t_{\rm coll} < dt$, the single star can continue to experience collisions until a new binary is formed;

• multiple interactions may occur within a single $dt$; after each interaction, $t_{\rm coll}$ and $t_{\rm bs}$ are re-evaluated, and the time is advanced by $\min(t_{\rm coll},t_{\rm bs})$ until the cumulative elapsed time exceeds $dt$;

• if both $t_{\rm coll}$ and $t_{\rm bs}$ exceed $dt$, an interaction may occur with a probability $P_{\rm int}=dt/\min(t_{\rm coll},t_{\rm bs})$;

• if no interaction takes place during a timestep, the simulation advances to the next one.

If an interaction occurs, we extract the mass of the third object from a set of $N_{\rm core}=M_{\rm core}/m_{\rm mean}$ masses distributed as a Kroupa initial mass function (IMF; Kroupa [99]), between $0.08-150\,M_\odot$, with $m_{\rm mean}=0.6\,M_\odot$. Here, $N_{\rm core}$ sets the reservoir of stars available for interactions, while the actual number of stars present at each timestep is updated consistently according to the evolving core mass. We assumed that in case of any interaction, the binary will interact more likely with the most massive stars. Hence, we initially drew the ZAMS mass of the perturber in the range $m_3=20-150\,M_\odot$. In case of a merger or collision, we removed the involved object (which can be either one of the binary components or the perturber) from the pool of stars in the core. Once all massive stars are consumed or have already collapsed, perturbers are drawn from the remaining lower-end of the IMF. Our analysis considers only stars that are initially present in the core. The impact of mass segregation will be investigated in future works. It is nevertheless expected to have a secondary effect, as it occurs on timescales $\gtrsim1$ Myr for a star with $m_\star=50\,M_\odot$ (Portegies Zwart & McMillan [131]) in the most compact clusters, i.e., after the collisional runaway has already been triggered, or even on longer timescales in less compact environments.

All simulations are halted when the primary star collapses.

Results

VMS mass growth

The onset and extent of the collision chain vary with the cluster’s properties, the frequency of dynamical interactions, and the evolutionary state of the primary star. Figure 2 shows the growth of the primary star mass and the evolution of the cluster core radius in clusters with total mass $M_{\rm cl,0}=10^6\,M_\odot$ for selected half-mass radii. In the two densest clusters ($R_{\rm hm,0}=0.05,0.1$ pc), the star can begin growing mass even before core collapse, followed by a steeper mass increase during the core-collapse phase. After this phase, the mass gained through dynamical interactions becomes smaller than the mass lost via stellar winds, leading to a gradual decrease in stellar mass until the end of its life at $\sim3$ Myr. In both clusters, core depletion leads to the formation of VMSs with masses of $\sim10^4\,M_\odot$, but on different timescales. In the first cluster, core collapse occurs before the onset of post-main-sequence stellar expansion. By the time the star expands, the cluster core is already depleted. As a result, the increased stellar radius does not lead to additional collisions or a further increase in the maximum mass. In the second one, instead, stellar expansion drives further mass growth phase before the direct collapse. However, in this case, the resulting puffed-up structure may lead to enhanced collision-induced mass loss compared to the estimates of Glebbeek et al. [64] adopted in this work. In contrast, considering $R_{\rm hm,0}=0.5$ pc, we observe only a few mergers/collisions before the onset of core collapse. In this case, the core-collapse time is roughly equal to $\tau_{\rm MS}$, i.e., the time after which the star expands (see Sect. 2.2 and Eq. (4)). The increased core density and the large star’s radius maximize the interaction rate, driving a fast increase of the VMS mass. The delayed mergers rejuvenate the star efficiently over a longer timescale, prolonging its life up to $\sim5.2$ Myr. In less dense clusters, the lower interaction rate prevents, in most cases, a runaway phase, with only few stars exceeding $300\,M_\odot$.

Figure 3 shows the statistical properties of the maximum mass reached by VMSs. A comparison between the two panels highlights the key roles of stellar dynamics. We find a clear correlation between the maximum VMS mass and both the initial core density and the total mass of the cluster. Metallicity primarily affects the width of the mass distribution and limits the maximum mass by a few percent, especially for clusters with $M_{\rm cl,0}\geq10^6\,M_\odot$.

Vector (PDF)

At low cluster mass ($M_{\rm cl,0}=10^5\,M_\odot$), the VMS mass remains relatively small, in the range $\sim100-1000\,M_\odot$. In more massive clusters (e.g., $M_{\rm cl,0}=10^8\,M_\odot$), the maximum mass can exceed $10^4\,M_\odot$ and grow up to $\sim10^6\,M_\odot$ for clusters with $M_{\rm cl,0}=10^8\,M_\odot$. Such extremely massive stars may undergo general relativistic instabilities that trigger the collapse of their inner core, leading to the formation of a very massive BH surrounded by a stellar envelope feeding the accretion disk (Hosokawa et al. [81]). Our current implementation does not include a dedicated prescription to model this phase. Therefore, for simplicity, these objects are treated as ordinary single stars. While the total cluster mass sets the overall reservoir of stars available for collisions and dynamical interactions, the core density determines the interaction rate, and therefore how much of that reservoir can be accreted. The central density, indeed, regulates how efficiently and how rapidly the VMS can grow. As the initial core density increases, the maximum VMS mass rises accordingly, reflecting the enhanced frequency of close encounters and mergers. However, for $\rho_{\rm core,0}\gtrsim10^9\,M_\odot\,{\rm pc}^{-3}$, the growth saturates and both the $50\%$ and $90\%$ quantiles converge toward the median. In this regime, the primary star effectively accretes all the mass available in the cluster core.

To assess the reliability of our approach, we compared our predictions with results from direct $N$-body (Arca Sedda et al. [17]; ); Vergara et al. [170]; ); Rantala et al. [138]) and Monte Carlo (Askar et al. [21], Askar et al. in prep; Giersz et al. in prep) simulations. This comparison allows us to evaluate how well our semi-analytical approach captures the key physical processes driving runaway growth. Despite the simplifying assumptions adopted in our model, our predictions are in overall good agreement with the $N$-body results of Arca Sedda et al. [15] ($M_{\rm cl,0}=0.7\times10^5\,M_\odot$, $R_{\rm hm,0}=0.7,1.75$ pc, $Z=0.0005$) and Rantala et al. [138] ($M_{\rm cl,0}=1.5\times10^5\,M_\odot$, $R_{\rm hm,0}=0.37$ pc, $Z=0.0002,0.002$), as well as Monte Carlo simulations by Askar et al. in prep., Giersz et al. in prep. ($M_{\rm cl,0}=(1.3-6.7)\times10^6\,M_\odot$, $R_{\mathrm{hm},0}=0.25-3.3$ pc, $Z=0.0005,0.0006$, for clusters with masses and half-mass radii comparable to those explored in this work. On the other hand, the masses reported by Vergara et al. [170], both with $N$-body and Monte Carlo simulations, are up to an order of magnitude higher than our estimates for clusters with $M_{\mathrm{cl},0}=10^{6}\,\mathrm{M}_{\odot}$. Their simulations focused on clusters with $M_{\mathrm{cl},0}=(0.58-4.38)\times10^{5}\,\mathrm{M}_{\odot}$, $R_{\mathrm{hm},0}=0.01-0.05$ pc and $Z=0.0001$. As also noted by Rantala et al. [138], significant differences may arise from the adopted prescriptions for stellar-wind mass loss and stellar radii. In this work, we employed the wind prescription from Vink [174], which leads to more substantial mass loss in VMSs than the prescription of Vink et al. [175], implemented by Vergara et al. [170] with an upper limit of $1.5\times10^{-4}\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$ for stars exceeding $6\times10^{5}\,\mathrm{L}_{\odot}$ [82, 27]. Nevertheless, even with these stronger winds our model predicts the formation of stars with masses exceeding $10^{3}\,\mathrm{M}_{\odot}$, largely consistent with IMBH progenitors (see also Appendix § 12 for a detailed assessment of the impact of stellar winds). In addition, radii from SEVN-Parsec tracks are smaller with respect to the prescription implemented in BSE (see, e.g., Hurley et al. [85, 84]). Other differences may stem from the underlying cluster models and, consequently, the mass of their core. This number may be enhanced by rapid, or even primordial, mass segregation. We leave a thorough exploration of the impact of such a mechanism to a forthcoming work.

Table 2 dissects the maximum VMS mass in our models into different ranges of cluster mass ($M_{\mathrm{cl},0}$) and radius ($R_{\mathrm{hm},0}$). The table further highlights the tight relation between cluster mass, density, and terminal VMS mass. Cluster models with sizes and masses comparable to the most massive GCs in the Milky Way can nurture the growth of a VMS as massive as $900\,\mathrm{M}_{\odot}$.

	$\mathit{R}_{\mathrm{hm},0}$ [pc]
$\log(M_{\mathrm{cl},0}/\mathrm{M}_{\odot})$	0.1 – 0.5	0.5 – 1	1 – 3
5 – 6	$900^{+11000}_{-440}$	$330^{+260}_{-190}$	$127^{+19}_{-12}$
6 – 7	$16000^{+56000}_{-3500}$	$3200^{+17000}_{-3000}$	$145^{+56}_{-23}$
7 – 8	$108000^{+240000}_{-29000}$	$29000^{+21000}_{-28000}$	$189^{+120}_{-59}$

As detailed in the previous section, in STARfall the VMS can grow either via merger of binary components, hierarchical collision, or their combination. In order to assess the mutual role of such processes in the VMS buildup, Figure 4 shows how cluster mass and initial radius affect the fraction of stellar collisions, fraction of stars formed with $m_{*}>300\,\mathrm{M}_{\odot}$, and mean VMS mass. We overlay typical observed properties of Galactic GCs [77], nuclear star clusters (NSC; [58]), young massive clusters (YC; [132], as reported by Paiella et al. [128], and the young clusters observed in the CGA at a redshift $z\sim10$ [4, 116]. To keep the figure readable, the young and massive clusters present in the galaxy Henize 2-10 (with $M_{\mathrm{cl}}\sim(0.2-2.3)\times10^{6}\,\mathrm{M}_{\odot}$ and $R_{\mathrm{hm}}\sim1.2-4$ pc; [140, 124]) were not included. However, they are consistent with the most massive GCs and lower-mass clusters in the CGA. Therefore, the discussion that follows also applies to the clusters in Henize 2-10. The plot highlights that stellar collisions are more efficient than stellar mergers at producing more massive VMS in the most massive and compact clusters. Indeed, a higher stellar collision efficiency is strictly correlated with a higher probability of VMS formation. Typically, the main bottleneck for triggering a collision runaway is the merger of the initial binary. In the densest systems ($R_{\mathrm{hm},0}<1$ pc and $M_{\mathrm{cl},0}>10^{6}\,\mathrm{M}_{\odot}$), once this first merger occurs, the remnant undergoes a prolonged sequence of collisions without forming a new binary. As the central density decreases, stellar mergers become the main trigger of VMS formation, but dynamical ejection and smaller interaction rates reduce the attainable mass, leading to lighter VMS. In addition, a VMS forms in $\sim15\%-50\%$ of the models even for $t_{\mathrm{cc}}\geq5$ Myr (a typical time threshold for the onset of runaway collisions; see, e.g., Portegies Zwart & McMillan [131]), provided that the cluster is sufficiently massive ($M_{\mathrm{cl}} \gtrsim10^{7}\,\mathrm{M}_{\odot}$). In systems with $R_{\mathrm{hm},0} \gtrsim1$ pc, dynamics do not play a significant role, and stars with $m_{*} > 300\,\mathrm{M}_{\odot}$ are formed only in the most massive clusters (i.e., $M_{\mathrm{cl},0} > 10^{7}\,\mathrm{M}_{\odot}$). Comparing our models to observed clusters suggests that massive and dense GCs were able to support the growth of VMSs of $\gtrsim300\,\mathrm{M}_{\odot}$, and probably beyond $1000\,\mathrm{M}_{\odot}$. Indeed, since these old systems likely underwent significant mass loss and expansion over their lifetime, their mass-radius distribution was originally likely skewed upward and leftward in the $M_{\mathrm{cl},0}-R_{\mathrm{hm},0}$ plane ([25, 29]), i.e. in a region of the phase-space where even heavier VMSs could form. Interestingly, some of the CGA clusters seem to fall exactly in such region, where VMS mass can reach $10^{3}-10^{4}\,\mathrm{M}_{\odot}$. Our results are consistent with recent findings from Williams et al. [181], who developed an analytic cosmological framework to explore the dynamical evolution of

Vector (PDF)

Notes. Median maximum stellar mass ($m_{1,\mathrm{max}}$) as a function of initial cluster half-mass radius (columns; $R_{\mathrm{hm},0}$) and cluster mass (rows; $\log(M_{\mathrm{cl},0}/\mathrm{M}_{\odot})$) at $Z = 0.0002$. The $\pm$ values give the $1\sigma$ spread. The masses are given in $\mathrm{M}_{\odot}$ with the respective uncertainties.

dense star clusters with and without a central massive BH. In their analysis, Williams et al. [181] find that stellar collisions in sufficiently dense environments can lead to the rapid growth of VMS in high-redshift and compact clusters, reaching masses consistent with the one we find in STARFORALL, despite the different adopted approach. We further discuss this aspect in the next section.

IMBH seeds

Depending on both the mass and the metallicity, a VMS may either undergo a PISN or directly collapse into an IMBH. Table 3 reports the fraction of models that successfully form an IMBH seed. For clusters with initial mass $M_{\mathrm{cl},0} = 10^{6}\,\mathrm{M}_{\odot}$ at metallicity $Z = 0.0002$ ($Z = 0.002$), the formation of IMBH seeds occurs in only $\sim8\%$ ($\sim37\%$) of the models for $R_{\mathrm{hm}} = 1$ pc, while this fraction increases to $\sim96\%$ ($96\%$) for clusters with $R_{\mathrm{hm}} = 0.05$ pc. These cases correspond to the IMBH mass distributions shown in Figure 5. VMSs that do not collapse into an IMBH, given our initial conditions, can either be disrupted in a PISN or undergo a pulsational pair-instability supernova (PPISN), forming a stellar mass BH. The maximum stellar mass and IMBH seed mass distribution for the other clusters we considered is presented in Appendix § 10. We discuss instead the impact of our astrophysical assumptions and the stellar mass BH distribution in Appendix § 11, while the impact of the mass-loss prescriptions we adopted is deepened in Appendix § 12.

$\log(M_{\mathrm{cl},0}/M_\odot)$	$R_{\mathrm{hm},0}$ [pc]	0.05	0.1	0.5	1
	$Z = 0.0002$
5		78%	67%	57%	6%
6		96%	92%	57%	8%
7		99%	99%	71%	30%
8		96%	99%	68%	66%
	$Z = 0.002$
5		78%	72%	55%	24%
6		96%	95%	65%	37%
7		99%	99%	85%	59%
8		99%	99%	70%	63%

The distribution of the maximum stellar mass reached at the end of the collision sequence is not significantly affected by metallicity. In both cases, the distribution spreads between $\sim100\,\mathrm{M}_{\odot}-4000\,\mathrm{M}_{\odot}$ ($\sim100\,\mathrm{M}_{\odot}-400\,\mathrm{M}_{\odot}$) for $R_{\mathrm{hm},0} = 0.5$ pc ($R_{\mathrm{hm},0} = 1$ pc), while it is peaked at $\sim10^{4}\,\mathrm{M}_{\odot}$ for $R_{\mathrm{hm},0} = 0.05,\,0.1$ pc.

In contrast, both metallicity and the formation time of the VMS play an important role in determining the mass of the final remnant. In the densest clusters, the VMS forms earlier and rapidly reaches its maximum mass, above which stellar wind mass loss dominates over stellar collisions (see also Figure 2). In less dense clusters, the VMS typically forms on longer timescales and therefore has less time to lose mass before collapsing, producing either more massive BHs or BHs with masses comparable to those produced in denser environments. In particular at $Z = 0.002$, while the maximum stellar mass reached during the collision sequence strongly depends on the dynamical properties of the cluster, the IMBH mass is primarily determined by the subsequent stellar evolution. Therefore, although the densest clusters produce the most massive stars, they do not necessarily produce the most massive IMBH seeds. Indeed, as shown in Figure 5, the maximum stellar mass reaches at most $m_* \sim2 \times10^4\,M_\odot$ in the densest cluster ($R_{\mathrm{hm},0} = 0.05$ pc),

Notes. Table divided into two sections corresponding to different metallicities: $Z = 0.0002$ (top) and $Z = 0.002$ (bottom). In each section, rows represent the initial cluster mass, $\log(M_{\mathrm{cl},0}/M_\odot)$, and columns correspond to the initial half-mass radius $R_{\mathrm{hm},0}$ in pc.

while the largest IMBH mass, $m_{\mathrm{IMBH}} \sim1.5 \times10^4\,M_\odot$, is obtained at lower densities ($R_{\mathrm{hm},0} = 0.1$ pc). This trend is particularly evident at $Z = 0.002$, due to the higher efficiency of stellar wind mass loss. In this case, the IMBH mass distribution for $R_{\mathrm{hm},0} = 0.05$ pc peaks at $2 \times10^3\,M_\odot$, comparable to that for $R_{\mathrm{hm},0} = 0.5$ pc, whereas the distribution for $R_{\mathrm{hm},0} = 0.1$ pc peaks at $4 \times10^3\,M_\odot$. At $R_{\mathrm{hm},0} = 1$ pc, instead, the conditions are not always favorable for triggering a runaway collision chain. In this case, the maximum stellar mass distribution exhibits a peak at $100$–$200\,M_\odot$, corresponding approximately to the initial stellar masses and therefore indicating negligible mass growth. These stars do not contribute to IMBH formation, as they are expected to undergo PISNe. A tail extending beyond $200\,M_\odot$ remains, increasing the probability of IMBH formation to up to $\sim60\%$ in more massive clusters (see Table 3).

While the probability of IMBH formation does not show any significant dependence on metallicity in the densest clusters, at $R_{\mathrm{hm},0} = 1$ pc IMBH formation appears to be more efficient at $Z = 0.002$ (see also Table 3). At that metallicity many IMBHs are formed around a few hundred solar masses while the same IMBH subpopulation disappears at $Z = 0.0002$ (see Figure 5). This can be explained by the effects of PISNe, which prevent stars with pre-SN masses $m_{\mathrm{pre-SN}} \lesssim240\,M_\odot$ at $Z = 0.0002$ from collapsing into an IMBH (see also Fig. E.1).

Figure 6 further illustrates the statistical properties of the IMBH mass distribution as a function of the density for the two metallicities. In both cases, the highest median IMBH mass is found at intermediate core densities, $\rho_{\mathrm{core}} = 10^8$–$10^9\,M_\odot\,\mathrm{pc}^{-3}$, beyond which it decreases as the density further increases. As discussed above, this trend is particularly pronounced at $Z = 0.002$, where the more efficient mass loss due to stellar winds significantly limits the final IMBH mass, capping it at approximately $5 \times10^3\,M_\odot$. This behavior is not observed in clusters with $M_{\mathrm{cl},0}=10^{5}\,M_{\odot}$, for which the median stellar mass increases only mildly with central density at both metallicities. In these systems, the lower masses attained during the runaway phase result in weaker stellar winds, reducing the impact of wind-driven mass loss. Interestingly, for the most compact clusters with $R_{\mathrm{hm},0}=0.01$ pc and initial masses of $M_{\mathrm{cl},0}=10^{7}-10^{8}\,M_{\odot}$, we again find the formation of massive BHs in a non-negligible fraction of the population. This behavior is related to the interplay between the runaway timescale and our rejuvenation prescription. At such high densities, the VMS accretes the cluster core over a very short timescale, leaving limited time for subsequent rejuvenation. In contrast, in the more massive clusters at intermediate densities, stellar collisions can proceed over a longer period, allowing repeated rejuvenation episodes. Although the mass contributed by individual stars is small, the cumulative effect of successive collisions repeatedly rejuvenate the VMS, extending its effective lifetime and, consequently, the period over which it loses mass through stellar winds. This prolonged wind-driven mass loss ultimately lead to the formation of less massive BHs. We provide a more detailed assessment of the impact of our mass-loss prescriptions in Appendix F.

Vector (PDF)

In the densest clusters, newly formed BHs are expected to undergo further mass growth through subsequent dynamical interactions and accretion processes. To understand which mechanisms are likely to lead to the mass growth of BHs, it is important to assess how many massive stars are still available in the cluster. Figure 7 shows the distribution of the number of massive stars ($m_{*}\geq20\,M_{\odot}$) not accreted by the forming VMS. Less compact clusters generally have a larger number of massive stars, while more compact systems tend to deplete them earlier. In models where zero or only a few massive stars survive in the core, subsequent IMBH growth is expected to proceed mainly through tidal disruption events (TDEs) involving low-mass stars. We discuss these aspects in more detail in the next section.

Vector (PDF)

Discussion

In the following, we discuss the implications of IMBH seed formation for different cluster populations, and the efficiency of subsequent growth channels.

To do so, in Figure 8, we show the IMBH seed masses predicted by our model after 5 Myr for clusters with masses and radii consistent with those observed in GCs and young massive clusters in CGA, in order to assess whether these systems may have experienced conditions conducive to runaway collisions in the past. To this end, we evolved the mass and half-mass radius of our clusters up to 30 Myr and 10 Gyr. Some of our densest clusters evaporate before reaching 10 Gyr (i.e., those with $\log(M_{\mathrm{cl},0}/M_{\odot})=5$ and $R_{\mathrm{hm},0}<0.5$ pc, and $\log(M_{\mathrm{cl},0}/M_{\odot})=6$ and $R_{\mathrm{hm},0}<0.1$ pc), and are therefore not included in our analysis at late times.

Aside from CGAs and GCs, we also discuss implications for LRDs, relying on their inferred central densities and structural parameters to evaluate their capability to nurture IMBH seeds, as suggested by other studies (e.g. [127]).

Implications for Little Red Dots

Little Red Dots (LRDs) are extremely compact, red, high–surface-brightness sources identified at high and medium redshift ([100 Matthee et al.] (Kocevski et al. 2023) 2024), whose small sizes and dense stellar components suggest they may host intense central activity. Their origin remains highly uncertain, and several explanations have been proposed, ranging from gas-enshrouded AGN (Inayoshi & Maiolino [87]) to quasi-stars (Begelman & Dexter [26]), to self-gravitating disks accreting onto a supermassive star (Zwick et al. 2026).

Guia et al. [73] suggested that the inferred stellar masses of LRDs imply extremely high central densities, with a significant tail of the distribution reaching $10^{5}-10^{8}\,M_{\odot}\,\mathrm{pc}^{-3}$, thereby favoring the onset of runaway collisions. As shown in Figure 3, at such densities our models predict the formation of VMSs with masses above $10^{3}\,M_{\odot}$ for clusters with initial masses of $M_{\mathrm{cl},0}=10^{6},10^{7}\,M_{\odot}$, which in turn produce IMBHs with masses between $\sim10^{3}$ and $10^{4}\,M_{\odot}$.

However, the typical stellar mass associated with LRDs is $\gtrsim10^{8}\,M_{\odot}$ ([100, 5]). The case of $M_{\mathrm{cl},0}=10^{8}\,M_{\odot}$ and $R_{\mathrm{hm},0}=1\,\mathrm{pc}$ is particularly interesting, as the growth of the VMS is largely driven by the intrinsic stochasticity of dynamical interactions. Although the median VMS mass points towards $\sim200\,M_{\odot}$ (Table 2), the distribution is remarkably broad, with a mean mass of a few $\times10^{3}\,M_{\odot}$ (Figure 4, but see also Appendix § 10) and tails extending beyond $10^{4}\,M_{\odot}$, naturally resulting in the formation of IMBHs in a fraction of the realizations (Figure 6). At even higher densities, however, the system enters again a regime of extremely rapid VMS growth. In this regime, the VMSs quickly reach very high maximum masses, after which their evolution becomes dominated by strong stellar-wind mass loss, ultimately leading to the formation of less massive BHs.

Then, our results suggest that, despite the high stellar densities, LRDs can mainly sustain the growth of relatively light VMSs ($\sim10^{2}-10^{3}\,M_{\odot}$). Other works focused on the same topics, instead, found that LRD-like systems could lead to VMS up to $10^{3}-10^{4}\,M_{\odot}$ (e.g. [127, 181, 162]). We note that our cluster models inherently assume that cluster contraction is halted by VMS growth, as VMS accretion is halted after the mass in the core has been completely consumed by the VMS itself. Our approach is thus slightly different from the one followed by Pacucci et al. [127], which assumes a Fokker-Planck approach to resolve the differential collapse of the innermost layer of a $10^{11}\,M_{\odot}$ Plummer sphere. Interestingly, the approach of Pacucci et al. [127] leads to central clusters with around $10^{6}$ stars and stellar densities $\sim10^{7}\,M_{\odot}\,\mathrm{pc}^{-3}$, which sustain the VMS growth up to $10^{3}-10^{4}\,M_{\odot}$. Under the same assumption of masses and stellar densities, STARFALL is able to retrieve similar VMS masses, as shown in our Figure 4.

We also note that, although gas is likely to be present in these environments, we do not consider gas inflows in this work. The interplay between gas accretion and the inward transport of stellar mass may play a key role in the formation of central massive objects (e.g., [106], but see also [159, 148] (Beckholt et al. 2018)).

Implications for clusters in the Cosmic Gems Arc

The clusters in the CGA have typical masses of $\sim10^{6}-10^{7}\,M_{\odot}$ and half-mass radii up to a few parsecs, indicating very dense environments, with ages spanning $\sim5-40\,\mathrm{Myr}$ and metallicities between $Z\sim0.0004-0.001$ ([4, 116]).

As already shown in Figure 4, about 60% of models consistent with clusters in the CGA can host a VMS with mass $>1000\,M_{\odot}$. The formation of such massive stars in high-redshift systems may be connected to the nitrogen-rich chemical signatures observed in early galaxies, which have been interpreted as evidence of efficient CNO processing in VMSs ([38, 181] (Gilles et al. 2018; Gilles et al. 2025)). In this context, the D.1 ($\log(M_{\mathrm{cl}}/M_{\odot})=5.89$, $R_{\mathrm{hm}}=0.93\,\mathrm{pc}$) and E.1 ($\log(M_{\mathrm{cl}}/M_{\odot})=5.97$, $R_{\mathrm{hm}}=0.53\,\mathrm{pc}$) clusters are particularly favorable environments for the formation of IMBH seeds with masses close to or heavier than $1000\,M_{\odot}$ via runaway collisions. Using a Fokker-Planck approach, Bocchi et al. (2026) found that clusters with properties similar to the one inferred for those in the CGA can host IMBH seeds with masses in the range $1600-2700\,M_{\odot}$, consistent with our results. The remaining compatible models fall really in or below the PISN mass gap. This may point to the importance of additional formation channels for IMBHs. In addition, the CGA is a strongly lensed system ([4, 32, 116]). The counter-images represent the lensed projections of the same sources, providing an independent insight on its structural and photometric properties. If we consider the counter-images of these systems (orange squares in Figure 8), these systems do not produce an IMBH seed. However, the C.2+D.2 system is blended and is treated as a single cluster, which may affect its inferred parameters. A more accurate determination of its intrinsic size would help to better constrain its position in parameter space.

In addition, these systems are not expected to fully deplete their massive stellar population. This regime may therefore favor the coexistence of an IMBH and a residual population of massive stars, enabling additional growth through mass transfer or subsequent BBH mergers. In this case, the number of massive stars still available could also favor the formation of additional VMSs through dynamical interactions, which may then collapse into IMBHs, promoting the formation of IMBH–IMBH binaries. Note that upon merger, the remnant would receive a recoil that, given the typical mass-ratios, hardly exceeds 20 km s$^{-1}$, i.e. smaller than the typical escape velocity from such environments (see, e.g., Islam & Wadekar [89]). Interestingly, in a recent work, Torniamenti et al. [165] showed that clusters with properties similar to those in the CGA can actually form and retain IMBHs with masses exceeding $10^3\,M_\odot$ through hierarchical BH mergers.

Implications for globular clusters

Observed GCs with half-mass radii $R_{\rm hm} \lesssim5$ pc and masses $M_{\rm cl} \gtrsim10^{4.5}\,M_\odot$ represent promising candidates for hosting IMBH seeds formed via runaway collisions. Indeed, for a significant fraction of GCs, our models show a $\sim50\%$ probability of forming VMSs $\gtrsim300\,M_\odot$ (Figure 4). These stars have also been proposed as potential contributors to the anti-correlations of nitrogen or helium with respect to oxygen in local globular clusters ([71, 36]). However, the BH masses estimated for most of these systems ([110, 93, 75, 22]), and shown in Figure 8, are up to almost one order of magnitude larger than those predicted by our model. This unsurprising discrepancy may be due to the lack of a treatment to model stellar accretion onto the IMBH, and indicates that runaway collisions alone are unlikely to account for the full IMBH mass budget. Considering the IMBH masses in Lützgendorf et al. [110] and the ages of the respective host clusters ([37, 123, 169, 36] (Gontcharov et al. 2026)), we estimated a growth rate in the range $\sim(0.05 - 0.9)\,M_\odot/\mathrm{Myr}$. However, several studies (see, e.g., [97, 17, 167] (Arca sedda et al. 2024)) have shown that hierarchical BBH mergers in GCs can be suppressed by cluster expansion, thereby limiting mass growth through this channel. In addition, during the runaway phase, massive stars are preferentially merged onto the forming VMS, due to their higher interaction probability in the cluster core. While this process efficiently builds up the VMS mass, at the same time it reduces the number of massive stars available to evolve into BHs and participate in later dynamical interactions (see also Figure 7). In the lower-mass clusters considered here ($M_{\rm cl,0} = 10^5$–$10^6\,M_\odot$), the reservoir of massive stars in the core is limited to at most a few hundreds. In compact systems ($R_{\rm hm,0} \lesssim0.5$ pc), this reservoir can be depleted within a few thousand years, thereby limiting the efficiency of subsequent hierarchical BBH mergers. A similar behavior is also found in more massive clusters ($M_{\rm cl,0} = 10^7$–$10^8\,M_\odot$) with small half-mass radii ($R_{\rm hm,0} \lesssim0.1$ pc), potentially progenitors of the densest and most massive GCs observed today. Although these systems host a larger initial population of massive stars, their extreme densities favor rapid and efficient accretion onto the VMS, again depleting the reservoir of BH progenitors and suppressing the efficiency of hierarchical BBH mergers.

As a consequence, once an IMBH seed has formed, alternative growth channels are expected to become relevant. TDEs involving low-mass stars may dominate the early growth, potentially increasing the BH mass by $\sim20\%$ within a few Myr (e.g., [138], but see also [158, 141]). On longer timescales, comparable to the ages of observed GCs, TDEs involving low-mass compact remnants (e.g. white dwarfs) may contribute both to mass growth and to observable high-energy transients ([62, 161] (Arca sedda et al. 2024)). IMBH–white dwarf systems have been proposed as potential progenitors of long gamma-ray bursts ([126, 183]), making them especially relevant in the context of multi-messenger astrophysics and IMBH mass growth.

On the other hand, when considering massive clusters with larger half-mass radii ($R_{\rm hm,0} \sim0.5 - 1$ pc) that form an IMBH, several hundred or even thousands of massive stars may remain. Given the large number of stars available in the core, further dynamical interactions can lead to the formation of additional VMSs, which may in turn collapse into massive BHs or IMBHs. If these BHs are retained, such clusters provide favorable environments for the formation of IMBH–BH or IMBH–IMBH binaries. In such a case, relativistic recoils could become comparable to or larger than the cluster escape velocity, lowering the retention probability.

Uncertainties and limitations

The present implementation of STARFALL relies on some simplifying assumptions that may affect the generality of our results.

First, the adopted stellar evolution treatment neglects the effects of binary stellar evolution, such as mass transfer and common envelope phases. Although these processes can significantly alter the binary properties, they typically occur on timescales longer than those explored in our simulations, and are therefore expected to have a limited impact on our conclusions.

When two stars merge or collide, the radius of the remnant is computed as presented in Sect. 2.2, where $m_*$ is the sum of the masses of the stars involved in the interaction, corrected for mass loss. This approximation is partially supported by hydrodynamical simulations of stellar collisions ([63]); however, the actual radius may depend on the collision geometry, the internal structure of the progenitors, their masses, and how rapidly the star relaxes to its compact state. Indeed, after a collision the remnant could be initially expanded, and Figure 2 shows that the primary star can experience multiple collisions within just a few thousand years in dense environments. This timescale may be shorter than the time required for the star to contract before undergoing a new interaction (see, e.g., [41]). Indeed, when multiple collisions occur out of equilibrium, the resulting VMS remnant is likely to be puffed up, increasing its cross-section for contact with other stars ([45]). This may lead to an underestimation of the interaction rate. However, a larger cross-section does not necessarily translate into a higher merger rate, as intruding stars may pass through the puffy envelope without being captured or merging with the VMS. This may therefore lead to an underestimation of the interaction timescales. On the other hand, Verger et al. (2026) recently showed that using the cross-section argument to compute interaction timescales can result in an overestimation when a massive VMS interacts with a few stars left in the core. At each timestep, nevertheless, we limit the number of mergers and collisions to the number of objects within the VMS’s radius of influence. In addition, the prescription for mass loss adopted here for stellar collisions and mergers may break down in the densest and most massive clusters we considered. In such environments, the velocity dispersion can approach the stellar escape velocity, reducing the amount of mass that remains bound to the remnant.

As mentioned in Sect. 2.2, we determine the final fate of a star based on its pre-SN mass, although the PISN regime is primarily determined by the He-core properties. Most interactions in the collision chain shown in Figure 2 involve two main-sequence stars ($\sim99\%$, $\sim60\%$ when $R_{\rm hm} = 0.05, 0.1$ pc respectively), for which assuming a merger remnant structure similar to that of the corresponding rejuvenated star is a reasonable approximation. A different situation can arise when a post-main-sequence primary collides with a main-sequence star. In the densest clusters ($R_{\rm hm} = 0.05 - 0.1$ pc), this may happen when the primary star is already above the PISN mass gap ($m_* \sim$ few $\times10^3\,M_\odot$). In the less dense cluster ($R_{\rm hm} = 1$ pc), the post-main-sequence star may undergo only one or two collisions at lower masses and subsequently enter the PISN mass gap. To first order, the relevant quantity for determining the final fate in this case is the mass of the primary’s He core. Given our initial conditions, the primary has a core mass that leads to a PISN, allowing us to reasonably assume that the remnant will undergo the same fate. Interactions between two post-main-sequence stars with a collision remnant entering the pair-instability regime, which would require a more detailed treatment of the remnant’s internal structure, are rare ($\lesssim0.1\%$) and will be investigated in future work.

From a dynamical point of view, STARfall includes binary-single interactions and direct stellar collisions, neglecting binary-binary encounters that may become dominant in environments with a high binary fraction (e.g., [120, 184, 16]). In addition, the code does not account for the formation and evolution of hierarchical triple systems arising from dynamical interactions. Such systems can drive large-amplitude eccentricity oscillations in the inner binary through the Kozai-Lidov mechanism ([96, 105]), potentially leading to stellar collisions. We note, however, that for dense systems stable triples are very improbable ([103]).

Stellar masses are sampled assuming a Kroupa IMF ([99]) up to 150 M$_\odot$. While this choice is reasonable for Population I and II stars, it remains uncertain for Population III, for which both the shape and the limits of the mass function, as well as the properties of binary systems, are still poorly constrained ([34, 66, 108, 133]).

Finally, we note that our results should be regarded as conservative. Our simulations focus on stars initially located in the cluster core and neglect contributions from the outer regions. However, strong mass segregation episodes and primordial mass segregation can funnel more massive stars toward the center ([137]), potentially leading to a more massive VMS and increasing the number of BHs in the core, with implications for the subsequent growth of the central BH.

Summary and conclusions

In this work, we investigated the formation of VMSs and IMBH seeds through runaway collisions and stellar mergers in dense stellar clusters.

We developed STARfall, a semi-analytic tool designed to follow the evolution of a massive stellar binary in the cluster core, accounting for binary-single interactions, stellar mergers, and collisions, as well as the dynamical evolution of the host cluster, with particular focus on core collapse (Sect. 2).

We showed that in star clusters with masses in the range $M_{\mathrm{cl},0} = 10^{6} - 10^{8}\,\mathrm{M}_\odot$ and half-mass radii $R_{\mathrm{hm},0} < 1\,\mathrm{pc}$, stellar collisions are more likely to occur and can lead to the formation of VMSs with masses $m_* \gtrsim10^{3}\,\mathrm{M}_\odot$ (Figure 2, Figure 4). The maximum VMS masses are strongly correlated with the initial mass and density of the initial cluster (Figure 3). In addition, the cluster structure, its dynamical evolution, and the relevant stellar evolutionary timescales play a key role in determining the probability of VMS formation. Even in clusters with core-collapse timescales longer than 5 Myr, the probability of forming a VMS in our model remains significant (between 15% and 50%), provided that the cluster is massive enough ($M_{\mathrm{cl}} \gtrsim10^{7}\,\mathrm{M}_\odot$).

We compared the maximum stellar masses reached at the end of the collisional chains with results from $N$-body simulations presented by Arca Sedda et al. [17]; Vergara et al. [170]; Rantala et al. [138], and found overall good agreement, with STARfall also probing the high-density tail of the parameter space, which is not accessible to $N$-body simulations.

We also showed that cluster models with properties similar to those observed in the CGA ([4, 32, 116]) can host VMSs with $m_* \geq10^{3}\,\mathrm{M}_\odot$ in $\sim60\%$ of the cases (Figure 8), indicating that such environments are particularly favorable for the formation of both VMSs and IMBHs. However, uncertainties in the cluster sizes can significantly affect their position in parameter space, potentially shifting the predicted IMBH seed masses toward, or even below, the PISN mass gap.

Globular clusters (GCs) also represent favorable environments for forming VMSs with $m_* \geq300\,\mathrm{M}_\odot$ and IMBH seeds (Figure 8; [110, 77, 93]). In dense, low-mass GCs, where the population of massive stars has already been largely depleted (Figure 7), alternative channels to BH-BH mergers may become dominant for subsequent mass growth of the IMBH. In particular, TDEs involving low-mass stars ([158, 141, 138]) are expected to play a key role in driving further mass growth.

In contrast, although LRDs may have high stellar masses and densities ([73, 5]), their structural properties, if interpreted as stellar clusters, do not always favor the onset of long stellar collision chains. For core densities between $10^{5}$ and $10^{8}\,\mathrm{M}_\odot\,\mathrm{pc}^{-3}$ and cluster masses of $10^{6}$–$10^{7}\,\mathrm{M}_\odot$, IMBHs with masses in the range $10^{3}$–$10^{4}\,\mathrm{M}_\odot$ can form. For more massive clusters ($M_{\mathrm{cl}} = 10^{8}\,\mathrm{M}_\odot$), the longer core-collapse timescale results in a lower median IMBH mass. Although the distribution is very broad, the probability of forming BHs with masses above $10^{3}\,\mathrm{M}_\odot$ remains low (Sect. 4.1). This suggests that high density alone is not sufficient to trigger collision chains. Resolving their internal structure will play a key role in determining their dynamical evolution. However, Pacucci et al. [127] suggest that, in these high density systems, mass segregation happens on very short timescales favoring the onset of runaway collisions and the formation of more massive VMSs. Nevertheless, when considering the same initial conditions ($M_{\mathrm{cl},0} \sim10^{6}\,\mathrm{M}_\odot$, $R_{\mathrm{hm},0} \sim0.1\,\mathrm{pc}$), we find consistent results.

We further deepened the impact of metallicity and stellar wind mass loss on the mass of IMBH seeds (Figure 5,Figure 6; but see also Appendix § 12). We found that, although denser clusters tend to produce more massive VMSs, they do not necessarily produce the most massive BHs, with clusters at intermediate densities ($\rho_{0,\mathrm{core}} \sim10^{7}$–$10^{9}\,\mathrm{M}_\odot\,\mathrm{pc}^{-3}$) producing IMBHs with masses higher than those produced in denser environments. This introduces a dependence on both the formation timescale of the VMS and its subsequent evolution (see also Figure 2).

While our work identifies the regimes in which VMSs with masses above $10^{3}\,\mathrm{M}_\odot$ can form, it also highlights the critical role of the subsequent stellar evolution once these objects reach their maximum mass. Future studies that accurately model the evolution of collision products, together with improved prescriptions for stellar winds, will be essential to further refine these results.

Acknowledgements. We thank the anonymous referee for their careful reading of the manuscript and their valuable comments, which improved the clarity and quality of this work. The authors gratefully acknowledge Antti Rantala, Michela Mapelli, Filippo Santoliquido, Gor Oganesyan and Dominik Schleicher for the helpful comments and suggestions. 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) and from the MERAC Foundation through the 2023 MERAC prize. MAS and MB acknowledge the ACME project which has received funding from the European Union’s Horizon Europe Research and Innovation programme under Grant Agreement No. 101131928. GI is supported by a fellowship grant from la Caixa Foundation (ID 100010434). The fellowship code is LCF/BQ/PI24/12040020. PB thank the support from the special program of the Polish Academy of Sciences and the U.S. National Academy of Sciences under the Long-term program to support Ukrainian research teams, grant No. PAN.BFB.S.BWZ.329.022.2023. AA acknowledges that this research was funded in part by National Science Centre (NCN), Poland, grant No. 2024/55/D/ST9/02585. For the purpose of Open Access, the author has applied a CC BY public copyright licence to any Author Accepted Manuscript (AAM) version arising from this submission. PB also appreciate the Polish high-performance computing infrastructure PLGrid (HPC Centre: ACK Cyfronet AGH – docs.hpc.cyfronet.pl) for providing Helios computer facilities and support within computational grant No. PLG/2026/019243. GC acknowledges financial support from the European Union–Next Generation EU, Mission 4, Component 2, CUP: C93C24004920006, project ‘FIRES’. MG was supported by the Polish National Science Center (NCN) through the grant 2021/41/B/ST9/01191. ST acknowledges financial support from the Alexander von Humboldt Foundation for the Humboldt Research Fellowship. MCV acknowledges financial support from the Center for Astrophysics and Associated Technologies (CATA; FB210003). We used B-POP (https://github.com/marcasedda/BPOP) to model star cluster evolution (Arca Sedda et al. [11, 17] (2026)). We used sevn (https://gitlab.com/sevncodes/sevn; Spera et al. [157]; Iorio et al. [88]) to interpolate Parsec stellar tracks (Costa et al. [44]). This research made use of NumPy (Harris et al. [76]), SciPy (Virtanen et al. [176]), Pandas (The pandas development Team [163]). For the plots we used Matplotlib (Hunter [83]).

Data availability

The STARfall code will be publicly available soon. In the mean-time, the code can be made available upon request. Interested users are encouraged to contact the authors.

Star cluster model

In STARfall, star clusters are modeled following the prescriptions in the B-POP code ([11, 17] (Arca Sedda et al. 2026)). Starting from a cluster mass $M_{\mathrm{cl},0}$, the evolution of the bound mass is described by:

where $t_{\mathrm{st}}=10\,\mathrm{Myr}$, $k=-0.1$, and $t_{\mathrm{rlx}}$ is the relaxation time. The parameter:

regulates the onset of relaxation-driven mass loss. We adopted $\alpha_{\mathrm{bnd}}=1000$ and $\beta_{\mathrm{bnd}}=0.6$. These parameters are calibrated to reproduce globular cluster evolution reported by Gieles & Gnedin [60].

Similarly, given a half-mass radius $R_{\mathrm{hm},0}$, its time evolution follows:

where $t_{0}=5\,t_{\mathrm{cc}}\left[1+\left(R_{\mathrm{hm},0}/1\,\mathrm{pc}\right)\right]$, $t_{\mathrm{cc}}$ is the core-collapse time, and $\alpha_{\mathrm{hm}}$ is randomly sampled at each time step from a uniform distribution between $0.15-0.25$. This prescription reproduces the half-mass radius evolution observed in the Dragon-II simulations ([15] (Arca Sedda et al. 2024; Arca sedda et al. 2024)).

The core properties are modeled through semi-analytical relations linking the core mass and core radius to $M_{\mathrm{cl},0}$:

where $\alpha_{\mathrm{core}}=0.03$. The exponent $\beta_{\mathrm{core}}$ depends on the evolutionary stage of the cluster and is set to $\beta_{\mathrm{core}}=0.4$ for $t<2t_{\mathrm{cc}}$, and to $\beta_{\mathrm{core}}=1/1.3$ at later times.

The ratio between the core radius and the half-mass radius evolves according to:

where $R_{\min}\propto(m_{\mathrm{mean}}/M_{\mathrm{cl}})^{1.3}$, and $m_{\mathrm{mean}}$ is the mean stellar mass ([30]).

Binary-single interactions

Binary-single interaction modeling is based on Heggie [79]; Hut & Bahcall [86]; Heggie & Hut [80]; Sigurdsson & Phinney [153].

The outcome of a binary–single encounter is set by the ratio of the transferred energy:

to the perturber kinetic energy in the three-body center-of-mass frame:

where $a$ is the binary semi-major axis, $r_{\mathrm{p,e}}$ the minimum separation between the binary and the incoming perturber, $m_{\mathrm{bin}}=m_{1}+m_{2}$, $m_{\mathrm{p}}$ the perturber mass, and $v_{\infty}$ its velocity at infinity ([79, 153]). If the transferred energy is small ($\Delta E_{\mathrm{bs}}<\Delta E_{\mathrm{c}}$), the binary outcome is set by its binding energy $E_{\mathrm{bin}}$. Specifically, if $\Delta E_{\mathrm{bs}<E_{\mathrm{bin}}}$, the binary either hardens or softens depending on the local velocity dispersion. Otherwise, the binary undergoes an exchange, replacing its least massive component if the perturber is heavier than the binary or its individual members.

If the binary preserves its components, the perturber will recede to infinity with velocity $v_{\mathrm{p}}<v_{\infty}$. If not, the former binary component with mass $m_{i}$ will escape with velocity

High-energy encounters ($\Delta E_{\mathrm{bs}}>\Delta E_{\mathrm{c}}$) are treated statistically by comparing $v_{\infty}$ to the critical velocity for ionization $v_{\mathrm{c}}$ ([86]):

For slow interactions ($v_{\infty}<v_{\mathrm{c}}$), the three bodies undergo a resonant encounter, typically resulting in an exchange of a binary component if the perturber is sufficiently massive or the energy transfer is large; otherwise, the lighter component is ejected. In fast interactions ($v_{\infty}>v_{\mathrm{c}}$), the binary either swaps a component if $\Delta E_{\mathrm{bs}}<E_{\mathrm{bin}}$ or is ionized if $\Delta E_{\mathrm{bs}}>E_{\mathrm{bin}}$.

The scattering is assumed to occur at $v_{\infty}=\sigma$, equal to the cluster velocity dispersion, which is a time-varying function that depends on the cluster mass and radius.

The interaction is assumed to be hyperbolic and fully characterized by the impact parameter and pericenter distance. The squared impact parameter, $b^{2}$, is drawn from a flat distribution between $0$ and $b_{\max}^{2}$, where $b_{\max}$ is defined as:

with $C=3$ (see e.g. Sigurdsson & Phinney [153], and reference therein). The pericenter, instead, is given by:

and the energy transferred in the interaction is:

The quantities $A$ and $B$ are functions of the binary orbital properties and scattering parameters (see Eqs. 2.12-2.13 in Sigurdsson & Phinney [153]). However, Eq. B.6 is valid only in the limit $r_{\mathrm{p}}/a_{\mathrm{bin}}\gg1$. Since, by construction, we consider hard binaries (see Sect. 2.3), for encounters that do not satisfy this condition, we instead sample the fractional energy change, $x\equiv\Delta E/E_{\mathrm{bin}}$, from the distribution derived by Heggie [79], who showed that the rate of close encounters producing a fractional change $x$ in the binary binding energy scales as $(1+x)^{-9/2}$.

If the recoil velocity imparted on the binary after the interaction exceeds the cluster escape velocity, $v_{\mathrm{esc}}\simeq(2GM_{\mathrm{cl}}/R_{\mathrm{hm}})^{1/2}$, the binary is ejected from the cluster.

Collision outcomes

In this Appendix we examine in more detail the assumption that a contact between two stars results in a merger. We note that, when computing the collision timescale, $t_{\mathrm{coll}}$, we always assume the pericenter distance, $r_{\mathrm{p}}$, to be equal to the sum of the radii of the two stars, i.e. $r_{\mathrm{p}}/(R_{1}+R_{2})=1$, where $R_i$ are the radius of the two star (see Sect. 2.3). Then, following Freitag & Benz [53], we define:

where $m_i$ are the masses of the two stars, and $v_{\mathrm{rel}}^{\infty}\sim\sigma$, with $\sigma$ the velocity dispersion of the cluster. In Fig. C.1 we show the distribution of $v_{\mathrm{rel}}^{\infty}/v_{\mathrm{contact}}$ as a function of the mass ratio for all collisions recorded in the simulated systems, for two representative cluster masses ($M_{\mathrm{cl},0}=10^{6}$ and $10^{8}\,M_{\odot}$) with $R_{\mathrm{hm},0}=0.5$ pc. We select these two cases because for $M_{\mathrm{cl},0}=10^{6}\,M_{\odot}$ the bulk of collisions is expected to occur during core collapse, while the more massive counterpart allows us to assess how the collision properties vary with cluster mass.

The vast majority of events concentrate at $v_{\mathrm{rel}}^{\infty}/v_{\mathrm{contact}}\lesssim0.2$, indicating that the relative velocity at infinity is small compared to the contact velocity and that the encounters are gravitationally focused. Furthermore, after the merger of the primordial binary (the necessary first step to initiate the collision chain; see also Sect. 3), the mass ratio rapidly decreases, driving the system into a regime where, following Freitag & Benz [53], a merger outcome is well justified. This is particularly relevant for $M_{\mathrm{cl},0}=10^{8}\,M_{\odot}$, where a non-negligible fraction of collisions occurs at $v_{\mathrm{rel}}^{\infty}/v_{\mathrm{contact}}\sim0.5-0.6$; however, these events have mass ratios below $0.05$, which still supports the merger assumption. Only a few isolated events fall in a region of parameter space where the distinction between a fly-by and a merger becomes ambiguous.

Finally, we note that, with respect to Freitag & Benz [53], we do not resolve the stellar structure, and in Eq. C.1 we therefore adopt the total stellar radius rather than the half-mass radius, which leads to an overestimation of $v_{\mathrm{rel}}^{\infty}/v_{\mathrm{contact}}$. In addition, as discussed in Sect. 2.1, core collapse is intrinsically halted by the presence of a VMS (or a binary system), preventing the velocity dispersion from reaching very high values. Moreover, after each collision the VMS radius is recomputed according to Eq. (3), yielding a relatively compact remnant. Taken together, these considerations justify treating contact events as mergers.

Mass distribution of very massive stars and intermediate-mass black hole seeds

Figures D.2 and D.3 show the distributions of the maximum stellar mass and IMBH seed mass for clusters with $M_{\mathrm{cl},0}=10^{7}-10^{8}\,M_{\odot}$. The overall trends are similar to those found for $M_{\mathrm{cl},0}=10^{6}\,M_{\odot}$ (Figure 5), but shifted toward higher masses.

For $M_{\mathrm{cl},0}=10^{7}\,M_{\odot}$ ($10^{8}\,M_{\odot}$) and $R_{\mathrm{hm},0}=0.05,0.1$ pc, the stellar mass distribution exhibits a peak at $\sim10^{5}\,M_{\odot}$ ($\sim6\times10^{5}\,M_{\odot}$), nearly independent of metallicity ($Z=0.0002,0.002$). At $R_{\mathrm{hm},0}=0.5$ pc, nearly half of the realizations produce VMSs with masses exceeding $10^{4}\,M_{\odot}$ for both cluster masses. At $R_{\mathrm{hm},0}=1$ pc, however, the two cluster masses exhibit markedly different distributions. For $M_{\mathrm{cl},0}=10^{7}\,M_{\odot}$, the VMS mass distribution is entirely confined between $10^{2}$ and $10^{3}\,M_{\odot}$. At the higher cluster mass, instead, the distribution becomes substantially broader, developing a prominent high-mass tail that extends beyond $10^{4}\,M_{\odot}$.

For clusters with initial mass $M_{\mathrm{cl},0}=10^{7}\,M_{\odot}$, the most massive IMBHs reach $10^{5}\,M_{\odot}$ for $R_{\mathrm{hm},0}=0.05$ pc and $\sim5\times10^{4}\,M_{\odot}$ for $R_{\mathrm{hm},0}=0.5$ pc, with higher metallicity reducing BH masses by nearly an order of magnitude. For more massive clusters, instead, $M_{\mathrm{cl},0}=10^{8}\,M_{\odot}$, the most massive IMBHs form at $R_{\mathrm{hm},0}=1$ pc, although only in a few realizations. In this case, the densest clusters do not produce BHs more massive than $\sim10^{4}\,M_{\odot}$, regardless of metallicity.

For lower-mass clusters ($M_{\mathrm{cl},0}=10^{5}\,M_{\odot}$, Fig. D.1), this bimodal behavior is not observed. The stellar mass distribution peaks at $\sim(5-7)\times10^{2}\,M_{\odot}$ for $R_{\mathrm{hm},0}=0.1$ and $0.5$ pc, whereas it extends to higher masses, reaching $\sim2\times10^{3}\,M_{\odot}$, for $R_{\mathrm{hm},0}=0.05$ pc. Accordingly, the IMBH seed mass can range from $\sim500\,M_{\odot}$ to $\sim2000\,M_{\odot}$, depending on the initial cluster density. For $R_{\mathrm{hm},0}=1$ pc, the smoother dynamical evolution largely suppresses runaway collisions, with only a few exceptions.

Impact of stellar evolution prescriptions

As discussed in Sect. 2, we extrapolated fitting formulae calibrated on stellar evolution tracks from SEVN-Parsec ([157, 88, 44]) to describe the initial stellar radius ($R_{\mathrm{ZAMS}}$; Eq. (3)), the time at which the star begins to expand ($\tau_{\mathrm{MS}}$; Eq. (4)) and the size of this expansion ($f_{\mathrm{rad}}=R_{\ast}(\tau_{\mathrm{MS}})/R_{\mathrm{ZAMS}}$; lower panel of Fig. E.1), the stellar lifetime ($\tau_{\mathrm{life}}$; Eq. (6)), and the remnant mass (upper panel of Fig. E.1).

We assumed a constant stellar radius, $R_{\mathrm{ZAMS}}$, for $t \leq\tau_{\mathrm{MS}}$. After this time, it is increased by a factor $f_{\mathrm{rad}}$. As noted when discussing Figure 2, we do not expect a more gradual growth of the stellar radius to significantly affect the final stellar mass, but rather to produce a smoother rise of the mass.

We computed the mass of the compact remnants extrapolating the relation between the pre-SN mass and the remnant mass, i.e., assuming that the remnant mass depends only on the supernova phase (rapid formalism; Fryer et al. [56]). This can have an important consequence on the final fate of the VMS. Within the range of initial stellar masses considered for the primary component of the binary (uniform distribution between $100$–$150\,M_{\odot}$), a large fraction of stars is expected to fall into the pair-instability regime (Woosley [182]; Spera & Mapelli [156]). Under these conditions, the evolution may lead to the formation of low-mass black holes in the case of a PPISN, or to no compact remnant at all in the case of a PISN, despite the high mass of the progenitor stars.

At $Z=0.0002$ ($Z=0.002$), the fraction of stars entering the PISN regime increases from $1\%$ ($4\%$) to $38\%$ ($63\%$) as the half-mass radius varies from $R_{\mathrm{hm},0}=0.05$ pc to $R_{\mathrm{hm},0}=1$ pc.

The fraction of stars producing stellar-mass black holes varies, instead, from $\sim3\%$ to $54\%$ at $Z = 0.0002$ for $R_{\mathrm{hm},0} = 0.05$ pc and $R_{\mathrm{hm},0} = 1$ pc respectively. For the case $Z = 0.002$, no stellar-mass BH are formed. All these values refer to clusters with an initial mass $M_{\mathrm{cl},0} = 10^6\,M_\odot$ (Figure 2; 5;E.2).

This behavior results from the combined effect of our extrapolation of the sevn-Parsec stellar tracks, initial conditions, and the mass loss induced by stellar winds and stellar collisions that allows stars to enter or exit the PISN regime depending on the number of interactions and their lifetime.

These values should not be interpreted as absolute rates, as they depend on the adopted initial conditions.

Impact of mass loss

In our work, we considered mass loss driven by both stellar winds and collision.

To model mass loss via stellar winds, we adopted the prescription by Vink [174]:

The mass loss in stellar collision and mergers is parametrized as in Glebbeek et al. [64]:

with $q = m_2/m_1 < 1$ and $D = 0.3$.

In addition, we ran some specific simulations to assess the effect of our assumptions.

Figure F.1 shows the distribution of the maximum primary mass $m_{1,\max}$ obtained in clusters with $M_{\mathrm{cl},0} = 10^5, 10^8 M_\odot$ for different metallicities ($Z = 0.0002, 0.002$) at $R_{\mathrm{hm},0} = 0.5$ pc.

Three different prescriptions for mass loss are compared: stellar wind and collision-driven mass loss, mass loss via stellar winds only, no mass loss.

The comparison highlights the strong impact of mass loss on the maximum stellar masses reached in the simulations. When no mass loss is included, the runaway growth driven by stellar collisions can produce extremely massive stars, reaching $> 10^6M_\odot$ in the most massive clusters (lower panels). Including stellar winds only already reduces the high-mass tail of the distribution, limiting the growth of the most massive objects. When both winds and collisional mass loss are taken into account, the maximum masses are further suppressed, shifting the distributions to significantly lower values ($\lesssim10^5 M_\odot$).

This effect becomes particularly evident in the most massive clusters ($M_{\mathrm{cl},0} = 10^8 M_\odot$), where runaway growth would otherwise lead to the formation of extremely massive stars in the absence of mass loss. The combined effect of winds and collision-induced mass loss efficiently regulates the growth of the most massive objects, preventing the formation of the most extreme masses predicted in the no-mass-loss scenario. The impact of stellar winds is metallicity dependent, with higher metallicity models ($Z = 0.002$) showing a stronger suppression of the high-mass tail due to the increased efficiency of wind-driven mass loss.

Vector (PDF)