EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 711 (July 2026) A&A, 711 (2026) A57 Full HTML
Open Access
Issue
A&A
Volume 711, July 2026
Article Number A57
Number of page(s) 13
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202659533
Published online 02 July 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1 Introduction

Young massive clusters (YMCs) are dense stellar systems with typical masses Mcl ≳ 104 M⊙, half-mass radii of a few parsecs, and ages ≲100 Myr (Portegies Zwart et al. 2010; Bastian et al. 2013). They are observed in a wide range of environments, from nearby star-forming galaxies to extreme starbursts, and may represent the progenitors of globular clusters (de Grijs 2009; Portegies Zwart et al. 2010; Bastian et al. 2013; Krumholz et al. 2019). Owing to their high stellar densities and large populations of massive stars, YMCs provide key laboratories for studying stellar dynamics, binary evolution, and the formation of compact objects, including stellar-mass black holes (BHs), gravitational-wave (GW) sources, and possible seeds of intermediate-mass black holes (IMBHs) with mass mBH ∼ 102–105 M⊙ (Portegies Zwart & McMillan 2002; Portegies Zwart et al. 2004; Mapelli 2016; Di Carlo et al. 2019, 2020; Kremer et al. 2020a; Di Carlo et al. 2021; Rastello et al. 2020, 2021; Arca-Sedda et al. 2021; Askar et al. 2024; Rantala et al. 2024, 2025; Vergara et al. 2025).

In sufficiently dense clusters, stellar collisions are expected to occur frequently, particularly among massive stars that rapidly segregate toward the cluster center (e.g., Hills & Day 1976; Spitzer 1988; Binney & Tremaine 2008). Early numerical studies showed that such collisions may trigger runaway growth, leading to the formation of very massive stars (VMSs) and, potentially, IMBHs. This process is especially relevant in clusters with short core-collapse times, where collisions are expected to begin before massive stars evolve into compact objects (Portegies Zwart & McMillan 2002; Portegies Zwart et al. 2004; Gaburov et al. 2008). Repeated stellar collisions have been shown to be effective in clusters with masses >105 M⊙ and central densities >5×106 M⊙ pc−3 (Rizzuto et al. 2021; González et al. 2021; Arca Sedda et al. 2023; González Prieto et al. 2024; Rantala et al. 2024, 2025, 2026), corresponding to structural parameters that are rarely observed in local YMCs. Whether stellar collisions can efficiently promote the growth of VMSs in YMCs more commonly found in the local Universe therefore remains uncertain, partly because of the computational challenge of simulating massive clusters with high primordial binary fractions (i.e., the fraction of binaries born with the cluster).

The dynamical evolution of YMCs has been investigated using a variety of numerical techniques, including Monte Carlo methods (e.g., Giersz et al. 2015; Hong et al. 2020; Rodriguez et al. 2022; Sharma & Rodriguez 2025), direct N-body simulations (e.g., Wang et al. 2016; Di Carlo et al. 2020; Rastello et al. 2021; Arca Sedda et al. 2024a; Banerjee 2025, 2026), and hybrid N-body approaches (e.g., Wang et al. 2020a; Barber & Antonini 2025), i.e., N-body codes based on Hamiltonian splitting of short- and long-range interactions. The most recent N-body simulations have mostly focused on exceptionally high densities at half-mass radius (>106 M⊙ pc−3, e.g., Rantala et al. 2025), attainable only under very extreme conditions such as those in the early Universe (e.g., Vanzella et al. 2023; Adamo et al. 2024). In these environments, repeated stellar collisions efficiently produce IMBHs with masses >104 M⊙.

In this work, we present the TITANS simulation set, a new suite of N-body simulations performed with the state-of-the-art N-body code PETAR1. Our models span wide ranges of cluster masses (8 × 104–9 × 105 M⊙) and half-mass densities (100–105 M⊙ pc−3), compatible with current globular clusters and YMCs (Portegies Zwart et al. 2010), while adopting observation-based primordial binary fractions (Moe & Di Stefano 2017). While the ultimate goal is to compare the properties of these clusters to those of local star clusters, we focus here on the first 20 Myr of cluster evolution (Teen TITANS), corresponding to the phase in which massive stars are still evolving into compact objects or pair-instability supernovae (PISNe). Using realistic initial conditions for YMCs (Portegies Zwart et al. 2010; Bastian et al. 2013; Krumholz et al. 2019) and stellar evolution prescriptions that include PISNe, we self-consistently study the efficiency of repeated stellar collisions, the formation and evolution of VMSs, and their impact on the BH mass spectrum in YMCs.

The paper is structured as follows. Section 2 presents the code used and the chosen initial conditions. Section 3 reports the main results on repeated stellar collisions and their final stellar products. Section 4 discusses the impact of our results on the formation of PISNe and on the BH mass distributions in our clusters. Finally, Section 5 summarizes our findings.

2 Methods

The TITANS set comprises 18 N-body simulations of star clusters in a Milky Way-like galaxy, taking into account single and binary stellar evolution, and considering a set of structural parameters compatible with observational limits in local YMCs. We performed this set of simulations with the hybrid N-body code PETAR (Wang et al. 2020a). In particular, PETAR employs the particle-tree particle-particle algorithm (P3T; Oshino et al. 2011), and the slow-down algorithmic regularization (SDAR; Wang et al. 2020b) scheme, which combines a fourth-order Hermite integrator with the SDAR method to accurately treat close encounters. SDAR reduces the number of required timesteps when integrating weakly perturbed few-body systems by effectively decoupling the perturbation timescale from the intrinsic dynamical timescale, enabling an efficient and accurate treatment of binaries. Thanks to this approach, PETAR can simulate massive clusters with large binary fractions while maintaining short integration times and small relative energy errors compared to direct N-body codes (e.g., NBODY6++GPU, Wang et al. 2015). The use of OpenMP for an efficient workload distribution across central processing units (CPUs), combined with the graphic processing unit (GPU) acceleration, contributes significantly to the strong parallel performance and scalability of PETAR. This code can include single and binary star evolution through population-synthesis algorithms. Here, we use the MOBSE population-synthesis code (Mapelli et al. 2017; Giacobbo & Mapelli 2018). PETAR can also model the effect of an external potential with GALPY (Bovy 2015).

We performed our simulations on two servers. The first, demoblack, hosted by the University of Padova, with 8 Nvidia Tesla V100 GPUs and 8 INTEL Xeon Platinum CPUs with 24 cores each. The second, bwForCluster Helix, on a partition hosting 4–8 Nvidia A100 GPUs and 2 AMD EPYC CPUs with 64 cores each.

2.1 Stellar and binary evolution

We modeled single and binary stellar evolution using the population synthesis code MOBSE (Mapelli et al. 2017; Giacobbo & Mapelli 2018, 2020), which is an upgraded custom version of BSE (Hurley et al. 2002). MOBSE includes an updated treatment of radiation-driven stellar winds (Chen et al. 2015), natal kicks (Giacobbo & Mapelli 2020), (pulsational) PISNe (Mapelli et al. 2020), and compact-object formation (Fryer et al. 2012). Specifically, the radiation-driven wind models in MOBSE account for Thomson scattering in radiation-pressure dominated stars. Following Chen et al. (2015), the mass loss scales as ˙M ∝ Zβ, where Z is the metallicity and β = 0.85 if the electron-scattering Eddington ratio of a star is Γ < 2/3, β = 2.45–2.4 Γ if 1 > Γ ≥ 2/3, and β = 0.05 if Γ > 1. This model self-consistently increases the wind mass loss in massive stars (≳100 M⊙).

Here, we assumed a metallicity of Z = 0.0002 for all clusters. Metal-poor YMCs are expected to contribute more significantly to GW sources than metal-rich YMCs (such as the ones in the Milky Way), as they are able to produce massive compact remnants (Mapelli 2016; Rantala et al. 2024; Arca Sedda et al. 2024a; Barber & Antonini 2025).

Binary evolution processes, such as tides, mass transfer, common-envelope evolution, and GW-driven orbital decay, were implemented following the prescriptions of Hurley et al. (2002). In this work, we modeled the common-envelope phase using the standard αCE–λCE formalism, with αCE the fraction of orbital energy used to unbind the envelope, and λCE a parameter characterizing the envelope’s binding energy and structure. In particular, we adopted αCE = 1, while λCE is calculated according to Claeys et al. (2014).

For core-collapse SNe, we relied on the delayed formalism by Fryer et al. (2012). This model adopts a convection-enhanced neutrino-driven mechanism for the SN explosion, and the revival of the shock wave happens ∼250 ms after the collapse. With this prescription the minimum BH mass is 3 M⊙.

After the SN explosion, compact remnants receive a natal kick modeled following Giacobbo & Mapelli (2020): υkick=fH05〈 mNS 〉mremmej〈 mej 〉,Mathematical equation: ${\upsilon _{{\rm{kick}}}} = {f_{{\rm{H}}05}}{{\left\langle {{m_{{\rm{NS}}}}} \right\rangle } \over {{m_{{\rm{rem}}}}}}{{{m_{{\rm{ej}}}}} \over {\left\langle {{m_{{\rm{ej}}}}} \right\rangle }},$(1)

where ⟨mNS⟩ and ⟨mej⟩ are the average neutron-star mass and ejecta mass, computed from a population of isolated neutron stars with metallicities representative of the Milky Way, mrem is the mass of the compact object, and mej is the ejected mass. The term fH05 is a number randomly drawn from a Maxwellian distribution with one-dimensional root mean square σkick = 265 km s−1, derived from the proper motions of young Galactic pulsars (Hobbs et al. 20052). This formalism implies that BHs forming through direct collapse receive a null natal kick. In this study we have not included a formalism for the natal spin of BHs.

The fitting formulas for PISNe and pulsational PISNe were taken from the appendix of Mapelli et al. (2020), based on Woosley (2017). When the helium core mass of a star, mHe,f, is between 32 and 64 M⊙, it undergoes a pulsational PISN. As a consequence of the adopted wind model and formalism for pulsational PISNe, MOBSE allows for the formation of BHs as massive as 65 M⊙ at Z = 0.0002 (Giacobbo & Mapelli 2018). A PISN is triggered instead when 64 ≤ mHe,f ≤ 135 M⊙. For isolated single stellar evolution at the metallicity adopted in our simulations, this corresponds to a range in the zero-age main-sequence mass 140 ≲ mZAMS ≲ 260 M⊙. When mHe,f > 135 M⊙, the star directly collapses into a BH.

2.2 Compact object mergers in PETAR

In dense star clusters with a high fraction of primordial binaries, collisions and mergers between stars and BHs can happen frequently (see e.g., Rastello et al. 2026). The outcome of such interactions is still unclear. On the one hand, we expect that stars with m∗ < 10 M⊙ would accrete little or no mass on the BH (Kremer et al. 2020b). On the other, simulations from Schrøder et al. (2020) show that common envelope events between a star and a BH can lead to the accretion of the stellar core on the BH and to the ejection of the envelope. Moreover, various dynamical studies (Kıroğlu et al. 2025a,b) have found that BHs can accrete mass from repeated collisions with stars in a dense cluster. As a consequence, previous works (Banerjee 2021; Rizzuto et al. 2021; Arca Sedda et al. 2024a) have assumed a fraction of stellar accreted mass on the BH fc > 0. In this work, we adopted a more conservative approach, assuming fc = 0.

In PETAR, GW–driven mergers are modeled by evolving the semi-major axis and eccentricity of the binary according to the formalism of Peters (1964). A binary is considered to merge when its GW timescale becomes shorter than its integration timestep, which depends on the strength of nearby perturbations and the binary’s slowdown factor. Once this condition is met, PETAR evolves the binary’s position to the merger time, allowing us to locate binary BH (BBH) mergers within the cluster.

2.3 External potential

We modeled the influence of an external galactic potential using the GALPY package (Bovy 2015). Specifically, we assumed that the star clusters resided in a Milky Way–like galaxy and followed a circular orbit around the Galactic center at a radius of 8 kpc with a velocity of 220 km s−1, corresponding to the Solar neighborhood. While this choice is not fully representative of the diverse orbital properties of Milky Way YMCs (Portegies Zwart et al. 2010; Bastian et al. 2013), it provides a useful approximation of the tidal effects experienced by clusters in a typical Galactic environment. For the galactic potential, we adopted the MWPotential2014 model, which includes a Miyamoto–Nagai potential for the disk (Miyamoto & Nagai 1975), a spherical bulge component based on a truncated power-law density profile (Hernquist 1990), and a Navarro–Frenk–White potential for the dark matter halo (Navarro et al. 1996).

2.4 Initial conditions

The initial conditions of our simulations were generated using MCLUSTER (Küpper et al. 2011). We adopted initial cluster masses in the range Mcl ∈ [8 × 104, 9 × 105] M⊙, with positions and velocities sampled from a King density profile (King 1966) characterized by a central potential parameter W0 = 6 and half-mass radii rh ∈ [1, 5] pc. These values are consistent with those measured for YMCs in the local Universe (Portegies Zwart et al. 2010; Bastian et al. 2013; Krumholz et al. 2019). The initial density at half-mass radius ρh of our clusters is between 100 and 105 M⊙ pc−3 and the half-mass relaxation time trh (Spitzer 1988) is between 100 Myr and 3 Gyr. We assumed that our clusters were not initially segregated. A summary of the initial conditions for our suite of simulations is provided in Table 1.

We sampled stellar masses from a Kroupa initial mass function (Kroupa 2001) between 0.08 and 150 M⊙. The global binary fraction fb is defined as fb = Nb/(Ns + Nb) and its value is ∼0.23 across all the models, corresponding to ∼38% of stars in a primordial binary. We assumed that the primordial binary fraction depends on the primary mass of the star, according to the fitting formulas by Moe & Di Stefano (2017), which are consistent with observations of stars across different environments, including massive stars in young clusters. With this assumption, the resulting binary fraction is fb ≈ 1 among massive stars and decreases to ≲0.4 among solar and subsolar mass stars.

We initialized binaries with distributions of mass ratio, semi-major axis, and eccentricity motivated by the intrinsic properties of massive binaries inferred by Sana et al. (2012). We extended the orbital period range up to 105 days to include wide, non-interacting binaries. In addition, we extrapolated the mass-ratio distribution to q < 0.1 to allow for extreme mass-ratio systems.

Figure 1 shows the location of our simulations in the N–ρh and N– fb planes, where N is the number of particles. We compare our models with state-of-the-art studies performed using Monte Carlo methods (e.g., Askar et al. 2017; Rodriguez et al. 2019; Maliszewski et al. 2022), direct N-body simulations (e.g., Mapelli et al. 2013; Mapelli 2016; Wang et al. 2016; Arca Sedda et al. 2024a; Rantala et al. 2024), and hybrid N body approaches (e.g., Wang et al. 2021; Rastello et al. 2026). The main novelty of the TITANS simulation set lies in the fact that, although it explores a region of the N − ρh parameter space already covered by previous studies, it does so adopting a global primordial binary fraction fb that has so far been investigated almost exclusively using Monte Carlo techniques (Askar et al. 2017; Maliszewski et al. 2022). Importantly, our simulations sample a region of parameter space that is representative of YMCs and globular clusters in the local Universe, whose typical densities are ≲106 M⊙ pc−3.

The total simulated time for each cluster, tsim, exceeds 20 Myr in all cases. For the first 20 Myr, we recorded the output snapshots every 0.125 Myr to track the evolution of each star cluster in greater detail. We choose to focus here on the first 20 Myr of cluster evolution, corresponding to the lifetime of a ∼12 M⊙ star; hence, we expect that all BHs have formed by this time and most SNe have already taken place. Moreover, within 20 Myr, the densest clusters contract and undergo core collapse (Spitzer 1988). After tsim > 20 Myr, output snapshots are recorded every 1 Myr.

Table 1

Initial conditions of the TITANS.

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Initial binary fraction (top row) and density at half-mass radius (bottom row) as a function of the initial number of stars, for Monte Carlo (green squares), direct N-body (violet circles), and hybrid N-body simulations (brown stars). The TITANS simulations are represented by a red star.

2.5 Properties of stellar collisions

Stellar collisions and repeated stellar collisions have been extensively investigated as formation channels of IMBH progenitors (Portegies Zwart & McMillan 2002; Portegies Zwart et al. 2004; Giersz et al. 2015; Mapelli 2016; Di Carlo et al. 2021; González et al. 2021; Rizzuto et al. 2021; Kritos et al. 2023; Reinoso et al. 2023; Arca sedda et al. 2024b; Rantala et al. 2024, 2025; Vergara et al. 2025; Paiella et al. 2026). In this work, we explore the properties of repeated stellar collisions in YMCs with a high fb, and features typical of local star clusters.

Here, we define a stellar collision chain as a sequence containing at least two stellar collisions or mergers. Throughout the paper, we use the term “stellar collision” in a broad sense, referring both to mergers of binary components and to direct collisions between unbound stars. We quantify the efficiency of repeated stellar collisions as ηrep = Nch,rep/Mcl, where Nch,rep is the number of stellar collision chains identified in each cluster. Figure 2 represents schematically a stellar collision chain, initiated either by the merger of a primordial or dynamical binary (i.e., a binary that is not born with the cluster), or by a stellar collision. The resulting star will itself either form a dynamical binary that will later merge, or collide with another star. We expect this process to reiterate itself for several times, especially in dense clusters with small relaxation times. When the remnants from two or more chains merge, we compute Nch,rep as the sum of the length of the chains.

We also examine how repeated stellar collisions affect the mass spectrum of VMSs and IMBHs. All through this study, we define VMSs as stars whose mass exceeds the upper limit of the adopted initial mass function (m∗ > 150 M⊙), and IMBHs as BHs with mBH > 100 M⊙.

Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

Schematic representation of a stellar collision chain. We show on the left a chain initiated by the merger of the components of a binary (either primordial or dynamically formed). On the right we show instead a chain initiated by an hyperbolic stellar collision.
Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

Efficiency of repeated stellar collisions ηrep as a function of the half-mass density ρh, normalized by ρh,3 and Mcl,5. The points are colored as a function of the cluster initial relaxation time trh.

3 Results

In the following Sections, we discuss the properties of repeated stellar collisions (Sect. 3.1), and their impact on the formation of VMSs (Sect. 3.2). We show our results considering only the first 20 Myr of evolution of each cluster.

3.1 Repeated stellar collisions

3.1.1 Impact of the host environment

Figure 3 shows the efficiency of repeated stellar collisions as a function of ρh/(ρh,3 Mcl,5), with ρh,3 = 103 M⊙ pc−3 and Mcl,5 = 105 M⊙. As defined in Sect. 2.5, the efficiency of repeated collisions is ηrep = Nch,rep/Mcl, where Nch,rep is the number of stellar collision chains identified in each cluster. As cluster density increases, so does the efficiency of repeated collisions. We fit the data with a linear relation in logarithmic space, obtaining log10(ηrep)≈0.8log10(ρhρh,3Mcl,5)−1.4 with ρh≳500M⊙ pc−3.Mathematical equation: ${\log _{10}}({\eta _{{\rm{rep}}}}) \approx 0.8{\log _{10}}\left( {{{{\rho _{\rm{h}}}} \over {{\rho _{{\rm{h}},3}}{M_{{\rm{cl}},5}}}}} \right) - 1.4{\rm{with}}{\rho _{\rm{h}}} \mathbin{\lower.3ex\hbox{$\buildrel>\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 500{\rm{}}{M_ \odot }{\rm{p}}{{\rm{c}}^{ - 3}}.$(2)

The intrinsic scatter around the best-fit relation is σint ∼ 0.3 dex. This trend indicates that a high cluster density is a key requirement for the onset of multiple stellar collisions. Moreover, we find that repeated collisions occur only for ρh ≳ 500 M⊙ pc−3.

An opposite trend is observed with respect to the relaxation time trh. At fixed density, clusters with shorter relaxation times exhibit a higher efficiency of repeated collisions. In fact, an efficient formation of merger chains requires a mass-segregation time tDF ∼ ⟨m⟩ trh/mbin shorter than the lifetime of the stars involved, tDF < tend, where the average stellar mass is ⟨m⟩ ∼ 0.6 M⊙ and the binary mass satisfies mbin ≤ 300 M⊙. For the most massive binaries in our clusters (mbin = 300 M⊙, m1 = m2 = 150 M⊙), the stellar lifetime is tend ∼ 3 Myr, corresponding to a maximum relaxation time trh ∼ 1.5 Gyr. Less massive binaries will require longer times to segregate to the cluster center.

As shown in Table 1, more than half of the TITANS models have trh ≲ 1.5 Gyr. As a consequence, despite their relatively long core-collapse times (tcc > 20 Myr), these clusters can form stellar collision chains because the most massive primordial binaries segregate rapidly to the center. Once in the core, these binaries merge, and the increasing central density of massive stars favors the formation of new dynamical binaries or direct stellar collisions. This result is consistent with the findings of Rantala et al. (2025) on the key role of primordial binaries in clusters with long core-collapse times. In contrast, clusters with trh ≳ 1.5 Gyr require more time for binaries to segregate and interact; by the time they reach the core, the most massive stars are often already evolved or have died as compact objects, reducing the efficiency of chain formation.

Figure 4 shows the number of stellar collision chains Nch,rep as a function of the normalized half-mass density ρh/ρh,3 and of the length of the chains themselves. The minimum density required for a chain to form increases with the number of collisions it contains: chains with two collisions appear at ρh ∼ 500 M⊙ pc−3, chains with three collisions require densities of ∼4000 M⊙ pc−3, and chains with four or more collisions form only at densities above ∼104 M⊙ pc−3. In all cases, Nch,rep increases approximately linearly with density. However, the slope decreases for longer chains, indicating that extended collision sequences are rare even in clusters with high densities and low relaxation times. As a result, the contribution of long collision chains to Nch,rep and ηrep remains limited. In the densest model of our suite (T16, ρh ∼ 105 M⊙ pc−3), 70% of the chains contain two collisions, 20% contain three collisions, and only 10% involve more than four collisions. We point out that models T2, T6, and T7, do not produce any stellar collision chains. Models T3, T10, T13, and T18 instead produce a single chain each (Nch,rep = 1), always containing two collisions.

Finally, we compute the fraction of stellar collision chains whose first merger involves the components of a primordial binary. This fraction is 100% up to ρh ∼ 104 M⊙ pc−3 and remains above ∼72% at higher densities. When the first collision does not involve a primordial binary merger, it results either from a hyperbolic encounter or from a highly eccentric collision between a low-mass main-sequence (MS) star and one component of a primordial binary. In the latter case, the perturber can be either a single star or a component of a low-mass primordial binary, which is disrupted during the interaction. Owing to the small mass ratio between the perturber and the massive primordial binary component, this first merger does not significantly alter the orbit of the primordial binary, whose components subsequently merge.

Overall, these results highlight that stellar collision chains are efficient only in clusters with ρh ≳ 104 M⊙ pc−3 and with trh ≲ 300 Myr, and that, even then, they mostly contain two collisions and are triggered by a primordial binary merger. In Sects. 3.2 and 4.3 we explore the effect of these findings on VMSs and IMBHs. Finally, we emphasize that the stellar collision chains identified in our simulations are not in the strong runaway regime, where a single growing object dominates the cluster core and alters the local potential. Instead, they involve only a small number of mergers and do not lead to the formation of a central massive object. As a result, the cluster potential and dynamics remain largely unaffected by individual collision products.

Thumbnail: Fig. 4 Refer to the following caption and surrounding text. Fig. 4

Number of stellar collision chains Nch,rep as a function of the normalized density at half-mass radius ρh/ρh,3. In violet we represent the number of chains containing two collisions, in blue the number of chains containing three collisions, and in green the number of chains containing four or more collisions. The dashed lines and annotations represent the linear fits of the points.
Thumbnail: Fig. 5 Refer to the following caption and surrounding text. Fig. 5

Primary stellar type of stars involved in repeated collisions. The x-axis shows the simulated models. The stellar phases on the y-axis are: MS, shell hydrogen burning (sHB), core helium burning (cHeB), shell helium burning (sHeB), and Wolf– Rayet star (WR). The color represents the relaxation time at half-mass radius trh.

3.1.2 Stellar properties

Figure 5 shows the stellar types of primary stars involved in repeated collisions in the simulated models, colored by the relaxation time of their host cluster. We do not show the stellar types of secondaries because in all models ≳93% stars are in their MS. Moreover, we do not show models T2, T6, and T7 since they do not produce stellar collision chains. In about half of the simulated models, the primary stars involved in repeated stellar collisions are preferentially in their MS. With increasing trh the fraction of evolved stars involved in repeated stellar collisions grows. In fact, longer relaxation times correspond to longer mass segregation timescales and, as a consequence, stars are more evolved by the time they reach the dense cluster center and take part in the collisions. Moreover, low densities lead to bigger intervals of time between consecutive collisions, allowing the stars to evolve. We move from model T1 (trh ∼ 100 Myr) with ∼77% primaries in the MS, to model T17 (trh ∼ 3 Gyr) with only ∼17% primaries in their MS. Models T3, T10, and T18, appearing as outliers, produce only one chain of stellar collisions. We emphasize that this result depends on stochastic fluctuations, as we ran a single realization for each simulation.

Figure 6 shows the primary masses of all the stars involved in repeated stellar collisions. The bars contain all the mass values, while the rectangles contain 50% of the data. We see that the bulk values of primary masses do not strongly depend on the properties of the host clusters. Dense clusters with trh ≤ 700 Myr (e.g., T1, T9, T14, T16) have larger bulk values (50 ≲ m1 ≲ 130 M⊙) and they cover broader mass intervals. In fact, repeated stellar collisions are more frequent in dense clusters with low relaxation times (Fig. 3), and they mainly involve stars in their MS (Fig. 5). Model T9 spans a primary mass interval extending up to 350 M⊙, due to the long chains of repeated stellar collisions forming in this cluster. Models T13 and T18, appearing as outliers, have Nch,rep = 1. We also represent here the maximum mass reached in our clusters through stellar collision chains; we see that the highest values are reached in models T1, T9, and T16, which have the lowest relaxation times. For secondary masses, we find that the bulk values are usually below m2 ≲ 25 M⊙ if trh ≳ 500 Myr.

3.2 Very massive stars

We summarize here the properties of the VMSs (m∗ > 150 M⊙) forming in our simulated clusters. Table 2 lists, for each model, the number of VMSs, their maximum mass, and their formation channels.

The number of VMSs generally increases with cluster mass, consistently with the fact that their dominant formation pathway is the merger of primordial binaries, whose abundance scales with Mcl. In the densest clusters with short relaxation times, stellar collision chains play a non-negligible role in VMS formation. This is especially evident in models T1, T9, and T16, where more than 25% of VMSs originate from repeated collisions. Models T9 and T16 form a larger number of VMSs than clusters of comparable mass because repeated stellar collisions are more efficient in these systems. This can be seen from the normalized number of VMSs as well, which is higher for these two models. In all the other models, the VMSs counted in NVMS,coll form predominantly (>82%) via primordial binary mergers. This result highlights again the key role of primordial binaries in driving the early dynamical and stellar evolution of YMCs with properties in line with low-redshift observations (Portegies Zwart et al. 2010; Rantala et al. 2025). We find that binary evolution also contributes to VMS formation: in particular, stable mass transfer can produce stars with masses up to ∼200 M⊙. Across all clusters, the median VMS mass shows only a weak dependence on cluster properties and remains in the range 160–200 M⊙. All VMSs identified in our simulations are still on the MS, as their formation typically occurs through stellar collisions involving young stars.

Table 2 also reports the maximum stellar mass attained in each simulation. Figure 7 shows it as a function of the initial half-mass density of the clusters. Since VMSs form early during cluster evolution (t < 5 Myr), this trend remains unchanged when considering the density at the time of their formation. For clusters with ρh ≲ 104 M⊙ pc−3 and trh > 300 Myr, the maximum stellar mass is nearly independent of density and relaxation time, lying between 215 and 300 M⊙, with stochastic variations reflecting the collisional nature of VMS formation. In all these cases, VMSs originate from primordial binary mergers. In models with ρh ≳ 104 M⊙ pc−3 and trh ≲ 300 Myr, the maximum stellar mass exceeds 330 M⊙. The three stars in this regime form in models T1, T9, and T16 through stellar collision chains.

Finally, the last column of Table 2 reports the number of VMSs that merge with a stellar-mass BH. Appendix A discusses the final fate of these stars and the impact of our initial prescriptions. The subsequent evolution of VMSs is instead addressed in Sect. 4.

Thumbnail: Fig. 6 Refer to the following caption and surrounding text. Fig. 6

Distribution of primary mass m1 of stars involved in repeated stellar collisions for each simulated model. The color represents the relaxation time at half-mass radius trh. We show with a red star the maximum mass of the product of a stellar collision chain in each model.
Table 2

Properties of VMSs in the simulated models.

4 Discussion

In the following, we discuss the impact of stellar collisions and repeated stellar collisions on PISNe (Sect. 4.1) and on the mass distribution of BHs (Sect. 4.2). We focus particularly on the formation and properties of IMBHs in Sect. 4.3. Also here, we consider only the first 20 Myr of evolution of each cluster, unless otherwise specified. In Sect. 4.5 we report the caveats and limitations of our models.

Thumbnail: Fig. 7 Refer to the following caption and surrounding text. Fig. 7

Relation between the maximum VMS mass mVMS,max and the initial density at half-mass radius ρh. The colors used represent the initial relaxation time trh.

4.1 Pair-instability supernovae

We expect that stellar collisions increase the number of stars entering the PISN mass range relative to predictions from the initial mass function, potentially enhancing the rates of PISNe and pulsational PISNe in metal-poor star clusters compared to the field. In our simulations, most PISNe originate either from stars that are born directly within the PISN mass range or from stars that enter it following a primordial binary merger. Since none of the PISN progenitors escape from the clusters, both formation channels scale with the cluster mass Mcl. This behavior is confirmed in Fig. 8, which shows that both the total number of PISNe and the fraction produced through stellar collisions correlate strongly with Mcl. The contribution of stellar collisions to the total PISN population is above ∼70%, with large scatter among clusters with similar properties, reflecting the stochastic nature of stellar collisions and massive binary evolution. Only four PISNe (3 in model T9, 1 in model T16) form via repeated stellar collisions. We find that ≳90% of our PISNe explode as single stars, while the remaining part is found in a binary system. The only exception is model T1, where 50% of the PISNe explode in a binary. This difference is driven by small-number statistics, as T1 is the least massive cluster in our sample and forms only two PISNe.

Figure 8 also shows the expected number of PISNe from an isolated population of single stars, assuming that the PISN mass gap starts at ∼140 M⊙ (Spera & Mapelli 2017; Giacobbo et al. 2018). We find that the number of stars born in our clusters within the PISN mass range is consistent with the isolated case. However, stellar collisions (primarily in the form of primordial binary mergers) and, to a lesser extent, repeated stellar collisions, allow additional stars to enter the PISN regime, increasing the total number of PISNe. As a result, in most of our models the total number of PISNe forming in star clusters is up to one order of magnitude larger than expected for an isolated single stellar population.

In our simulations, the PISN mass gap corresponds to 140 ≲ m∗ ≲ 400 M⊙ considering the final stellar mass before the explosion, with the bulk of the population <289 M⊙ (e.g., Arca Sedda et al. 2024a). We expect these results to be robust against uncertainties in stellar-wind prescriptions owing to the low metallicity adopted in this work (Simonato et al. 2025).

Finally, we find that ≳70% of the VMSs formed in the clusters explode as PISNe. The remaining VMSs typically lose a substantial fraction of their mass and collapse into stellar-mass black holes with mBH < 100 M⊙ as a result of pulsational PISNe (see Sect. 4.2), or are removed through collisions with stellar-mass black holes (see Appendix A). The only exception is model T1, in which ∼33% of the VMSs die as PISNe.

Thumbnail: Fig. 8 Refer to the following caption and surrounding text. Fig. 8

Number of PISNe as a function of the cluster mass Mcl. The filled red points represent the total number; the empty black points represent instead the number of PISNe due to collisions. As a comparison we plot the number of PISNe expected from single stellar evolution NPISN,single.

4.2 Black hole mass distribution

To assess the impact of single and repeated stellar collisions on the BH mass spectrum, we focus on three models with different values of density and relaxation time (T2, T13, T16). Figure 9 shows the corresponding mass distributions of all the BHs forming in the simulations within 20 Myr, highlighting the contributions from stellar collisions and stellar collision chains. As expected, the total number of BHs increases with the mass of the host cluster, ranging from 721 to 2892 BHs in the most massive model (T16). In parallel, the number of BHs originating from single stellar collisions also grows with cluster mass, accounting for ∼28–35% of the total BH population (Di Carlo et al. 2019; Arca Sedda et al. 2024a), reflecting the increasing number of massive stars and primordial binaries in more massive systems. A qualitatively different behavior is observed only in model T16, the densest cluster in our suite (ρh ∼ 105 M⊙ pc−3) with a short relaxation time (trh = 300 Myr). In this case, a non-negligible fraction (3%) of BHs forms from the collapse of the final products of stellar collision chains.

The BH mass distributions are characterized by a pronounced primary peak at mBH ∼ 5 M⊙ and a secondary pile-up at ∼35 M⊙. The latter is primarily driven by the adopted prescriptions for pulsational PISNe (Woosley 2017; Spera & Mapelli 2017), which, at the low metallicity of our models, induce strong mass loss and effectively funnel BH remnants toward masses below mBH ≲ 50 M⊙. In addition, the large primordial binary fractions adopted in our simulations, combined with the preference for short initial orbital separations (Sana et al. 2012), promote efficient envelope stripping and mass loss during binary evolution, further contributing to the accumulation of BHs in this mass range.

Black holes that formed from single stellar collisions broadly follow the overall mass distribution, with a comparatively smaller contribution to the low-mass peak. In contrast, BHs originating from stellar collision chains predominantly populate the ∼35 M⊙ peak. This behavior indicates that, even when stellar progenitors undergo substantial mass growth through repeated collisions, the combined effects of pulsational PISNe and binary interactions efficiently limit the final BH mass. Finally, the number of BHs with mBH ≳ 60 M⊙ is small, marking the lower edge of the PISN mass gap in our models. Nevertheless, we find that the population of BHs in the upper-mass gap grows with cluster mass. In model T13 a single IMBH forms; this case is discussed in detail in Sect. 4.3. In model T16, instead, the enhanced impact of stellar collisions and collision chains leads to the formation of 85 BHs with masses 60 ≤ mBH ≲ 90 M⊙.

Thumbnail: Fig. 9 Refer to the following caption and surrounding text. Fig. 9

Mass distribution of all BHs forming within 20 Myr for models T2, T13, and T16. We show in red the mass distribution of all BHs; the black dashed line shows the BHs formed through single stellar collisions; the black dotted line shows instead the BHs formed through repeated stellar collisions.

4.3 Intermediate-mass black holes

In the following, we consider two families of IMBHs (mBH > 100 M⊙): those that arise from stellar collisions and those that originate from BBH mergers. These two classes are not mutually exclusive, as IMBHs formed through BBH mergers may have stellar progenitors that themselves experienced stellar collisions (Arca Sedda et al. 2023; Arca sedda et al. 2024b; Rantala et al. 2025).

Among the 18 TITANS simulations that reached tsim = 20 Myr, five produce an IMBH through a purely stellar channel. All these IMBHs have masses ≤107 M⊙ and share the same formation pathway: a MS star merges with a cHeb companion, resulting in a star whose helium core remains below the threshold for pulsational PISN (Di Carlo et al. 2020; Kremer et al. 2020a; Arca-Sedda et al. 2021), and that will later collapse into an IMBH. In four models (T12, T13, T14, T18), this merger happens only once and between the components of a primordial binary. Model T9, however, exhibits a more complex evolution: the initial merger product (between a MS and a cHeb star) subsequently undergoes two additional collisions with MS stars (see Fig. 10).

By contrast, none of the VMSs formed in our clusters collapses into an IMBH. Instead, most VMSs either explode as PISNe (Sect. 4.1), leave behind stellar-mass BHs (Sect. 4.2), or are stripped of mass through collisions with compact remnants (Appendix A). This means that, unless substantial material is accreted during BH-star collisions ( fc = 0), and within the model uncertainties, the seeding of IMBHs through the collapse of a VMS is not possible in star clusters with initial ρh ≲ 105 M⊙ pc−3 and trh > 100 Myr.

We also assess the impact of stellar collisions on BBH mergers that produce IMBHs. For this analysis we consider the entire simulated duration of each model, and not just the first 20 Myr. This choice introduces a bias toward simulations evolved for longer times, but it enables us to capture the full role of primordial binaries, stellar collisions, and stellar collision chains across all clusters. Throughout the following discussion, we describe only IMBHs resulting from a single BBH merger, as relativistic kicks are not taken into account in PETAR (see Sect. 4.5). We identify twelve IMBHs formed through BBH mergers, including four produced outside their parent cluster, all with masses ≲140 M⊙. Among these, five have one progenitor produced by a primordial stellar merger, one has a progenitor that experienced repeated collisions, and one has both progenitors originating from primordial stellar mergers. As a consequence, we find that stellar mergers (T10, T12, T17) and mass transfer (T8) in primordial binaries can facilitate the formation of IMBHs through BBH mergers even before core collapse, when it is expected that high cluster densities promote the BBH coupling and eventual merger.

We note that most of the IMBHs formed in our simulations, whether originating from the stellar channel or from BBH mergers, lie within the PISN mass gap (60 ≲ mBH ≲ 120 M⊙; e.g., Banerjee 2022). The high primordial binary fractions adopted in this study enable YMCs to produce IMBHs in, or just above, the PISN range through both stellar interactions and BBH mergers. However, the relatively low escape velocities of these clusters, combined with the absence of mechanisms to grow a seed above the PISN mass range, prevent the formation of substantially larger IMBHs. As a result, our simulations indicate that YMCs with properties consistent with those observed in the local Universe are unlikely to produce large IMBHs even in the presence of high primordial binary fractions.

4.4 Comparison with observations

Very massive stars have been observed in YMCs in the Milky Way (e.g., NGC 3603) and in the Large Magellanic Cloud (e.g., R136, Crowther et al. 2010; Bestenlehner et al. 2020). Owing to its mass (∼5 × 104–105 M⊙) and central density (∼104–105 M⊙ pc−3), R136 provides the closest observational counterpart to our simulated models, although its metallicity (Z ≈ 0.008) is higher than that of our clusters (Z = 0.0002). At least three stars in R136 with age <1.5 Myr have inferred masses in the range 150 ≲ m∗ ≲ 220 M⊙ (Bestenlehner et al. 2020). Among our models, T1 and T9 are the most similar clusters to R136, in terms of mass and density. In T1, no VMSs form within the first 2 Myr, whereas model T9 produces four VMSs with masses 165 ≲ m∗ ≲ 360 M⊙. This comparison confirms that dynamical interactions enhance the formation of VMSs. As discussed in Sect. 3.2, while the total cluster mass Mcl regulates the number of VMSs, the density ρh and relaxation time trh are the key parameters controlling whether repeated stellar collisions can efficiently trigger their formation.

Our models do not form IMBHs with masses comparable to those inferred for IMBH candidates in local star clusters (∼104 M⊙), such as the one reported in ω Centauri (Häberle et al. 2024). As discussed in Sect. 4.3, the VMSs formed in our simulations typically remain within the PISN regime and therefore do not collapse into IMBHs. Alternative formation pathways, such as hierarchical cluster assembly, more extreme initial conditions, or secular growth via tidal disruption of stars, may enhance the seeding and growth of IMBHs (see Sect. 4.5). Moreover, ω Centauri is likely the nuclear star cluster of a disrupted dwarf galaxy (Ibata et al. 2019). In such case, its early dynamical evolution would be remarkably different from that of a globular cluster (e.g., Antonini et al. 2019).

We find that most IMBHs formed in our simulations lie within the PISN mass gap and originate from a single BBH merger, with at least one component already in the gap. Primordial binary evolution, together with repeated stellar collisions, enables the formation of first-generation BHs with masses ≳65 M⊙, which can subsequently merge to produce IMBHs. In this respect, our simulated YMCs naturally reproduce the mass range of BBH mergers associated with GW events such as GW190521 (Abbott et al. 2019) and GW231123 (Abac et al. 2025).

In addition to their impact on BH and IMBH formation, our results have implications for the expected rate of PISNe. Observationally, only a handful of candidate PISNe have been identified to date (e.g., SN 2007bi, Gal-Yam et al. 2009). As a consequence, empirical constraints on the PISN rate are still highly uncertain. In this context, our results suggest that dense, metal-poor star clusters may significantly enhance the production of PISNe compared to isolated stellar populations, by enabling additional stars to enter the PISN mass range through binary mergers and stellar collisions. This effect is particularly relevant at low metallicity, where stellar winds are weaker and massive stars can more easily grow core masses >30 M⊙. If a non-negligible fraction of star formation occurs in such clustered environments, the cosmic PISN rate could therefore be higher than predicted by models based solely on single stellar evolution (Gabrielli et al. 2024; Simonato et al. 2025; Gabrielli et al., in prep.).

Thumbnail: Fig. 10 Refer to the following caption and surrounding text. Fig. 10

Formation pathway of the IMBH in model T9.

4.5 Caveats

The initial conditions adopted in this work (see Table 1) are broadly consistent with observations of YMCs in the local Universe. We also assume high primordial binary fractions (Moe & Di Stefano 2017) and employ a state-of-the-art treatment of single and binary stellar evolution. Nevertheless, our results rely on a number of assumptions that are worth discussing.

First, we assume that all clusters form through the monolithic collapse of a molecular cloud. In reality, star clusters may form from clouds with significant departures from spherical symmetry or through the hierarchical assembly of smaller subclusters (e.g., Sánchez & Alfaro 2009; Ballone et al. 2020; Torniamenti et al. 2021). Recent studies have shown that hierarchical assembly can significantly enhance the rate of repeated stellar collisions, favoring the formation of large (∼103 M⊙) IMBH seeds (Rantala et al. 2024, 2025, 2026). Our results should therefore be regarded as conservative with respect to scenarios in which cluster substructure persists during the early stages of evolution.

Second, the properties of YMCs may evolve with redshift. Recent observations with JWST have begun to reveal extremely compact and massive stellar systems in the early Universe (Claeyssens et al. 2026, and references therein), characterized by high stellar densities and short relaxation times. Such environments would likely experience an enhanced efficiency of stellar collisions and collision chains, potentially allowing the formation of VMSs with helium cores above the PISN regime and increasing the probability of IMBH seed formation. For example, Vergara et al. (2025) recently showed that a cluster with ρh ∼ 108 M⊙ pc−3 and trh ∼ 7 Myr forming in the early Universe can produce a VMS as massive as 50 000 M⊙ collapsing in an IMBH.

Third, our conclusions depend on the adopted prescriptions for stellar and binary evolution. In particular, the stellar mass range associated with the PISN and pulsational PISN regimes is uncertain and model dependent. Alternative prescriptions based on Spera & Mapelli (2017) predict a narrower PISN mass range, which could allow IMBH seeds to form even under initial conditions similar to those explored here. Additional uncertainties arise from the treatment of stellar winds, mass transfer, SN natal kicks, and envelope stripping in massive binaries, all of which can affect the final masses of collision products and compact remnants.

In our simulations, we take the mass of the product of a star–star collision to be equal to the sum of the masses of the two stars involved. This assumption may lead to an overestimation of the mass of the collision product by roughly ∼10–30% (Gaburov et al. 2010; Ballone et al. 2023; Ramírez-Galeano et al. 2025). Moreover, collision products are rejuvenated according to the simplified assumptions reported by Hurley et al. (2002).

We also assume that the fraction of stellar mass accreted by a BH during a collision is zero ( fc = 0). Adopting a higher accretion efficiency can significantly alter the outcomes of stellar-BH collisions and would lead to the formation of up to nine additional IMBHs in our simulations (see Appendix A).

Finally, PETAR does not include relativistic recoil kicks following GW mergers. As a consequence, in Sect. 4.3 we consider only first-generation BBH mergers as viable IMBH formation channels. Given the relatively low escape velocities of our clusters, we expect most successive-generation merger remnants to be ejected. A follow-up study will explore the long-term evolution of compact objects and their merger products in greater detail. Moreover, due to the computational cost of these simulations, we have run here a single realization for each set of initial conditions. We will perform multiple realizations of each model in the future to better quantify the impact of stochasticity on our results.

5 Summary

We have investigated the impact and properties of repeated stellar collisions in YMCs over the first 20 Myr of evolution using the new set of direct N-body simulations TITANS. These models span a wide range of cluster masses and densities consistent with low-redshift YMCs, and include a high, observation-based fraction of primordial binaries (Fig. 1).

We find that the efficiency of repeated stellar collisions increases with increasing half-mass density ρh and decreasing half-mass relaxation time trh (Fig. 3). Clusters with ρh ≳ 104 M⊙ pc−3 and dynamical friction times trh ≲ 300 Myr efficiently form stellar collision chains, although these typically involve only two collisions and are generally triggered by primordial binary mergers. Only a minority of chains in the densest clusters undergo three or more collisions (Fig. 4).

The evolutionary stage of the primary stars involved in repeated collisions depends strongly on the cluster relaxation time (Fig. 5). In particular, the fraction of primaries on the MS decreases from ∼77% in the model with the shortest relaxation time (T1) to ∼17% in the model with the longest relaxation time (T17). The masses of primaries in collision chains vary accordingly, with typical values 50 ≲ m1 ≲ 130 M⊙ in dense clusters with trh ≤ 700 Myr (Fig. 6).

Very massive stars are predominantly formed via primordial binary mergers and typically remain below 300 M⊙. Only in the densest clusters with the shortest relaxation times (T1, T9, T16) do repeated stellar collisions produce VMSs with masses up to ∼400 M⊙ (Fig. 7). As a consequence, the number of VMSs generally scales with cluster mass. Most of these stars either explode as PISNe or leave behind stellar-mass BHs with mBH < 50 M⊙, and none represent viable IMBH seeds.

The number of PISNe in our clusters grows with cluster mass (Fig. 8), as the dominant formation channels involve massive stars and primordial binary mergers, whose frequency increases with cluster mass. Compared to an isolated population of single stars, which is expected to produce PISNe only for zero-age MS masses above 140 M⊙, our clusters generate up to an order of magnitude more PISNe owing to the contribution of stellar collisions.

The BH mass distribution exhibits a primary peak at ∼5 M⊙ and a secondary peak at ∼35 M⊙ (Fig. 9). Stellar collision chains contribute mainly to the secondary peak, indicating that pulsational PISNe and binary interactions effectively regulate the final BH masses, even when stellar progenitors experience substantial mass growth through repeated collisions.

Finally, we find that IMBHs form through two main channels: (i) collisions, either single or repeated, between a cHeb star and a MS star (Fig. 10), producing five IMBHs, and (ii) single BBH mergers, producing twelve IMBHs. In both cases, the resulting IMBHs remain relatively small (mBH < 140 M⊙). Both formation pathways are strongly enhanced by the high primordial binary fraction, which promotes stellar collisions, mass transfer episodes, and the assembly of merging BBHs.

In summary, adopting initial conditions compatible with low-redshift YMCs and state-of-the-art N-body modeling, we find that repeated stellar collisions can be efficient when mass-segregation times are sufficiently short. However, collision chains typically involve only two events and are largely driven by primordial binaries. Within our conservative modeling assumptions ( fc = 0), VMSs do not act as viable IMBH seeds, and only collisions between stars at different evolutionary stages can produce small IMBH progenitors. While this work focuses on the early evolution of YMCs, our simulations are ongoing and have already reached tsim > 1 Gyr for the least dense clusters. A follow-up study will investigate the long-term evolution of compact objects and compact binary mergers in these systems.

Acknowledgements

MM acknowledges financial support from the European Research Council for the ERC Consolidator grant DEMOBLACK, under contract no. 770017. MM also acknowledges financial support from the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC 2181 – 390900948) STRUCTURES. MAS acknowledges funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 101025436 (project GRACE-BH) and from the MERAC Foundation. MB and MAS acknowledge support from the Astrophysics Center for Multi-messenger Studies in Europe (ACME), funded under the European Union’s Horizon Europe Research and Innovation Program, Grant Agreement No. 101131928. SR acknowledges financial support from the Beatriu de Pinós postdoctoral fellowship program under the Ministry of Research and Universities of the Government of Catalonia (Grant Reference No. 2021 BP 00213). The authors acknowledge support by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grants INST 35/1597-1 FUGG and INST 35/1503-1 FUGG. This research made use of NUMPY (Harris et al. 2020), SCIPY (Virtanen et al. 2020), PANDAS (pandas development team 2020). For the plots we used MATPLOTLIB (Hunter 2007). Our simulations made use of the N-body code PETAR (Wang et al. 2020a, and https://github.com/lwang-astro/PeTar) and of the population-synthesis code MOBSE (Mapelli et al. 2017; Giacobbo et al. 2018, and https://gitlab.com/micmap/mobse_open).

References

  1. Abac, A. G., Abouelfettouh, I., Acernese, F., et al. 2025, ApJ, 993, L25 [Google Scholar]
  2. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, ApJ, 882, L24 [Google Scholar]
  3. Adamo, A., Bradley, L. D., Vanzella, E., et al. 2024, Nature, 632, 513 [NASA ADS] [CrossRef] [Google Scholar]
  4. Antonini, F., Gieles, M., & Gualandris, A. 2019, MNRAS, 486, 5008 [NASA ADS] [CrossRef] [Google Scholar]
  5. Arca-Sedda, M., Rizzuto, F. P., Naab, T., et al. 2021, ApJ, 920, 128 [NASA ADS] [CrossRef] [Google Scholar]
  6. Arca Sedda, M., Kamlah, A. W. H., Spurzem, R., et al. 2023, MNRAS, 526, 429 [NASA ADS] [CrossRef] [Google Scholar]
  7. Arca Sedda, M., Kamlah, A. W. H., Spurzem, R., et al. 2024a, MNRAS, 528, 5119 [NASA ADS] [CrossRef] [Google Scholar]
  8. Arca Sedda, M., Kamlah, A. W. H., Spurzem, R., et al. 2024b, MNRAS, 528, 5140 [CrossRef] [Google Scholar]
  9. Askar, A., Szkudlarek, M., Gondek-Rosińska, D., Giersz, M., & Bulik, T. 2017, MNRAS, 464, L36 [CrossRef] [Google Scholar]
  10. Askar, A., Baldassare, V. F., & Mezcua, M. 2024, in Black Holes in the Era of Gravitational-Wave Astronomy, eds. M. Arca Sedda, E. Bortolas, & M. Spera, 149 [Google Scholar]
  11. Ballone, A., Costa, G., Mapelli, M., et al. 2023, MNRAS, 519, 5191 [NASA ADS] [CrossRef] [Google Scholar]
  12. Ballone, A., Mapelli, M., Di Carlo, U. N., et al. 2020, MNRAS, 496, 49 [NASA ADS] [CrossRef] [Google Scholar]
  13. Banerjee, S. 2021, MNRAS, 500, 3002 [Google Scholar]
  14. Banerjee, S. 2022, A&A, 665, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Banerjee, S. 2025, Phys. Rev. D, 112, 063016 [Google Scholar]
  16. Banerjee, S. 2026, arXiv e-prints [arXiv:2602.09694] [Google Scholar]
  17. Barber, J., & Antonini, F. 2025, MNRAS, 538, 639 [Google Scholar]
  18. Bastian, N., Cabrera-Ziri, I., Davies, B., & Larsen, S. S. 2013, MNRAS, 436, 2852 [NASA ADS] [CrossRef] [Google Scholar]
  19. Bestenlehner, J. M., Crowther, P. A., Caballero-Nieves, S. M., et al. 2020, MNRAS, 499, 1918 [Google Scholar]
  20. Binney, J., & Tremaine, S. 2008, Galactic Dynamics, 2nd edn. Bovy, J. 2015, ApJS, 216, 29 [Google Scholar]
  21. Chen, Y., Bressan, A., Girardi, L., et al. 2015, MNRAS, 452, 1068 [Google Scholar]
  22. Claeys, J. S. W., Pols, O. R., Izzard, R. G., Vink, J., & Verbunt, F. W. M. 2014, A&A, 563, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Claeyssens, A., Adamo, A., Kokorev, V., et al. 2026, A&A, submitted [arXiv:2601.16281] [Google Scholar]
  24. Crowther, P. A., Schnurr, O., Hirschi, R., et al. 2010, MNRAS, 408, 731 [Google Scholar]
  25. de Grijs, R. 2009, in Globular Clusters – Guides to Galaxies, eds. T. Richtler & S. Larsen, 353 [Google Scholar]
  26. Di Carlo, U. N., Giacobbo, N., Mapelli, M., et al. 2019, MNRAS, 487, 2947 [NASA ADS] [CrossRef] [Google Scholar]
  27. Di Carlo, U. N., Mapelli, M., Giacobbo, N., et al. 2020, MNRAS, 498, 495 [NASA ADS] [CrossRef] [Google Scholar]
  28. Di Carlo, U. N., Mapelli, M., Pasquato, M., et al. 2021, MNRAS, 507, 5132 [NASA ADS] [CrossRef] [Google Scholar]
  29. Disberg, P., & Mandel, I. 2025, ApJ, 989, L8 [Google Scholar]
  30. Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91 [Google Scholar]
  31. Gabrielli, F., Lapi, A., Boco, L., et al. 2024, MNRAS, 534, 151 [NASA ADS] [CrossRef] [Google Scholar]
  32. Gaburov, E., Gualandris, A., & Portegies Zwart, S. 2008, MNRAS, 384, 376 [NASA ADS] [CrossRef] [Google Scholar]
  33. Gaburov, E., Lombardi, J. C., Jr., & Portegies Zwart, S. 2010, MNRAS, 402, 105 [NASA ADS] [CrossRef] [Google Scholar]
  34. Gal-Yam, A., Mazzali, P., Ofek, E. O., et al. 2009, Nature, 462, 624 [NASA ADS] [CrossRef] [Google Scholar]
  35. Giacobbo, N., & Mapelli, M. 2018, MNRAS, 480, 2011 [Google Scholar]
  36. Giacobbo, N., & Mapelli, M. 2020, ApJ, 891, 141 [NASA ADS] [CrossRef] [Google Scholar]
  37. Giacobbo, N., Mapelli, M., & Spera, M. 2018, MNRAS, 474, 2959 [Google Scholar]
  38. Giersz, M., Leigh, N., Hypki, A., Lützgendorf, N., & Askar, A. 2015, MNRAS, 454, 3150 [Google Scholar]
  39. González, E., Kremer, K., Chatterjee, S., et al. 2021, ApJ, 908, L29 [CrossRef] [Google Scholar]
  40. González Prieto, E., Weatherford, N. C., Fragione, G., Kremer, K., & Rasio, F. A. 2024, ApJ, 969, 29 [CrossRef] [Google Scholar]
  41. Häberle, M., Neumayer, N., Seth, A., et al. 2024, Nature, 631, 285 [CrossRef] [Google Scholar]
  42. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [NASA ADS] [CrossRef] [Google Scholar]
  43. Hernquist, L. 1990, ApJ, 356, 359 [Google Scholar]
  44. Hills, J. G., & Day, C. A. 1976, Astrophys. Lett., 17, 87 [NASA ADS] [Google Scholar]
  45. Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974 [Google Scholar]
  46. Hong, J., Askar, A., Giersz, M., Hypki, A., & Yoon, S.-J. 2020, MNRAS, 498, 4287 [NASA ADS] [CrossRef] [Google Scholar]
  47. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
  48. Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [Google Scholar]
  49. Ibata, R. A., Bellazzini, M., Malhan, K., Martin, N., & Bianchini, P. 2019, Nat. Astron., 3, 667 [Google Scholar]
  50. King, I. R. 1966, AJ, 71, 64 [Google Scholar]
  51. Kıroğlu, F., Kremer, K., Biscoveanu, S., González Prieto, E., & Rasio, F. A. 2025a, ApJ, 979, 237 [Google Scholar]
  52. Kıroğlu, F., Kremer, K., & Rasio, F. A. 2025b, ApJ, 994, L37 [Google Scholar]
  53. Kremer, K., Spera, M., Becker, D., et al. 2020a, ApJ, 903, 45 [NASA ADS] [CrossRef] [Google Scholar]
  54. Kremer, K., Ye, C. S., Rui, N. Z., et al. 2020b, ApJS, 247, 48 [NASA ADS] [CrossRef] [Google Scholar]
  55. Kritos, K., Berti, E., & Silk, J. 2023, Phys. Rev. D, 108, 083012 [NASA ADS] [CrossRef] [Google Scholar]
  56. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  57. Krumholz, M. R., McKee, C. F., & Bland-Hawthorn, J. 2019, ARA&A, 57, 227 [Google Scholar]
  58. Küpper, A. H. W., Maschberger, T., Kroupa, P., & Baumgardt, H. 2011, MNRAS, 417, 2300 [Google Scholar]
  59. Maliszewski, K., Giersz, M., Gondek-Rosinska, D., Askar, A., & Hypki, A. 2022, MNRAS, 514, 5879 [NASA ADS] [CrossRef] [Google Scholar]
  60. Mapelli, M. 2016, MNRAS, 459, 3432 [NASA ADS] [CrossRef] [Google Scholar]
  61. Mapelli, M., Giacobbo, N., Ripamonti, E., & Spera, M. 2017, MNRAS, 472, 2422 [NASA ADS] [CrossRef] [Google Scholar]
  62. Mapelli, M., Spera, M., Montanari, E., et al. 2020, ApJ, 888, 76 [NASA ADS] [CrossRef] [Google Scholar]
  63. Mapelli, M., Zampieri, L., Ripamonti, E., & Bressan, A. 2013, MNRAS, 429, 2298 [NASA ADS] [CrossRef] [Google Scholar]
  64. Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533 [NASA ADS] [Google Scholar]
  65. Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15 [Google Scholar]
  66. Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [Google Scholar]
  67. Oshino, S., Funato, Y., & Makino, J. 2011, PASJ, 63, 881 [NASA ADS] [Google Scholar]
  68. Paiella, L., Arca Sedda, M., Mestichelli, B., & Ugolini, C. 2026, A&A, 708, A200 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. pandas development team, T. 2020, https://doi.org/10.5281/zenodo.20127038 [Google Scholar]
  70. Peters, P. C. 1964, Phys. Rev., 136, 1224 [Google Scholar]
  71. Portegies Zwart, S. F., & McMillan, S. L. W. 2002, ApJ, 576, 899 [Google Scholar]
  72. Portegies Zwart, S. F., Baumgardt, H., Hut, P., Makino, J., & McMillan, S. L. W. 2004, Nature, 428, 724 [Google Scholar]
  73. Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431 [Google Scholar]
  74. Ramírez-Galeano, L., Charbonnel, C., Fragos, T., et al. 2025, A&A, 699, A223 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Rantala, A., Lahén, N., Naab, T., Escobar, G. J., & Iorio, G. 2025, MNRAS, 543, 2130 [Google Scholar]
  76. Rantala, A., Naab, T., & Lahén, N. 2024, MNRAS, 531, 3770 [NASA ADS] [CrossRef] [Google Scholar]
  77. Rantala, A., Naab, T., Lahén, N., et al. 2026, MNRAS, submitted [arXiv:2601.07917] [Google Scholar]
  78. Rastello, S., Mapelli, M., Di Carlo, U. N., et al. 2020, MNRAS, 497, 1563 [CrossRef] [Google Scholar]
  79. Rastello, S., Mapelli, M., Di Carlo, U. N., et al. 2021, MNRAS, 507, 3612 [CrossRef] [Google Scholar]
  80. Rastello, S., Iorio, G., Gieles, M., & Wang, L. 2026, A&A, 707, A217 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Reinoso, B., Klessen, R. S., Schleicher, D., Glover, S. C. O., & Solar, P. 2023, MNRAS, 521, 3553 [CrossRef] [Google Scholar]
  82. Rizzuto, F. P., Naab, T., Spurzem, R., et al. 2021, MNRAS, 501, 5257 [NASA ADS] [CrossRef] [Google Scholar]
  83. Rodriguez, C. L., Zevin, M., Amaro-Seoane, P., et al. 2019, Phys. Rev. D, 100, 043027 [Google Scholar]
  84. Rodriguez, C. L., Weatherford, N. C., Coughlin, S. C., et al. 2022, ApJS, 258, 22 [NASA ADS] [CrossRef] [Google Scholar]
  85. Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444 [Google Scholar]
  86. Sánchez, N., & Alfaro, E. J. 2009, ApJ, 696, 2086 [Google Scholar]
  87. Schrøder, S. L., MacLeod, M., Loeb, A., Vigna-Gómez, A., & Mandel, I. 2020, ApJ, 892, 13 [CrossRef] [Google Scholar]
  88. Sharma, K., & Rodriguez, C. L. 2025, ApJ, 983, 162 [Google Scholar]
  89. Simonato, F., Torniamenti, S., Mapelli, M., et al. 2025, A&A, 703, A215 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Spera, M., & Mapelli, M. 2017, MNRAS, 470, 4739 [Google Scholar]
  91. Spitzer, L. S. 1988, Dynamical Evolution of Globular Clusters [Google Scholar]
  92. Torniamenti, S., Ballone, A., Mapelli, M., et al. 2021, MNRAS, 507, 2253 [NASA ADS] [CrossRef] [Google Scholar]
  93. Vanzella, E., Claeyssens, A., Welch, B., et al. 2023, ApJ, 945, 53 [NASA ADS] [CrossRef] [Google Scholar]
  94. Vergara, M. C., Askar, A., Kamlah, A. W. H., et al. 2025, A&A, 704, A321 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Virtanen, P., Gommers, R., Burovski, E., et al. 2020, https://doi.org/10.5281/zenodo.4406806 [Google Scholar]
  96. Wang, L., Spurzem, R., Aarseth, S., et al. 2015, MNRAS, 450, 4070 [NASA ADS] [CrossRef] [Google Scholar]
  97. Wang, L., Spurzem, R., Aarseth, S., et al. 2016, MNRAS, 458, 1450 [Google Scholar]
  98. Wang, L., Iwasawa, M., Nitadori, K., & Makino, J. 2020a, MNRAS, 497, 536 [NASA ADS] [CrossRef] [Google Scholar]
  99. Wang, L., Nitadori, K., & Makino, J. 2020b, MNRAS, 493, 3398 [NASA ADS] [CrossRef] [Google Scholar]
  100. Wang, L., Fujii, M. S., & Tanikawa, A. 2021, MNRAS, 504, 5778 [NASA ADS] [CrossRef] [Google Scholar]
  101. Woosley, S. E. 2017, ApJ, 836, 244 [Google Scholar]
2

Disberg & Mandel (2025) have recently revised the fit by Hobbs et al. (2005), reporting natal kicks that are a factor of ∼2 lower. This difference has negligible impact on our results, because all the BHs we consider here have negligible natal kicks, as they form from direct collapse.

Appendix A Impact of the chosen accretion fraction fc

In the following we explore the impact of the assumed accretion fraction fc in case of a collision or merger between a star and a BH.

Table 2 shows that nine VMSs in our set of simulations collide with a BH. Because of our assumption of fc = 0, no mass is accreted by a BH in case of a collision with a star. A clear example of this scenario, where the VMS is the result of a stellar collision chain, is shown in Fig. A.1. In most cases, the VMSs form either via binary evolution, accreting mass from the companion which will later become a BH, or via a primordial binary merger. Once the companion is a BH, the star transfers its mass, which is lost, bringing away angular momentum and widening the orbit. The merger finally happens when the star expands enough.

If we assumed fc = 0.5 (e.g., Arca Sedda et al. 2024a), we could obtain nine more IMBHs up to mBH ∼ 215 M⊙ in our simulations, deriving from the merger of binaries formed by a VMS and a BH. We also expect that the accretion of stellar mass from a VMS onto a stellar-mass BH might lead to a significant spin-up of the BH (Kıroğlu et al. 2025a,b).

To evaluate the impact of our prescription on the BH population in our models, we have also looked into all the collisions involving a primary star and a secondary BH. We show in Fig. A.2 the distribution of the mass ratio q = mBH/m∗ for all the simulated models. We see that usually q is smaller than ∼0.4 and that it shows a peak at q ≲ 0.2. Even though the amount of spin-up via stellar collisions is still largely unknown, we can expect that, in case of small mass ratios such as these, both the mass and spin of the BH would increase. Consequently, our adopted choice of fc = 0 effectively represents a conservative scenario, as it minimizes the mass growth of BHs via stellar collisions in our models. Figure A.2 also shows the distribution of BH masses after collisions assuming fc = 0.5. Most BHs in this case have final masses ≲50 M⊙, with a peak around mBH,f ∼ 20 M⊙. The only BHs reaching the IMBH regime originate from collisions in which the primary was a VMS; these cases are listed in Table 2.

Thumbnail: Fig. A.2 Refer to the following caption and surrounding text. Fig. A.2

Mass-ratio q distribution of merging systems with a stellar primary and a BH secondary (top) and final BH mass distribution assuming fc = 0.5 (bottom), for all the simulated models.
Thumbnail: Fig. A.1 Refer to the following caption and surrounding text. Fig. A.1

Stellar collision chain leading to a VMS-BH merger in model T1. In this case fc = 0 and the final BH does not accrete mass from the VMS.

All Tables

Table 1

Initial conditions of the TITANS.

In the text
Table 2

Properties of VMSs in the simulated models.

In the text

All Figures

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Initial binary fraction (top row) and density at half-mass radius (bottom row) as a function of the initial number of stars, for Monte Carlo (green squares), direct N-body (violet circles), and hybrid N-body simulations (brown stars). The TITANS simulations are represented by a red star.
In the text
Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

Schematic representation of a stellar collision chain. We show on the left a chain initiated by the merger of the components of a binary (either primordial or dynamically formed). On the right we show instead a chain initiated by an hyperbolic stellar collision.
In the text
Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

Efficiency of repeated stellar collisions ηrep as a function of the half-mass density ρh, normalized by ρh,3 and Mcl,5. The points are colored as a function of the cluster initial relaxation time trh.
In the text
Thumbnail: Fig. 4 Refer to the following caption and surrounding text. Fig. 4

Number of stellar collision chains Nch,rep as a function of the normalized density at half-mass radius ρh/ρh,3. In violet we represent the number of chains containing two collisions, in blue the number of chains containing three collisions, and in green the number of chains containing four or more collisions. The dashed lines and annotations represent the linear fits of the points.
In the text
Thumbnail: Fig. 5 Refer to the following caption and surrounding text. Fig. 5

Primary stellar type of stars involved in repeated collisions. The x-axis shows the simulated models. The stellar phases on the y-axis are: MS, shell hydrogen burning (sHB), core helium burning (cHeB), shell helium burning (sHeB), and Wolf– Rayet star (WR). The color represents the relaxation time at half-mass radius trh.
In the text
Thumbnail: Fig. 6 Refer to the following caption and surrounding text. Fig. 6

Distribution of primary mass m1 of stars involved in repeated stellar collisions for each simulated model. The color represents the relaxation time at half-mass radius trh. We show with a red star the maximum mass of the product of a stellar collision chain in each model.
In the text
Thumbnail: Fig. 7 Refer to the following caption and surrounding text. Fig. 7

Relation between the maximum VMS mass mVMS,max and the initial density at half-mass radius ρh. The colors used represent the initial relaxation time trh.
In the text
Thumbnail: Fig. 8 Refer to the following caption and surrounding text. Fig. 8

Number of PISNe as a function of the cluster mass Mcl. The filled red points represent the total number; the empty black points represent instead the number of PISNe due to collisions. As a comparison we plot the number of PISNe expected from single stellar evolution NPISN,single.
In the text
Thumbnail: Fig. 9 Refer to the following caption and surrounding text. Fig. 9

Mass distribution of all BHs forming within 20 Myr for models T2, T13, and T16. We show in red the mass distribution of all BHs; the black dashed line shows the BHs formed through single stellar collisions; the black dotted line shows instead the BHs formed through repeated stellar collisions.
In the text
Thumbnail: Fig. 10 Refer to the following caption and surrounding text. Fig. 10

Formation pathway of the IMBH in model T9.
In the text
Thumbnail: Fig. A.2 Refer to the following caption and surrounding text. Fig. A.2

Mass-ratio q distribution of merging systems with a stellar primary and a BH secondary (top) and final BH mass distribution assuming fc = 0.5 (bottom), for all the simulated models.
In the text
Thumbnail: Fig. A.1 Refer to the following caption and surrounding text. Fig. A.1

Stellar collision chain leading to a VMS-BH merger in model T1. In this case fc = 0 and the final BH does not accrete mass from the VMS.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.