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

abstract · abstract

Intermediate-mass black holes (IMBHs, $10^2 < m_{\rm BH}/\mathrm{M}_{\odot} < 10^5$) are elusive objects bridging the gap between stellar-mass and supermassive black holes, whose origin remains uncertain. They may form from the collapse of very massive stars (VMSs, $m_*\geq150\,\mathrm{M}_{\odot}$) assembled through repeated stellar collisions in compact and fast evolving star clusters. We investigate the formation of VMSs and IMBH seeds via repeated mergers and collisions as a function of cluster initial conditions, assessing the implications for globular clusters and dense high-redshift system. We develop STARfall, a semi-analytic tool that follows the evolution of a massive stellar binary via binary-single interactions, mergers, collisions, and stellar evolution, while consistently accounting for the host cluster evolution. Our results show that clusters with central densities above $10^6\,\mathrm{M}_{\odot}\,\mathrm{pc}^{-3}$ can undergo a catastrophic sequence of stellar collisions, or "starastrophe", within $3$ Myr. Through this process, VMSs rapidly grow in mass and collapse into IMBHs with masses of up to $10^3-10^4\,\mathrm{M}_{\odot}$, despite the effects of intense stellar winds. Applying our models to observed systems, we find that compact clusters in the Cosmic Gems Arc and Milky Way globular clusters could have hosted VMSs that collapsed into IMBHs in this mass range. Despite their longer core-collapse timescales, the inferred central densities of Little Red Dots may also support IMBH formation in the same mass interval. Our models show that VMS assembly requires two essential ingredients: short core-collapse times and large central densities. Compared to observed systems, our models suggest that clusters in the Cosmic Gems Arcs and several Galactic globular clusters could have nurtured VMSs that directly collapsed to IMBHs with masses $10^3-10^4\,\text{M}_{\odot}$.

S5: 5 Summary and conclusions / block 1 · html

We showed that in star clusters with masses in the range M cl , 0 = 10 6 − 10 8 ​ M ⊙ M_{\rm cl,0}=10^{6}-10^{8}\,{\rm M}_{\odot} and half-mass radii R hm , 0 < 1 R_{\rm hm,0}<1 pc, stellar collisions are more likely to occur and can lead to the formation of VMSs with masses m ∗ ≳ 10 3 ​ M ⊙ m_{*}\gtrsim 10^{3}\,{\rm M}_{\odot} (Figs. 2 , 4 ). The maximum VMS masses are strongly correlated with the initial mass and density of the host cluster (Fig. 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 cl ≳ 10 7 ​ M ⊙ M_{\rm cl}\gtrsim 10^{7}\,{\rm M}_{\odot} ).

S5: 5 Summary and conclusions / block 1 · html

We also showed that cluster models with properties similar to those observed in the CGA ( Adamo et al., 2024 ; Bradley et al., 2025 ; Messa et al., 2026 ) can host VMSs with m ∗ ≥ 10 3 ​ M ⊙ m_{*}\geq 10^{3}\,{\rm M}_{\odot} in ∼ 60 % \sim 60\% of the cases (Fig. 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.

S4.SS4: 4.4 Uncertainties and limitations / block 1 · html

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 ( Rantala et al., 2024 ) , 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.

S5: 5 Summary and conclusions / block 1 · html

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

S5: 5 Summary and conclusions / block 1 · html

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 ).

S5: 5 Summary and conclusions / block 1 · html

We compared the maximum stellar masses reached at the end of the collisional chains with results from N N -body simulations presented by Arca Sedda et al. 2023c ; Vergara et al. 2026 ; Rantala et al. 2026 ; Askar et al. in prep.; Giersz et al., in prep. (Fig. 3 ). We find overall good agreement, with STARfall also probing the high-density tail of the parameter space, which is not accessible to N N -body simulations.

S5: 5 Summary and conclusions / block 1 · html

Globular clusters (GCs) also represent favorable environments for forming VMSs with m ∗ ≥ 300 ​ M ⊙ m_{*}\geq 300\,{\rm M}_{\odot} and IMBH seeds (Fig. 8 ; Lützgendorf et al. 2013 ; Harris et al. 2015 ; Kızıltan et al. 2017 ). In dense, low-mass GCs, where the population of massive stars has already been largely depleted (Fig. 7 ), alternative channels to BH-BH mergers may become dominant for subsequent mass growth of the IMBH. In particular, TDEs involving low-mass stars ( Stone et al., 2017 ; Rizzuto et al., 2023 ; Rantala et al., 2026 ) are expected to play a key role in driving further mass growth.

S5: 5 Summary and conclusions / block 1 · html

In contrast, although LRDs may have high stellar masses and densities ( Guia et al., 2024 ; Akins et al., 2025 ) , 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 10^{5} and 10 8 ​ M ⊙ ​ pc − 3 10^{8}\,{\rm M}_{\odot}\,\text{pc}^{-3} and cluster masses of 10 6 10^{6} – 10 7 ​ M ⊙ 10^{7}\,{\rm M}_{\odot} , IMBHs with masses in the range 10 3 10^{3} – 10 4 ​ M ⊙ 10^{4}\,{\rm M}_{\odot} can form. For more massive clusters ( M cl = 10 8 ​ M ⊙ M_{\rm cl}=10^{8}\,{\rm 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 ​ M ⊙ 10^{3}\,{\rm 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. (2025) 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 cl , 0 ∼ 10 6 ​ M ⊙ M_{\rm cl,0}\sim 10^{6}\,{\rm M}_{\odot} , R hm , 0 ∼ 0.1 R_{\rm hm,0}\sim 0.1 pc), we find consistent results.

S5: 5 Summary and conclusions / block 1 · html

We further deepened the impact of metallicity and stellar wind mass loss on the mass of IMBH seeds (Figs. 5 , 6 ; but see also Appendix F ). 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 ( ρ 0 , core ∼ 10 7 − 10 9 ​ M ⊙ ​ pc − 3 \rho_{\rm 0,core}\sim 10^{7}-10^{9}\,{\rm M}_{\odot}\,{\rm 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 Fig. 2 ).

S2.SS2: II.2 Analytical Calculation / block 1 · html

for f c = 0.8 f_{c}=0.8 , which is the empirical result found in [ 26 ] . The mass accretion rate onto the central VMS is estimated as ∼ 3 % \sim 3\% of the star formation rate in the stellar cluster. This result holds across a wide range of stellar densities and metallicities and is the outcome of mass segregation and binary-induced collisions in the evolution of a compact cluster.

S3.SS2: III.2 Analytical Calculation: Results / block 1 · html

To perform the calculation, we need to estimate the star formation rate and the metallicity of a typical LRD. Regarding the former, estimates for the LRDs vary widely. For example, [ 36 ] estimate values between ∼ 12 ​ M ⊙ ​ yr − 1 \sim 12\,\,{\rm M_{\odot}}\,\mathrm{yr}^{-1} and ∼ 49 ​ M ⊙ ​ yr − 1 \sim 49\,\,{\rm M_{\odot}}\,\mathrm{yr}^{-1} . [ 75 ] obtain values between ∼ 11 ​ M ⊙ ​ yr − 1 \sim 11\,\,{\rm M_{\odot}}\,\mathrm{yr}^{-1} and ∼ 30 ​ M ⊙ ​ yr − 1 \sim 30\,\,{\rm M_{\odot}}\,\mathrm{yr}^{-1} , although instantaneous values may reach ∼ 150 ​ M ⊙ ​ yr − 1 \sim 150\,\,{\rm M_{\odot}}\,\mathrm{yr}^{-1} . [ 43 ] infer values in an ultraluminous LRD at z = 4.47 z=4.47 that are significantly higher. Given such a range of estimates, we assume a conservative value of ∼ 30 ​ M ⊙ ​ yr − 1 \sim 30\,\,{\rm M_{\odot}}\,\mathrm{yr}^{-1} . Regarding the metallicity, estimates from [ 36 ] and [ 75 ] provide values between 6 % ​ Z ⊙ 6\%\,\rm Z_{\odot} and 32 % ​ Z ⊙ 32\%\,\rm Z_{\odot} , while [ 43 ] provides a value of 10 % ​ Z ⊙ 10\%\,\rm Z_{\odot} . Hence, we reasonably assume a reference value of Z = 10 % ​ Z ⊙ Z=10\%\,\rm Z_{\odot} . From the star formation rate, we calculate a constant accretion rate of M ˙ acc ∼ 1 ​ M ⊙ ​ yr − 1 \dot{M}_{\rm acc}\sim 1\,\,{\rm M_{\odot}}\,\mathrm{yr}^{-1} .

S2.SS3: II.3 Collisions / block 1 · html

We integrate this expression analytically and provide the full formulae in Appendix A . For black hole systems, the equilibrium condition results in the [ 7 ] stellar power law of ρ ∝ r − α \rho\propto r^{-\alpha} with α = 7 / 4 \alpha=7/4 [ 55 , 78 , e.g.,] . Given that we also consider self-gravitating stellar systems, and systems which may be out of equilibrium, we allow for a range of power law indices between 1-2. Figures are presented for α = 6 / 5 \alpha=6/5 , a more cored profile than the [ 7 ] profile, which will result in a more conservative estimate (as demonstrated by the lower left panel of Figure 3 ). We discuss the effect of varying α \alpha in the text.

S2.SS5: II.5 Collision outcomes / block 1 · html

where v esc , ∗ v_{\rm esc,*} is the escape velocity from the surface of the stars, then the collision will result in significant mass loss–typically destroying the stars involved. Such collisions release material of the inner layers of the stars into the local region.

S1: 1 Introduction / block 1 · html

Astrophysical models in which a compact object such as a neutron star (NS) or a black hole (BH) is embedded within the gaseous layers of a star have for a long time been constructed out of theoretical curiosity and invoked to interpret various unexplained observational results (e.g. Clayton et al. 1975 ; Oncins et al. 2022 ; Bellinger et al. 2023 ) especially in the context of binary stellar evolution (e.g. Sparks and Stecher 1974 ; Paczynski 1976 ; MacLeod and Ramirez-Ruiz 2015 ; Cruz-Osorio and Rezzolla 2020 ) and supermassive black hole (SMBH) seed formation. Such more or less exotic astrophysical objects include red giant stars with NS cores (Thorne-Żytkow objects, Thorne and Zytkow 1977 ; Cannon et al. 1992 ; Hutilukejiang et al. 2018 ; Farmer et al. 2023 ; Everson et al. 2024 ; Nathaniel et al. 2025 ; Williams et al. 2025 ) or quasi-stars (QSs; Begelman et al. 2008 ; Volonteri and Begelman 2010 ; Ball et al. 2011 ; Schleicher et al. 2013 ; Coughlin and Begelman 2024 ) in hydrostatic equilibrium with their central BH cores. Quasi-stars have been considered as relatively short-lived black hole gas envelope phases of the low metallicity direct collapse gas cloud scenario channel (DCBH; Bromm and Loeb 2003 ; Begelman et al. 2006 ; Dijkstra et al. 2008 ; Begelman 2010 ; Agarwal et al. 2012 ; Latif et al. 2013 ; Habouzit et al. 2016 ; Regan et al. 2017 ; Wise et al. 2019 ; Chon et al. 2021 ; Jeon et al. 2025 ) for SMBH seed formation ( Rees, 1978 ; Volonteri et al., 2008 ; Natarajan, 2014 ; Inayoshi et al., 2020 ) . Recently, attention for such embedded BH scenarios has been reinvigorated by the discovery of a novel and enigmatic class of objects now known as little red dots (LRDs; Matthee et al. 2024 ) by the James Webb Space Telescope ( JWST ). LRDs are characterized, among other properties, by their abundance at z ≳ 4 z\gtrsim 4 ( Kokorev et al., 2024 ; Kocevski et al., 2025 ) , compact ( ≲ 100 \lesssim 100 pc) sizes (e.g. Labbe et al. 2025 ), and peculiar spectral energy distributions (SEDs; Barro et al. 2024 ; Harikane et al. 2023 ; Williams et al. 2024 ; Yue et al. 2024 ; D’Eugenio et al. 2026 ; Inayoshi et al. 2026 ) including broad emission lines, X-ray

S4: 4 Discussion: the formation and the early evolution of faint multiple LRD systems / block 1 · html

Observationally, the expected bolometric luminosity of a ∼ 10 4 ​ M ⊙ \sim{10^{4}}\>\mathrm{M_{\odot}} SMS or QS is of the order of L QS ∼ 10 8 ​ L ⊙ L_{\mathrm{QS}}\sim 10^{8}L_{\odot} – 10 9 ​ L ⊙ 10^{9}L_{\odot} ( Santarelli et al., 2026 ; Roman-Garza et al., 2026b ) corresponding to log 10 ⁡ ( L QS / erg s -1 ) ∼ 42.5 \log_{\mathrm{10}}(L_{\mathrm{QS}}/\textrm{erg\>s${}^{-1}$})\sim 42.5 – 43.5 43.5 . As such, the QS like SMSs in our simulations are less luminous than typical LRDs with log 10 ⁡ ( L / erg s -1 ) ∼ 44 \log_{\mathrm{10}}(L/\textrm{erg\>s${}^{-1}$})\sim 44 ( L ∼ 3 × 10 10 ​ L ⊙ L\sim 3\times 10^{10}L_{\odot} ) (e.g. Akins et al. 2025 ; Greene et al. 2026 ; Umeda et al. 2026 ). However, the inferred QS luminosities are consistent with luminosities of the potential emerging population of faint LRDs ( Yanagisawa et al., 2026 ) that can be detected in strongly magnified star forming galaxies at high redshift.

S4.SS1: IV.1 SEDs of Individual Component / block 1 · html

Figures 3 , 4 , and 5 , show the SED-fitting results for each component. The photometric redshifts of the individual components are generally consistent at z ∼ 7.4 z\sim 7.4 . The BAGPIPES fits that include a blackbody component reproduce both the blue UV continuum and the red optical continuum observed in the LRD images. Without this component, the best-fit solutions require dusty star-forming galaxy templates to match the red optical continuum, which in turn leads to a significant underestimate of the UV flux. In contrast, the SFG components are well-explained by young star-forming galaxy models. We note that the fits to the SFG images with and without the blackbody component yield nearly identical best-fit SEDs, implying little contribution from the blackbody component. For the LRD images, the best-fit T BB T_{\mathrm{BB}} values are ∼ 4000 ​ – ​ 4500 ​ K \sim 4000\text{--}4500\,\mathrm{K} , which are slightly lower than theoretical predictions of ∼ 5000 ​ K \sim 5000\,\mathrm{K} from the black hole envelope models [ 45 , 40 ] but consistent with the JWST/PRISM spectroscopy demographic study of de Graaff et al. [20] .

S5: V Summary / block 1 · html

In this work, we present the detection of a very close pair of LRDs found in a strongly lensed galaxy and separated by ∼ 70 ​ pc \sim 70\,\mathrm{pc} , which are embedded in the typical SFGs at z ∼ 7 z\sim 7 . The main findings in this work are summarized below:

S7: VII Gravitational-wave background / block 1 · html

The resulting SGWBs are shown in Fig. 3 (thick lines). The WD capture background peaks at deciHz and are observable with an SNR > 10 \rm SNR>10 when compared with the corresponding LGWA PLS [ 1 ] and LISA PLS [ 47 ] for T obs = 4 ​ yr T_{\rm obs}=4\,\rm yr . For comparison, we also show the background from massive black-hole mergers under two seeding scenarios [ 9 ] . While the BH and NS backgrounds have SNR < 10 \rm SNR<10 within four years, there are still thousands of captures within z < 6 z<6 with LISA. Finally, the SGWB from inspirals is undetectable.

S9: IX Conclusions / block 1 · html

We modeled stellar feeding onto SMBHs in NSCs, including TDEs and compact-object captures. Stellar plunges contribute only a sizable fraction to SMBH mass growth; however, gas accretion still dominates at all redshifts.

S2.SS2.SSS3: II.2.3 Critical collisional parameters to retain the BH / block 1 · html

Whether the sBH can be captured also depends critically on the impact parameter b b . Taking a representative sBH with M ∙ = 10 ​ M ⊙ M_{\bullet}=10\,M_{\odot} , we conduct a grid search over b b and v ∞ v_{\infty} to map the retention threshold. The results are displayed in the right panel of Figure 2 , where the red dashed line represents the retention boundary, which is well-fitted with v ∞ / v 0 = − 0.352 ​ ( b / R ⋆ ) 2 − 0.572 ​ ( b / R ⋆ ) + 1.35 v_{\infty}/v_{0}=-0.352(b/R_{\star})^{2}-0.572(b/R_{\star})+1.35 , and the solid red line marks the geometric collision boundary. At low v ∞ v_{\infty} , the critical impact parameter for retention is comparable to the geometric collision threshold. However, as v ∞ v_{\infty} increases, the critical retention impact parameter deviates rapidly. Ultimately, when v ∞ ∼ 1.4 ​ v 0 v_{\infty}\sim 1.4\,v_{0} , only head-on ( b = 0 b=0 ) collisions can marginally retained the sBH inside the star. However, for most sBHs retained in the star, the stopping time is still ∼ 10 \sim 10 – 100 ​ h 100\,\rm h .

S3: III Hydrodynamical Simulations / block 1 · html

The results in Section 3 indicate that energy dissipation during the sBH’s orbital decay can lead to changes in the massive main-sequence star’s interior. However, this 1-D MESA model does not fully capture the rapid 3-D dynamical response of the stellar envelope. In this section, we perform a suite of hydrodynamical simulations of star–sBH collisions, using the Lagrangian code GIZMO [ 21 ] . The simulations solve the standard hydrodynamical equations 1 1 1 Specifically: (1) ∂ t ρ + ∇ ⋅ ( ρ ​ 𝒗 ) = 0 \partial_{t}\rho+\nabla\cdot(\rho\bm{v})=0 ; (2) ∂ t ( ρ ​ 𝒗 ) + ∇ ⋅ ( ρ ​ 𝒗 ​ 𝒗 + p ​ 𝐈 ) = ρ ​ 𝒈 \partial_{t}(\rho\bm{v})+\nabla\cdot(\rho\bm{v}\bm{v}+p\mathbf{I})=\rho\bm{g} ; (3) ∂ t ( ρ ​ e + ρ ​ v 2 / 2 ) + ∇ ⋅ [ ( ρ ​ e + ρ ​ v 2 / 2 + p ) ​ 𝒗 ] = ρ ​ 𝒈 ⋅ 𝒗 \partial_{t}\left(\rho e+\rho v^{2}/2\right)+\nabla\cdot\left[\left(\rho e+\rho v^{2}/2+p\right)\bm{v}\right]=\rho\bm{g}\cdot\bm{v} ; (4) p = ( γ − 1 ) ​ ρ ​ e p=(\gamma-1)\rho e . with the code’s Meshless Finite-Mass (MFM) mode, and self-gravity is solved with tree-gravity [ 20 ] .

S4.SS3: IV.3 M BH − M ∗ M_{\rm BH}-M_{*} relation / block 1 · html

One source has a modest M BH / M ∗ ∼ 0.1 % M_{\rm BH}/M_{*}\sim 0.1\% . This value is considerably lower and aligns more closely with the relation observed in the local universe [ 65 , ∼ 0.5 % \sim 0.5\% ; ] . This result corroborates recent discoveries of low-luminosity AGNs in massive galaxies that also follow the local trend [ 69 ] . Four sources (ID 11026, ID 46938, ID 12248, and ID 39803) have a significantly higher BH-to-stellar mass ratio of M BH / M ∗ ≥ 2 % M_{\rm BH}/M_{*}\geq 2\% , which are consistent with faint AGNs identified in previous studies [ 47 , 80 , e.g., ] .

A2: Appendix B All fitting results / block 1 · html

Figure B.2: Fitting results. Similar to that shown in Figure B.1 .

S2.SS2: 2.2 Stars and black holes / block 1 · html

Stars with masses in the type II (11 - 20 M ⊙ ) supernovae (SNe) or the hypernova (20 - 40 M ⊙ ) mass bands will experience an explosion and undergo some degree of mass loss. The explosion energy and resulting black hole remnant mass are calculated from the stellar mass as taken from Nomoto et al. (2006) , and the explosion is modelled with an isotropic injection of thermal energy into the surrounding 10 pc. The particle is then converted into a black hole particle.

S5: 5 Conclusions / block 1 · html

Here we have presented a new cosmological simulation using the SEEDZ framework, with enhanced resolution compared to the original suite. The simulation was performed with the hydrodynamics code Arepo2 and probes scales of ∼ \sim 0.5 physical pc. On top of the base code, the SEEDZ framework includes subgrid models for PopIII star formation, supernovae explosions, metal enrichment, PopII star formation, black hole formation, and heavy seed black hole formation. Heavy seeds form in regions where the mass inflow into 10 pc regions surrounding a gravitationally unstable sub-pc sized cell exceeds 1 M ⊙ yr -1 . By splitting our catalogue of early galaxies into two groups – those that form heavy seeds and those that do not – we investigate the conditions necessary for high inflow environments to occur and by extension the conditions necessary for SMS and/or dense stellar cluster formation. We make the following conclusions: