Intermediate-mass Black Hole Progenitors from Stellar Collisions in Dense Star Clusters

Very massive stars (VMSs) formed via a sequence of stellar collisions in dense star clusters have been proposed as the progenitors of massive black hole seeds. VMSs could indeed collapse to form intermediate-mass black holes, which would then grow by accretion to become the supermassive black holes observed at the centers of galaxies and powering high-redshift quasars. Previous studies have investigated how different cluster initial conditions affect the formation of a VMS, including mass segregation, stellar collisions, and binaries, among others. In this study, we investigate the growth of VMSs with a new grid of Cluster Monte Carlo star cluster simulations—the most expansive to date. The simulations span a wide range of initial conditions, varying the number of stars, cluster density, stellar initial mass function (IMF), and primordial binary fraction. We find a gradual shift in the mass of the most massive collision product across the parameter space; in particular, denser clusters born with top-heavy IMFs provide strong collisional regimes that form VMSs with masses easily exceeding 1000 M ⊙. Our results are used to derive a fitting formula that can predict the typical mass of a VMS formed as a function of the star cluster properties. Additionally, we study the stochasticity of this process and derive a statistical distribution for the mass of the VMS formed in one of our models, recomputing the model 50 times with different initial random seeds.


Introduction
Although the dynamical evolution of dense star clusters has been studied extensively for decades, the details of the cluster formation stages and their initial conditions remain highly uncertain.Efforts to tackle these open questions are taking place both theoretically and observationally.In particular, the latest cosmological simulations are approaching the resolution necessary to robustly track the formation of bound clusters in a range of galaxy types and redshifts (e.g., Ma et al. 2020;Grudić et al. 2023;Rodriguez et al. 2023).Even so, resolving cluster formation in these large cosmological simulations remains challenging because of their multiscale nature.On the observational side, the James Webb Space Telescope (JWST) has opened a new window into star cluster formation, with several studies reporting observations of candidate protoglobular clusters at high redshifts (e.g., Mowla et al. 2022;Vanzella et al. 2023).
Globular clusters (GCs), being some of the densest environments in the Universe, host numerous exotic objects and transient phenomena arising from strong dynamical interactions, including direct physical collisions.Previous studies have shown that young star clusters, the likely progenitors of GCs, may produce stars with masses greatly exceeding the maximum mass in the stellar initial mass function (IMF) through successive stellar collisions and mergers (e.g., Sanders 1970;Lee 1987;Quinlan & Shapiro 1990;Ebisuzaki et al. 2001;Portegies Zwart & McMillan 2002).These so-called very massive stars (VMSs) have been the focus of considerable previous theoretical work because they are natural progenitors for intermediate-mass black holes (IMBHs).
The collisional process to form a VMS begins immediately after the segregation of the most massive stars deep into the core of the cluster.Due to the Spitzer instability (Vishniac 1978), these stars cannot achieve energy equipartition, and the core develops a high-velocity dispersion, which promotes stellar collisions.For sufficiently dense clusters, growth will start as a result of stellar collisions and mergers.As the mass of a merger product grows, so does its collision cross section, resulting in even more collisions.This results in a positive feedback loop that can rapidly produce a VMS of hundreds to thousands of solar masses.
To study VMS formation in star clusters, Gürkan et al. (2004) investigated the effects of mass segregation and core collapse.This was accomplished through the implementation of Monte Carlo simulations in systems where parameters such as the cluster density profile, stellar IMF, and initial star count were systematically varied.This study found that the mass of the collapsing core was always close to ∼10 −3 times that of the total cluster mass.Remarkably, this follows the observed correlation between central BH mass and total host mass in many astrophysical environments (Ferrarese & Merritt 2000).Note, however, that Gürkan et al. (2004) did not account for the effects of stellar evolution.Freitag et al. (2006aFreitag et al. ( , 2006b) performed the first cluster simulations that included precise treatment of stellar collisions and followed the evolution of the cluster and formation of a collisional runaway.Particularly, Freitag et al. (2006b) studied runaway collisions in young star clusters by varying physical parameters such as cluster mass, size, and initial concentration.
These studies confirmed that when the core collapse timescale is shorter than the main-sequence (MS) evolution timescale for the most massive stars (t ≈ 3 Myr), the cluster will undergo a collisional runaway.However, these studies did not incorporate the role of binaries in the runaway process, which have an important role in the evolution of the cores of star clusters.
VMSs are often assumed to be progenitors of IMBHs (e.g., Ebisuzaki et al. 2001;Portegies Zwart et al. 2004;Gürkan et al. 2006;Giersz et al. 2015;Mapelli 2016).In an early study, Ebisuzaki et al. (2001) introduced the collisional runaway formation scenario for IMBHs and discussed the possibility that these IMBHs will eventually sink to the Galactic center and be the seeds for supermassive BHs (SMBHs).The possibility that massive collision products could avoid the pair-instability regime and directly collapse into a massive BH (e.g., Di Carlo et al. 2019;Spera et al. 2019;Di Carlo et al. 2020) has recently been confirmed via hydrodynamic simulations of stellar collisions and the product's ensuing evolution (Costa et al. 2022;Ballone et al. 2023).The possible presence of IMBHs at the centers of GCs has also been studied for many years (see Greene et al. 2020 for a review).Tentative evidence of massive BHs at the cores of GCs includes velocity dispersion signatures in nearby GCs (e.g., Noyola et al. 2010;Jalali et al. 2012;Baumgardt 2017), accretion signatures from radio observations (e.g., Maccarone 2004;Paduano et al. 2024), hypervelocity stars (e.g., Edelmann et al. 2005;Gualandris & Portegies Zwart 2007), observations of ultraluminous X-ray sources (e.g., Colbert & Mushotzky 1999;Farrell et al. 2009), and pulsar acceleration measurements (Kızıltan et al. 2017).
The likelihood of a cluster forming a VMS depends on various physical properties, among which is the IMF, which remains poorly constrained to this day, especially at high stellar masses.Although many studies of GCs assume a canonical Kroupa IMF (Kroupa 2001), observations suggest that it may not be universal (e.g., De Marchi et al. 2007;Bartko et al. 2010;Haghi et al. 2017;Wirth et al. 2022).Furthermore, several studies have shown that the IMF strongly impacts the dynamical evolution and survival of GCs (e.g., Chernoff & Weinberg 1990;Chatterjee et al. 2017;Giersz et al. 2019;Weatherford et al. 2021).
In particular, Weatherford et al. (2021) explored the impact of the slope of the IMF (at the high-mass end) on the compact object population.This study found that in addition to producing more BHs, clusters with a top-heavy IMF also produce substantially more binary BH (BBH) mergers, especially those involving (or resulting in) production of upper-mass-gap BHs (e.g., Spera & Mapelli 2017;Takahashi et al. 2018;Farmer et al. 2019;Marchant et al. 2019) and IMBHs.The latter is due to three factors: top-heavy IMFs produce heavier stars and therefore heavier BHs, but also lead to several times more stellar collisions-due to scaling of stellar radii and gravitational focusing with mass-and more hierarchical mergers.
Another physical parameter that influences the rate of dynamical interactions in star clusters is the primordial binary fraction (e.g., Heggie & Hut 2003;Fregeau & Rasio 2007;Chatterjee et al. 2010).Since binaries have a larger interaction cross section than single stars, they offer a larger effective area for encounters to take place.As shown in previous studies by González et al. (2021) and González Prieto et al. (2022), increasing the binary fraction for high-mass stars (M > 15M e ) to 100%, more in line with observed binary fractions in the Galactic field (e.g., Sana et al. 2012;Moe & Di Stefano 2017), dramatically increases the number of massive stellar collisions and thus results in more massive BHs.
In this paper, we re-examine the formation of VMSs while self-consistently modeling stellar and binary evolution.We systematically cover the parameter space, extending boundaries in cluster size, density, and mass.Furthermore, we fully explore the stochasticity of this process and derive statistical distributions for the masses of the VMSs.In Section 2, we describe the physical prescriptions and parameters varied in this study.In Section 3, we analyze the formation of VMSs in our models, while Section 4 presents a simple equation to estimate the most massive star formed through collisions in a cluster.We discuss the resulting BH population in Section 5 and present a statistical study of our models in Section 6.Finally, in Section 7, we discuss the implications and caveats of this study.

Cluster Simulations
We perform our simulations using Cluster Monte Carlo (CMC), a Hénon-type Monte Carlo code that models the evolution of star clusters (see Rodriguez et al. 2022 for the most recent overview).CMC incorporates prescriptions for various physical processes such as two-body relaxation (Joshi et al. 2000), treatment for stellar collisions (Fregeau & Rasio 2007), and direct integration of small N-body strong encounters using Fewbody (Fregeau et al. 2004).Finally, the population synthesis code COSMIC is fully integrated into CMC to treat stellar and binary evolution (Breivik et al. 2020).
We run a set of 324 simulations (listed fully in Table 1) that systematically investigate a broad spectrum of initial cluster properties.First, the grid varies the initial number of objects in the cluster-both singles and binaries-in the range N3 = (4, 8, 16, 32) × 10 5 .Second, to examine the impact of cluster density, we vary the cluster's initial virial radius r v /pc = (0.5, 1, 2).Both the values for N and r v are motivated by earlier work demonstrating that clusters with these initial conditions evolve into GCs similar to those observed in the Milky Way (Kremer et al. 2020b;Rui et al. 2021).
While past studies have explored the role of the IMF and binary fraction independently, the present work examines their combined effect on the cluster.We assume a typical primordial binary fraction of f b = 0.05 for stars born less massive than 15 M e and vary the high-mass binary fraction f b,high = (0.05, 0.25, 1.0) for stars born more massive than 15 M e .We sample primary stellar masses from the Kroupa (2001) To vary the IMF, we choose three different values for α 3 = (1.6,2.3, 3.0), corresponding to the approximate 95% confidence interval around the canonical value, α 3 = 2.3  1.63 (Kroupa 2001).For each set of initial conditions, we run three statistically independent realizations of the same cluster.See Section 6 for a detailed discussion of the number of realizations necessary to resolve key behavior.

Physical Prescriptions
Modified Radii Prescriptions for Massive Stars: Stellar evolution is an active area of research, with many uncertainties, especially in the high-mass regime.Agrawal et al. (2020) studied the uncertainties in massive stellar evolution models by comparing different stellar evolution codes, finding a notable disparity (see their Figure 8) between the current extrapolation of maximum radius for massive stars in the Single Stellar Evolution code (SSE) and Modules for Experiments in Stellar Astrophysics (MESA).In particular, for stars with a zero-age main sequence (ZAMS) mass > 40 M e , the predicted maximum stellar radius in SSE (used in COSMIC) is 1 order of magnitude higher than the one predicted by more detailed stellar evolution models such as MESA.Note.Columns (2)-( 6) list the initial physical parameters of our clusters, including the initial number of objects, virial radius, the absolute value of the high-mass stellar IMF slope (α 3 ), high-mass binary fraction, and cluster mass.Columns (7)-( 9) list the mass of the most massive star formed in each realization of the same model.The dagger indicates masses that result from stellar IMF alone (i.e., those that do not experience any collisional growth).Column (10) lists the mass of the most massive star as predicted by the fitting formula described in Section 4, with the error bars indicating the 95% confidence interval.Finally, column (11) records the number of binary-single and binary-binary interactions involving the merger of at least two stars that are both more massive than 15 M e , normalized by the initial number of massive stars sampled in the cluster.This gives a rough sense of the number of collisions per massive star in each model.
To correct for this likely overestimation of the stellar radius, we have truncated the radius of any star with a ZAMS mass M 40 M e to a maximum value of 10 3 R e .So if a star in our simulations reaches a stage in its evolution where it is assigned a stellar radius >10 3 R e , we simply rescale the radius-and the radii of the core and convective envelope, proportionately-to this prescribed limit.While more precise extrapolations of stellar radii are currently under investigation, this provisional change prevents an artificially large collision cross section, thereby ensuring a more accurate collision rate.We have rigorously tested these new prescriptions on thousands of stars, confirming that it does not alter their stellar evolution from default SSE assumptions.
Stellar Collision Products: The properties of a collision product depend on the details of the collision and internal structure of the stars.Due to the large uncertainties in isolated high-mass stellar evolution-let alone the hydrodynamic complexities of postcollision evolution-we adopt the conservative choice of setting the total mass of the collision product M 3 equal to the sum of the masses of the colliding stars (M 3 = M 1 + M 2 ) for collisions involving two MS stars.This assumption of mass conservation is motivated by hydrodynamic simulations of stellar collisions in globular cluster-like environments (e.g., Lombardi et al. 1996;Sills et al. 2001;Costa et al. 2022;Ballone et al. 2023).In the case of a collision between a giant star and an MS star, we assume that the resulting object has the core of the giant star (M c3 ) embedded in the envelope of both stars The product of all stellar collisions must be rejuvenated since new gas is introduced into the envelope and potentially the core of the new star-giving opportunity to burn more fuel.We assign the rejuvenated effective age of the merger product to be where (t MS1 , t MS2 , t MS3 ) are the MS lifetimes of the two collision components and the collision product while (t 1 , t 2 ) are the stellar ages of the two collision components.f rejuv is a coefficient that determines the level of rejuvenation experienced by the collision product.We adopt f rejuv = 1 by default and refer the reader to Breivik et al. (2020) for a discussion of these rejuvenation prescriptions as well as the choice for f rejuv .

VMS Formation
To study how the formation of a VMS depends on the initial conditions of a star cluster, we closely examine the formation process of the most massive star in each cluster, denoted as M ,max


. Figure 1 shows the mean mass of the most massive star formed through stellar collisions across the three realizations performed at each point in the model grid.The plot reveals a consistent trend: as the number of initial objects increases and the cluster becomes more compact (indicated by a smaller r v value), M ,max  also rises.Furthermore, clusters born with a topheavy IMF (α 3 = 1.6) feature a collision rate in the first 10 Myr that is ≈4 times higher than those born with a canonical IMF (for typical GCs born with N = 8 × 10 5 ).As a consequence, for a given combination of N and r v within the grid, a more topheavy stellar IMF (smaller α 3 ) results in higher values of M ,max  , since a higher number of stellar collisions facilitates the growth of the VMS.
A comparable correlation occurs in the case of the high-mass binary fraction, where a higher primordial fraction of massive binaries increases the mass of the VMS.This is a result of an increased rate of massive star collisions due to the presence of massive binaries in the cluster, which ultimately facilitate the formation of a more massive star.Furthermore, since massive stars are rare compared to low-mass stars, increasing the highmass binary fraction does not significantly increase the total binary fraction.Consequently, clusters will not experience significant heating from the addition of these primordial massive binaries alone and the process to form the VMS can proceed uninterrupted.The trend is less pronounced for lower-N runs and lower-density clusters, which tend to yield more diffuse clusters where the collisional rate is reduced.Overall, it is evident that the slope of the IMF and the virial radius have the strongest effect on the collisional formation of a VMS in a star cluster.
By examining each section within the parameter space more closely, we can learn more about the physical processes driving the formation of the massive star.Focusing on models with a virial radius of 2 pc, we find that M ,max  never exceeds 400 M e across all initial conditions.The formation of the most massive star typically involves a few stellar interactions, or in most cases, one binary-single or binary-binary interaction resulting in the collision of more than two stars.These less-concentrated models, unlike their denser counterparts, do not have as strong of a correlation between M ,max  and the slope of the IMF or the high-mass binary fraction.This can be explained by massive binaries taking a longer time to segregate toward the cluster center, thus limiting their ability to significantly increase the collision rate and trigger a runaway process.
Models with a virial radius of 1 pc are initially more dense, allowing us to begin observing the onset of the formation of a VMS via a sequence of stellar collisions.The trend becomes apparent as the slope of the IMF becomes shallower, which is equivalent to a leftward movement within each 3 × 3 box in Figure 1.Most notably, in models with α 3 = (1.6,2.3) and f b,high = 1, the formation of M ,max  primarily occurs through a few collisions (typically between 2 and 3) that involve massive stars.In contrast with the formation channel described for the 2 pc models, the stars involved in these collisions are slightly more massive than those initially sampled from the IMF.These unusually massive stars tend to be products of binary coalescences, which increase as the primordial fraction of binaries increases.When combined with a slightly denser cluster, this mechanism facilitates the formation of a VMS.
Finally, models with initial virial radii of 0.5 pc represent the densest clusters in the grid.Within this subset, we observe a greater diversity of pathways leading to the formation of massive stars.In models with a top-light IMF (α 3 = 3.0), the progenitors of massive stars experience multiple collisions both during their MS stage and their giant phase.These collisions typically result in a mass gain of approximately 16 M e per collision.Most importantly, these collisions tend to occur at later cluster ages, roughly t > 4 Myr.For models with a canonical IMF, the average mass gain is ≈37 M e per collision.A detailed analysis of the collisional histories of M ,max  reveals that the significant (and fast) mass accumulation primarily arises from a series of massive binary-mediated interactions resulting in the merger of more than two stars.This enables a very rapid increase in mass and forms a single VMS on timescales shorter than 4 Myr.Finally, in models with a top-heavy IMF (α 3 = 1.6), we typically observe more than one VMS forming on a timescale <3 Myr.In these clusters, multiple massive stars promptly scatter with each other in the center of the cluster and merge, triggering the start of a more extreme runaway process.The average mass gained per collision in these models is ≈71 M e .
In Figure 2 we show the collisional history for M ,max  across models spanning from the least dense and least massive to the most dense and most massive (represented by the four corners of Figure 1).Among each set of three realizations, we selected the one yielding the median VMS mass as a representative case.From the figure, we see that the star that experiences the most growth is not necessarily the most massive star initially sampled from the stellar IMF (shown by the initial mass of the VMS in model l8, which is below the IMF's upper bound).Furthermore, the star in model c8 (indicated by a dot) does not experience any growth via collisions.The star in model l8 also does not experience significant growth, as that model only differs from c8 in the initial number of objects.As we move toward denser and initially more top-heavy IMF models shown in yellow and blue, we see the formation of a VMS via successive stellar collisions.Particularly, we begin to see the exponential growth of a star in <3 Myr.
To further illustrate how the channels for forming M ,max  vary across different regions of the parameter space, Figure 3 depicts the cumulative distribution of the total number of stellar collisions leading to the formation of M ,max  in each model.In the leftmost panel, models with lower r v consistently exhibit a Figure 1.The mean mass of the most massive star formed in a cluster based on its initial conditions.Specifically, the horizontal axis specifies the cluster's initial virial radius r v , and the vertical axis specifies its initial number of objects N (both singles and binaries).Within each box in the r v -N grid, we present a 3 × 3 subgrid to distinguish models with different high-mass stellar IMF slope α 3 = (1.6,2.3, 3.0), from left to right, and high-mass binary fraction f b,high = (0.05, 0.25, 1.0), from bottom to top.The size and color of the circles reflect the mass of the most massive star formed at each set of initial conditions (averaged over all three realizations).
richer dynamical history, aligning with the expectations that denser clusters will experience higher collision rates.The middle panel reveals a subtler trend, but the overarching message remains that a top-heavy IMF enhances the collision rate.In the rightmost panel, although the trend appears less distinct for varying binary fraction values (as expected from the results shown in Figure 1), we still observe that a higher binary fraction enhances the collision rate.It is important to emphasize that we only fix one physical parameter per panel, so the range of the total number of stellar collisions exhibited by each line is a consequence of the diverse formation channels across the entire grid.
To gain a more detailed understanding of the dynamics within the core of the clusters that form a VMS, we closely analyze in the top panel of Figure 4 the evolution of the Lagrange radii (enclosing the specified percentages of the cluster's mass).We pay particular attention to the evolution of objects within the 1% Lagrange radius in models a1, a4, and a7, which differ only in α 3 (and thereby cluster mass).Notably, the behavior of the innermost particles demonstrates that a topheavy IMF results in a more pronounced core contraction, which facilitates the formation of a VMS.Although we show only the α 3 -dependence of the Lagrange radii evolution in Figure 4 (with a fixed value for N, r v , and f b,high ), similar trends exist when comparing models with different N, r v , or f b,high .In particular, a deeper collapse is typically seen for models with smaller r v .
It can also be seen in Figure 4 that in models where a very massive star forms, the initial core collapse occurs on a timescale shorter than the MS lifetime for the most massive stars (t å ≈ 3 Myr, marked as a vertical dashed line).This agrees with previous studies that found a runaway process only occurs when t cc < t å (e.g., Freitag et al. 2006a).After 3 Myr, mass loss due to supernovae and corresponding remnant ejections from natal kicks halts core collapse.This leads to a re-expansion of the core, resulting in a decrease in density, effectively preventing a runaway process.To account for this, we have included the time of the first supernova for each cluster in Figure 4 (shown as vertical dotted lines).Except for the model with α 3 = 1.6, the core begins to re-expand as the first supernova event occurs, averting an extreme collisional runaway.
For the model with α 3 = 1.6 (and many dense clusters), the initial re-expansion is due to a combination of factors.First, as the core contracts, many objects sink toward the center of the cluster.This increase in central density forms new binaries through three-body binary formation, which in turn heats the core (Cohn et al. 1989).Furthermore, the collisions that formed the massive object have extracted some of the (negative) potential energy of the cluster.Thus, once the series of collisions stops, the cluster begins to re-expand.In cases where a very massive star does not form, like model a4 with a canonical stellar IMF (purple curve), it is a combination of mass loss due to supernovae and three-body binary formation that drives the core re-expansion.
To demonstrate how the initial core collapse contributes to the formation of a VMS, we present the collisional history of M ,max  in run a1 in the lower panel of Figure 4.Each distinct colored curve within the diagram represents a separate branch in the evolutionary process.At about 2.3 Myr, three separate massive stars form, each weighing approximately 700 M e .This occurs precisely when the core undergoes its initial contraction.The stars formed in the yellow and magenta pathways merge, giving rise to a star of roughly 1000 M e .While the core begins to re-expand, the merger product undergoes gradual mass loss due to stellar winds.At approximately 4.5 Myr, this star collides with a star of mass 400 M e , resulting in a collision product with a mass of approximately 1100 M e .
It is important to note that at any given time in a cluster, more than one VMS might be present, which is not depicted in the figures in this paper.Our study focuses on assessing the extent of runaway phenomena occurring within the initial 10 Myr of the cluster's lifetime.This takes into account whether or not the massive stars formed in the cluster have enough time to sink to the cluster center due to dynamical friction and merge.
In general, although VMSs form during the initial core collapse (core contraction) in our cluster simulations, the runaway process halts once the core re-expands.As a result, we are not in a regime (unlike previous studies, e.g., Freitag et al. 2006aFreitag et al. , 2006b) ) where an extreme collisional runaway scenario causes the entire cluster core to collapse and form a VMS with a mass ∼10 −3 times the cluster mass.Instead, we observe smaller-scale runaways that enable the cluster to re-expand and continue its evolution.This can be primarily attributed to the delay in the core collapse timescale because of updates in stellar evolution prescriptions and the role of binaries in heating up the cluster core.Consequently, none of our models reach core collapse during the initial 10 Myr.The future fate of such clusters (beyond the 10 Myr modeled in this paper) falls outside the scope of our current study but will be investigated in a subsequent publication.

Fitting Formulae
Carefully mapping the different evolutionary outcomes of clusters across a broad physical spectrum requires an extensive, high-resolution grid of simulations.However, this task is rather computationally impractical, so we develop a simple fitting formula that can be used to estimate M ,max  for a cluster, based upon the model grid explored in this paper.We begin with a simple power-law formula for the dependence of the maximum Figure 2. The collision history for the most massive star formed in models c8, a0, l8, and j2, listed in Table 1.The number following the underscore indicates the realization (chosen to be the one that yields the median VMS mass).The massive star formed in model c8 (shown as a dot) does not experience any growth.stellar mass on N, r v , α 3 , and f b,high : To determine the values of the coefficient A and the power-law exponents, we perform a Bayesian inference technique known as nested sampling (see Skilling 2006, for a review of the method), computing the parameters constraining the model using the nestle package (Mukherjee et al. 2006;Shaw et al. 2007;Feroz et al. 2009).This method restricts mass priors by sampling within likelihood contours, demonstrating excellent efficacy in high-dimensional parameter estimation.We refer the reader to the Appendix to learn more about the performance of the fit.
We use a uniform prior from 10 −3 to 10 3 for A and from 10 −3 to 10 1 for each of the exponents.We find that the data is best described by the values in Table 2, which shows that the parameters that have the most significant effect on the estimation of M ,max  are r v and α 3 .This is expected since a smaller virial radius leads to a higher overall collision rate and a shorter mass segregation timescale for the massive stars, promoting earlier and more frequent massive collisions.Upper panel: the time evolution of the Lagrange radii-from top to bottom enclosing the 99%, 50%, 10%, and 1% of the cluster's total mass-for simulations a1, a4, and a7 in Table 1.These models all have N = 4 × 10 5 , r v = 0.5 pc, and f b,high = 0.25.We vary α 3 = (1.6,2.3, 3.0) shown in black, purple, and green curves, respectively.The dashed vertical line represents the main-sequence lifetime for the most massive stars in the cluster (t = 3 Myr).The dotted vertical lines show the time of the first supernova in each cluster model.Lower panel: the collisional history of the three most massive bodies in the simulation with a top-heavy IMF (α 3 = 1.6), distinguished by color.Two of these massive stars (blue and purple) merge at t ≈ 2.3 Myr, and the remnant (colored blue) merges with the third massive star (yellow) at t ≈ 4.4 Myr to form the final VMS.Note.The best-fit values for the fitting formulae shown in Equation (3) obtained using the nested sampling method (Skilling 2006).With this model to predict the onset of VMS formation, we aim to extend predictions across a broader range of cluster initial conditions.Figure 5 shows the predicted maximum stellar-mass distribution across models of diverse f b,high , r v , and α 3 parameters.To obtain predictive model outcomes, we extract 1000 samples for every data point across the parameter space by randomly drawing from the parameter distribution listed in Table 2. Subsequently, we calculate the mean of the 1000 samples as well as the 95% credibility region, illustrated in Figure 5 as a solid line and shaded region, respectively.
To qualitatively demonstrate the performance, we overplot the data obtained from the CMC simulations, revealing strong agreement as the majority of the data falls within the 95% confidence interval.Discrepancies between the model and data occur in regions where the initiation of the collisional runaway becomes highly stochastic, leading to substantial variance in the mass of the most massive star.This is the case for the model with N = 32 × 10 5 , r v = 0.5 pc, α 3 = 1.6, and f b,high = 1 (model j2 in Table 1), where the stochasticity is apparent in the three data points for M ,max  .
To assess the accuracy of the predictive model, we computed additional simulations with r v = 4 pc, shown in Figure 5.These simulations were not utilized in the parameter estimation process for determining the best-fit model.Instead, they are only used to demonstrate the performance of the model.We see that the mean masses from our CMC simulations fall within the 95% confidence interval of the predictive model.We also note that the model underperforms in some clusters.This discrepancy arises from the fact that in our simulations, if the cluster does not experience a significant number of collisions, the most massive star at any given time typically corresponds to the most massive star initially sampled from the Kroupa IMF, often around 150 M e .Thus, it is to be expected that CMC models with low densities and top-light IMFs will not form collisional products with masses much higher than the IMF upper limit.This is shown by the flat evolution of the maximum stellar mass as a function of virial radius of the CMC data plotted in Figure 5. Currently, our fitting formulae do not take into account the assumed IMF, so the prediction for the maximum mass is allowed to go below the IMF upper limit.
It is important to note that this fitting formula is intrinsic to our assumed stellar evolution prescriptions.As such, it may not hold for clusters with different stellar treatments and physical assumptions (e.g., an Elson profile), or for clusters that deviate significantly from the parameter space covered by our model suite (e.g., very low-mass or very high-mass clusters).Moreover, in the denser clusters, the mechanism through which a massive star forms is stochastic, so we expect considerable variance in mass within those regimes.

Massive BH Formation
An essential question stemming from this research is the fate of the massive stars in these clusters, and whether the formation Figure 6.Normalized BH mass distribution for all models listed in Table 1.Models with f b,high = 1 are depicted in black, those with f b,high = 0.25 in magenta, and models with f b,high = 0.05 in green.The solid lines represent models with α 3 = 1.6, dashed lines denote models with α 3 = 2.3, and dotted lines correspond to models with α 3 = 3.0.The blue-shaded region indicates the upper mass gap (defined here between 40 and 120 M e ), while an arrow marks the start of the IMBH regime (M > 120 M e ).
of the runaway object has an impact on the compact object population.In Figure 6, we show the spectrum of BH masses formed across the first 10 Myr of our cluster models.For all values of N, as r v decreases (moving left in each row), there's a notable increase in the number of BHs formed within or beyond the upper mass gap (assumed here to be between 40 and 120M e , but boundaries are uncertain; e.g., Spera & Mapelli 2017;Takahashi et al. 2018;Farmer et al. 2019).This agrees with the findings of González Prieto et al. (2022), who showed that denser clusters exhibit higher rates of stellar collisions, facilitating the formation of BHs in the upper-massgap and IMBH regimes.
When increasing the initial number of objects while maintaining a constant r v (moving downward in each column), we observe that the total number of massive BHs formed increases for dense models.This is due to the increased "mass budget" as the number of initial objects grows, which allows more stars that were not previously massive BH progenitors to merge and populate the massive BH region.This is also apparent in the decrease in BHs with masses in the range 4-10 M e .Furthermore, within each panel, a higher binary fraction often correlates with the formation of more massive BHs alongside a lower α 3 value.This is also consistent with the general trends we observe in Figure 1.
Crucially, these BH spectra solely represent the initial 10 Myr of the cluster's evolution and do not constitute a complete sample of the full BH population.In fact, due to the rejuvenation prescriptions outlined in Section 2 and lower-mass stars, there are BH progenitors left in most clusters at that time.We defer more detailed analyses of the long-term population and retention of compact objects to future studies.Nonetheless, the presence of such a diverse population of massive compact objects in these clusters hints at the possibility they could significantly contribute to numerous gravitational-wave events (Rodriguez et al. 2019;Kremer et al. 2020a;González et al. 2021;Weatherford et al. 2021;González Prieto et al. 2022).

CMC Statistical Analysis
Given the inherently statistical nature of the Monte Carlo algorithm, we now investigate whether simulations with identical macroscopic initial conditions but different random seeds (which set the exact initial stellar positions and velocities) yield a statistical mass distribution for M ,max  .In Figure 7, we present the distribution of stellar masses from a set of 50 simulations with an initial population of 16 × 10 5 objects, a virial radius of 1 pc, α 3 = 1.6, and a high-mass binary fraction of 1.For qualitative comparison, we also overplot the first three runs in orange.These specific initial conditions were chosen because they represent one of the densest regions in our current grid, often resulting in the formation of a VMS with a mass of a few hundred M e .While these models are to some extent stochastic in nature, an examination of Figure 7 reveals that the stellar masses roughly follow a Gaussian distribution centered at 440 M e , with a spread of 50 M e (the Gaussian is added for illustrative purposes, using the mean and standard deviation of the data).While the distribution appears to cover a broad range of stellar masses, it is much narrower than the spread of masses shown across the entire grid explored in this study (see Figure 1).
In all 50 realizations, the formation of the massive object results from a series of stellar collisions occurring during binary-single and binary-binary interactions.Figure 8 illustrates the cumulative distribution of the total number of collisions that contributed to the formation of M ,max  in each of the 50 cluster models.To emphasize the collisions that significantly contribute to mass buildup, we show in orange those collisions where the colliding star is more massive than 15 M e .As depicted, for most runs, the massive star forms after a few (≈3-5) stellar collisions.In numerous cases, more than one star merges during a binary-mediated interaction.To account for this, we also plot in black the number of interactions that resulted in collisions.It is evident that the number of interactions follows a narrow distribution with a mean of approximately 2.7 and a standard deviation of 1.68.Thus, across all 50 simulations, there is a high level of consistency in the number of interactions that M ,max  undergoes.This marks the first time we have been able to characterize the realization-to-realization variability of CMC models with high resolution, even in the densest and most stochastic regimes.The consistent agreement in the low number of collisions required for the formation of a VMS  emphasizes the importance of studies using precise hydrodynamic simulations to understand and model the properties of collision products (Costa et al. 2022;Ballone et al. 2023).

Summary
This paper presents findings from an extensive grid of CMC simulations tracking the dynamical evolution of star clusters during the first 10 Myr of their lifetime across a spectrum of initial conditions.We particularly focus on the formation process of VMSs, which are likely progenitors of IMBHs.The results from this study can be condensed into three principal findings: 1. Clusters that start with sufficiently high densities experience a phase of core contraction at early times.If this contraction precedes the first supernovae in the cluster, it leads to the formation of a VMS through a collisional runaway instability.In order of importance (see Table 2), the maximum mass reached depends strongly on the high-mass slope of the stellar IMF, the initial cluster density, and the high-mass primordial binary fraction.
2. We have derived a fitting formula that can be used to estimate the mass of the VMS as a function of initial cluster conditions.While this equation depends on specific assumptions regarding stellar evolution and collision prescriptions, it serves as a useful tool to evaluate the potential for a collisional runaway before performing computationally expensive N-body simulations.
3. At the end of our simulations, some of the VMSs have collapsed to form a BH in the upper mass gap or an IMBH.These massive BHs will sink to the center of the cluster and participate in dynamical encounters that will result in interesting signatures such as tidal disruption events and binary mergers.In particular, BBH mergers containing a more massive component are potentially very important LIGO/ LISA sources (Rodriguez et al. 2019;Kremer et al. 2020a;González et al. 2021;Weatherford et al. 2021;González Prieto et al. 2022).
The clusters modeled in this work could be similar to the massive star clusters that are believed to be proto-GCs.These clusters are thought to be massive and dense, with low metallicity.Our models can thus help constrain the initial properties of clusters that might be the birthplace of the seeds for the very high-redshift quasars observed by many telescopes, including JWST.

Caveats and Future Work
The study of the evolution of single massive stars is currently an active area of research.As an added layer of complexity, most massive stars are observed in close binary systems, making their modeling even more challenging (e.g., Sana et al. 2012;Moe & Di Stefano 2017).As a consequence, we must make assumptions when modeling the evolution of these massive stars.A crucial parameter is their stellar radii, which is poorly constrained for massive stars.As outlined in Section 2.2, this has been partially addressed in this study by the rescaling of the stellar radii, but more accurate modeling is needed.
Even more uncertain are the properties of the collision products.This aspect is particularly relevant to this study since we investigate the formation of massive stars resulting from numerous stellar collisions.Here, we adopt the simple "sticky sphere" prescription for stellar collisions, assuming there is no significant mass loss.Freitag et al. (2006b) showed that this assumption is a good approximation in old clusters with lowvelocity dispersion of the type considered here.Nevertheless, because this prescription is the most "optimistic" scenario, the results of this study are upper limits in the formation of VMSs in star clusters.
A major source of uncertainty concerns the interior structure and radius of the collision product, particularly the effective size of the product.In cases where the VMS experiences exponential growth, the timescale between collisions is shorter than the Kelvin-Helmholtz timescale.This implies that the collision product does not have time to relax back into equilibrium before the next interaction.This is known as the "transparency problem" (Lightman & Shapiro 1978).The hydrodynamics of an interaction that involves this kind of object and another star are not well understood and need to be explored in future studies.
There have been remarkable strides in the field of modeling stellar collision products.In particular, recent studies led by Ballone et al. (2023) simulated the collision of two massive stars using the smooth particle hydrodynamics (SPH) code Starsmasher.The study found that the resulting stellar remnant only experienced 12% mass loss during the merger.Costa et al. (2022) modeled the evolution of this collision product using PARSEC and MESA, concluding that a BH in the upper mass gap was formed as a product of the collision.These studies represent some of the first steps toward carefully modeling collision products.They also shine light on the fact that this process is very intricate and the outcome depends on the properties of the colliding objects.Thus, detailed modeling of the hydrodynamics of stellar encounters and the properties of stellar merger products is essential to better understand VMS growth and IMBH formation from VMSs.
In future studies, we plan to research the prolonged impact and eventual collapse of these VMSs in star clusters.By running a subset of the models to a Hubble time, we aim to study cluster morphology, hypervelocity stars, and tidal tails.

Figure 4 .
Figure 4.Upper panel: the time evolution of the Lagrange radii-from top to bottom enclosing the 99%, 50%, 10%, and 1% of the cluster's total mass-for simulations a1, a4, and a7 in Table1.These models all have N = 4 × 10 5 , r v = 0.5 pc, and f b,high = 0.25.We vary α 3 = (1.6,2.3, 3.0) shown in black, purple, and green curves, respectively.The dashed vertical line represents the main-sequence lifetime for the most massive stars in the cluster (t = 3 Myr).The dotted vertical lines show the time of the first supernova in each cluster model.Lower panel: the collisional history of the three most massive bodies in the simulation with a top-heavy IMF (α 3 = 1.6), distinguished by color.Two of these massive stars (blue and purple) merge at t ≈ 2.3 Myr, and the remnant (colored blue) merges with the third massive star (yellow) at t ≈ 4.4 Myr to form the final VMS.

Figure 3 .
Figure 3.The cumulative distribution of the total number of collisions that contributed to the formation of the most massive star in each of the models.Each panel shows the distribution across all model realizations, differentiated by virial radius (left panel), high-mass stellar IMF slope (center), and high-mass binary fraction (right).

Figure 5 .
Figure 5.The maximum stellar mass as a function of initial cluster parameters, as given by the fitting formula Equation (3).The shaded region represents the 95% credibility region.The mean of the CMC data is overplotted, with the error bars indicating the maximum and median values obtained in our models.Additional models that were not utilized in the parameter estimation process but are used to demonstrate the performance of the predictive model are shown in triangles.Supplementary versions of this figure with the axes instead showing α 3 and f b,high are available as a figure set.(The complete figure set of 3 images is available.)

Figure 7 .
Figure 7. Distribution of the maximum VMS mass across all of the 50 realizations of the model with N = 1.6 × 10 6 , r v = 1 pc, α 3 = 1.6, and f b,high = 1.0.This captures the stochasticity in the outcomes of the collisional runaway.The values from the first three runs are shown in orange.A Gaussian with the same mean and standard deviation is shown in black.

Figure 8 .
Figure 8.The cumulative distribution of the total number of collisions that contributed to the formation of the most massive VMS in each of the 50 realizations of the model is shown in blue.The orange line specifically denotes the total number of collisions where the colliding star had a mass 15 M e .The total number of interactions the star experienced is shown in black.

Table 1
List of Cluster Models