title: The mass budget for intermediate-mass black holes in dense star clusters
authors: Yanlong Shi Michael Y. Grudić Philip F. Hopkins
text: 
  Vehicle Telematics Via Exteroceptive Sensors:
   A Survey
      Fernando Molano Ortiz, 
          Matteo Sammarco, 
          Luís Henrique M. K. Costa, Senior Member, IEEE
          and Marcin Detyniecki
  F. M. Ortiz and L. H. M. K. Costa are with the GTA/PEE-COPPE/DEL-Poli, Federal University of Rio de Janeiro, Brazil, (e-mail: {fmolano, luish}@gta.ufrj.br)
  M. Sammarco and M. Detyniecki are with Axa, France (e-mail: {matteo.sammarco, marcin.detyniecki}@axa.com).
      December 30, 2023
  =================================================================================================================================================================================================================================================================================================================================================================================================================
  
  
  
  
  Intermediate-mass black holes (IMBHs) could form via runaway merging of massive stars in a young massive star cluster (YMC). We combine a suite of numerical simulations of YMC formation with a semi-analytic model for dynamical friction and merging of massive stars and evolution of a central quasi-star, to predict how final quasi-star and relic IMBH masses scale with cluster properties (and compare with observations). The simulations argue that inner YMC density profiles at formation are steep (approaching isothermal), producing some efficient merging even in clusters with relatively low effective densities, unlike models which assume flat central profiles resembling those of globular clusters (GCs) after central relaxation. Our results can be approximated by simple analytic scalings, with M_ IMBH∝ v_ cl^3/2 where v_ cl^2 = G M_ cl/r_ h is the circular velocity in terms of initial cluster mass M_ cl and half-mass radius r_ h. While this suggests IMBH formation is possible even in typical clusters, we show that predicted IMBH masses for these systems are small, ∼ 100-1000 M_⊙ or ∼ 0.0003 M_ cl, below even the most conservative observational upper limits in all known cases. The IMBH mass could reach ≳ 10^4 M_⊙ in the centers nuclear star clusters, ultra-compact dwarfs, or compact ellipticals, but in all these cases the prediction remains far below the present observed supermassive BH masses in these systems.
  
  
  
  stars: black holes – quasars: supermassive black holes – globular clusters: general – galaxies: star clusters: general – galaxies: formation – stars: formation
  
  
  
  
  
  § INTRODUCTION
  
  
  
  Intermediate massive black holes (IMBHs), which typically weigh 10^2–10^5 M_⊙, are believed to be the missing link between stellar mass black holes and super massive black holes (SMBHs). These objects, if they exist, are expected to play an important role in multiple astrophysical processes, e.g., affecting the evolution of globular star clusters and powering off-nuclear ultraluminous X-ray sources (ULXs). More importantly, they are potentially the progenitors of SMBHs which are known to live in most galaxies <cit.>. Observations of ULXs and stellar kinematics argued that there may be some evidence for such objects in galaxies <cit.>, and globular star clusters <cit.>. But these claims remain controversial.
    
  Theoretically, several different IMBH formation channels have been proposed. Major ideas include: direct collapse of hyper-mass quasi-stars in isolation <cit.>, runaway hyper-Eddington accretion onto stellar mass black holes <cit.>, and runaway mergers in globular (star) clusters <cit.>.  All these mechanisms have challenges. For the direct collapse channel, the fragmentation of molecular clouds may not form quasi-stars instantly, but star clusters <cit.>. For the hyper-Eddington channel, one important question is whether such high-efficiency accretion is sustainable or even possible. Finally, for the runaway merger channel, gravitational recoil due to merging stellar-mass BHs will likely “kick” the IMBH to a high velocity <cit.>, sufficient to  make it escape the star cluster. Only relatively massive IMBHs (≳ 10^3 M_⊙) could remain in the galactic field even after the star cluster dissolves and survive such a merger without too-large a “kick” <cit.>, which means the IMBH must be that massive before stellar-mass BH mergers occur. To solve this problem, the runaway merging process must be rapid enough such that massive stars merge together before they evolve off the main sequence and become black holes individually.
    
  
    
  Previous studies of runaway mergers in star clusters have provided us with a possible scenario: relatively massive main-sequence stars sink to the cluster's center due to dynamic friction from the background, the stars then merge into a supermassive quasi-star which then self-collapses to an IMBH after ∼3 Myr. In <cit.>, the authors showed N-body simulations of the process, and found that star clusters with initial half-mass relaxation time scale t_ rlx≲ 25  Myr can form IMBHs. 
  
  More precise simulations in <cit.> drew a similar conclusion and predicted that the quasi-star's mass could account for ∼0.1% of the total cluster mass. The studies focused on cluster dynamics, while the evolution of the quasi-star is another important key step of the runaway merger scenario. Studies have found that quasi-stars' mass (M_ q hereafter) can reach up to 10^6 M_⊙ in principle (given infinite “fuel”), and the remnant BH mass is M_ BH∼ 0.1M_ q <cit.>. These models also find the quasi-stars' lifetime to be ∼ 3  Myr, with only a weak dependence on their masses.
    
  These models, taken at face value, however, would actually imply almost no IMBHs in GCs or other dense stellar systems. The problem is that most GCs (let alone nuclear stellar clusters or galaxy bulges) have half-mass relaxation time scale much longer than 100 Myr (e.g., t_ rlx∼ 2.5 Gyr for M15). However, the studies cited above assumed the initial mass profile of GCs was essentially identical to the mass profiles of nearby relaxed clusters observed today (e.g., with a flat King-type central density profile). In short, if one were to assume that the GCs' present-day mass distribution reflects their mass distribution at formation, this would rule out the runaway merger channel in most globular clusters. However, calculations following the dynamical evolution of globular clusters over cosmological timescales unanimously find that this is not a good assumption <cit.>. Rather, the combined effects of stellar evolution and mass loss, dynamical ejections, mass segregation and “binary burning", and tidal heating/stripping all tend to puff up and flatten the central mass profile slope of dense stellar systems (usually on time scales far shorter than the N-body relaxation time), implying that many presently-observed clusters once had much denser inner cores.
    
  Indeed, the closest observable cousins to proto-globular clusters, young massive clusters (YMCs), are generally found to have density profiles that are significantly different from old globular clusters of comparable mass. Their half-mass radii are generally smaller, with a typical half-mass radius of ∼ 1  pc that has no clear correlation with mass <cit.>, and thus their relaxation times are generally shorter. They also have a relatively compact density profile with an outer asymptotic power-law slope ρ∝ r^-η_1, where η_1 is typically in the range 2-3 <cit.>. <cit.> further found that hydrodynamical simulations of YMC formation were able to reproduce this density profile robustly, and proposed that these density profiles arise from the star cluster assembly process.
    
  In this article we revisit the basic physical processes involved in the assembly a massive stellar object in the centre of a dense star cluster, using the results of the <cit.> simulations, which successfully reproduce observed YMC outer density profiles <cit.>, as well as a range of giant molecular cloud (GMC) properties including their turbulent structure, magnetic field strengths, and stellar auto-correlation functions or stellar clustering <cit.>. These simulations attempt to capture (to the extent possible with state-of-the-art simulations) the cluster properties as they form, which is the most relevant time for potential IMBH formation. Using these simulation results to guide our space of cluster models, we perform a set of Monte Carlo (MC) simulations to track the mass segregation process and study the evolution of the central mass (the quasi-star). Using these methods we predict the mass growth history of the quasi-star and its dependence on the properties of the progenitor cloud or host cluster.
    
    
  The article is organized as follows: in Section <ref> we introduce the analytical and numerical methods used to study mass segregation and the runaway growth of massive objects in star clusters; in Section <ref>, we show the numerical results from the MC simulations and discuss some secondary effects; in Section <ref>, we expand the discussion to observational aspects and make predictions; finally, in Section <ref> we summarize our main findings. 
  
  
  
  
  § MODELS AND METHODS
  
  
  
  
  
  
  
  
   §.§ Initial Conditions from Cluster Formation Simulations
  
  
  For our initial conditions, we extract catalogues of star clusters as they form in the simulations from <cit.>. These are N-body plus magneto-hydrodynamic (MHD) simulations of cloud collapse and star formation, including detailed models for radiative cooling and chemistry, star formation, and “feedback” once stars form in the form of radiation (e.g. radiation pressure and HII regions), stellar winds, and supernovae. The simulations follow the collapse of giant molecular clouds, the assembly of star clusters, and the eventual dispersal of gas due to stellar feedback. One such simulation of e.g. a massive complex can produce many independent clusters: we identify gravitationally-bound star clusters remaining after gas dispersal[We have also compared the results extracting clusters at the time of peak star formation; the time difference is sufficiently small that it has little effect on our results.] using a group-finder which associates stars belonging to a common potential well which are also gravitationally bound within that well (see Appendix A in <cit.> for details). We restrict to clusters which form >100 bound star particles.
  
  This gives us an ensemble of ∼1000 clusters “at formation,” one of which is shown in Fig. <ref> Note that this sample of clusters should not be considered statistically representative of a cluster population that would form in a real galaxy: the initial conditions of the simulations were simply uniformly sampled on a logarithmic grid in mass-radius parameter space, which is ideal for our study here. 
  
  These simulations are designed to (1) sample an enormous parameter space, and (2) simulate large complexes through the entirety of star formation and stellar evolution: as such, the numerical resolution is such that individual “star particles” in the original simulations represent an IMF-averaged ensemble of stars. To properly evolve stellar dynamics, we therefore re-sample each star particle into an ensemble of individual stars, drawing probabilistically from the stellar initial mass function (IMF) conserving total stellar mass. By default (since it is the same used for the original simulation stellar evolution models) we adopt a <cit.> IMF with an upper mass limit of m_ max=100 M_⊙.[We have tested and adopting instead a <cit.> makes negligible difference to our conclusions. Likewise we find that varying the upper “cutoff” mass of the IMF makes only weak (logarithmic) corrections to our predictions (because these stars contribute negligibly to the total massive-star stellar mass budget).] For this assumption the median stellar mass is ⟨ m ⟩≈ 0.38 M_⊙, and the mean is 1.5 M_⊙.
  
  
  
  
   §.§ Analytic Models from the Simulations
  
  
  Although our simulation suite is extensive, it is still limited by (1) finite sampling of parameter space and (2) finite resolution. Especially in cluster centers (particularly important here), the original simulation will always produce finite-resolution effects. Moreover, although the simulated clusters have some non-axisymmetric structure, we generally find this is small and generates torques which are weak compared to dynamical friction (discussed below). Therefore, it is especially useful to also consider general analytic models for the initial conditions, motivated by the cluster catalogue from our simulation suite. 
  
  We consider three simple, spherically-symmetric analytic density profiles, which we will show allows us to capture almost all of the key behaviors we study. These are shown in Fig. <ref> as fits to one example profile.[We have experimented with a variety of different methods for fitting the analytic profiles to the simulation outputs, and find the most robust results fitting directly to the spherically-averaged ρ(r) in log-log space with uniform weights but constraining the analytic fit to reproduce the total mass and half-mass radius (specifying ρ_c and r_ c) exactly, so only the slopes η_1,2 are “free.”] First, a “Modified Jaffe”  model (from , Eq. 2.64): 
  
      ρ(r) = ρ_ c(r/r_ c)^-η_1(1+r/r_ c)^-η_2+η_1.
  
  with inner power-law slope η_1, outer slope η_2, turnover radius r_ c, and normalization ρ_c (given by e.g. the total mass). We also consider a similar “double power” law model: 
  
      ρ(r) =   ρ_ c/(r/r_ c)^η_1+(r/r_ c)^η_2,
  
  which has the same qualitative features as the “Modified Jaffe” model but features a much sharper turnover around r_ c, which is useful in what follows as it dramatically reduces the covariance between the parameters η_1 and η_2. 
  
  Finally, we also consider a “single power” law model: ρ(r)=ρ_ c(r/r_ c)^-η for r<r_ c. This obviously cannot fit any mass profile over the entire dynamic range of r with finite mass; we therefore restrict the fit only to radii smaller than the half-mass radius (so η≈η_1). This is included here because it allows us to derive some simple analytic expressions in regimes where the inner profile dominates the behavior. 
  
  
  
  
  
  
  
   §.§ Sinking and “Merging” Stars
  
  
  
  
  Even with the simplifications above, integrating the full N-body dynamics of massive stars through a cluster into merger with a central object is computationally impossible (both given our large parameter space of models and sample of “clusters” reaching ∼ 10^8 M_⊙, let alone uncertainties in the actual size/evolution of the central object). However, full N-body studies of a small number of smaller clusters <cit.> have shown that dynamical friction is an excellent approximation to rate of sinking and merger of massive stars with m≫⟨ m⟩ (which are those that dominate the buildup of a central quasi-star on the timescales of interest). This quickly circularizes the orbits of the massive stars and leads to orbital decay with 
  
      ṙ =    -4π G^2 m lnΛ v_c^-3 ρ_b r
  
  where m is the mass of the sinking star, v_c^2≡ G M(<r)/r reflects the enclosed mass M(<r) inside r, ρ_b is the density of the background stars at radius r (e.g. the ρ(r) in the profiles above), and Λ is a Coulomb logarithm which we take to be Λ(r) ≈ 0.1 M(<r)/⟨ m ⟩. Less massive stars will not sink: we approximate this (conservatively, for now) by simply applying a cutoff ignoring any ṙ below m_ min=8 M_⊙.[This ignores back-reaction causing lower-mass stars to migrate outwards, but this is a small effect on the timescales we consider, and we show below the exact choice of m_ min also has relatively weak effects on our conclusions.] 
  
  From a MC realization of the ICs (either directly from the formation simulations, or analytic fits), we then evolve the system forward in time. As massive stars approach the center, the first to reach the center (region interior to which there are no other m>m_ min stars) becomes the “seed” quasi-star (with mass M_c). Because we are interested in mergers while the massive stars are on the main sequence, we subsequently “merge” into this any massive star which (1) has not yet reached the end of its main sequence lifetime (adopting the relation from ; typically ∼3 Myr for the most massive stars), (2) approaches the quasi-star within a radius r < r_q which represents some “interacting binary” or “common envelope” radius (for which we take the value quoted by <cit.> for models of a rapidly-accreting protostar: r_q≈ 2600 R_⊙ (M_q/100 M_⊙)^1/2),[This is essentially an extrapolation from “normal” pre-main sequence stars. Of course the sizes of quasi-stars are purely theoretical and uncertain: however varying this by factors of several has very little effect on our conclusions, as the “sinking” times around these radii are relatively small. But we need to include some finite “merger radius” since we do not model effects like gravitational wave emission which could merge point-mass-like particles.] and (3) reaches before the quasi-star itself has reached the end of its lifetime. We simply add the merged mass to M_c, neglecting e.g. mass-loss associated with the merger.
  Note that during collisions the central object mass M_c(t) will contribute to the total enclosed mass M(<r) as an additional point mass, which is included as a correction when solving Eq. (<ref>) in our full version of MC simulations. This “point mass correction” (PMC) is not included in our analytical calculation and the corresponding MC simulations.
  
  The quasi-star “lifetime” essentially sets the end of our simulation, and the final mass of the quasi-star. Studies of quasi-star structure <cit.> have found that because these stars are approximately Eddington-limited, they have lifetimes ∼3 Myr akin to massive stars. We have therefore considered simply taking the mass M_q = M_c(t=3  Myr). We have also considered a more sophisticated model motivated by the same pre-main-sequence models described above: some accretion rate Ṁ_ acc = d M_c/d t from mergers sustains the quasi-star lifetime and keeps it “puffed up” (allowing efficient mergers) as long as it is larger than the fuel consumption/loss rate from a combination of nuclear burning and stellar mass-loss, which occurs on a characteristic timescale t_0∼ 3  Myr. We therefore take the final M_q to be M_c at the first time where Ṁ_ acc falls below M_c/t_0. In practice, because the merging stars also have lifetimes ∼ 3Myr, it makes very little difference which of these assumptions we adopt.
  
  
  
  
  
  
  
  § RESULTS
  
  
  
  
  
   §.§ Density Profiles
  
  
  
  <cit.> showed that the simulated clusters here produce a distribution of density profile shapes after relaxation in good agreement with observations; however, no analysis of the inner density profiles at formation was performed. In Fig. <ref>, we show the distribution of the inner (η_1) and outer (η_2) mass profile slopes fit to all clusters. In both “modified Jaffe” and “double-power-law” models, the outer slopes are typically in the range of 3.5-4.5 as found in <cit.>. The inner slopes cluster around 1-2.5 (as compared to post-relaxation profiles, which broadly follow a “flat”  distribution). The best fit distribution for η_1 is most narrowly-peaked (around η_1≈ 2) for the “single-power-law” fits (fit to just r within the half-mass radius), and most broad for the “modified Jaffe” fit. Our extensive experimentation with different fitting methods indicates that this directly traces the covariance between η_1 and η_2. The single-power fit, with only one slope, has no η_1-η_2 covariance. The double-power fit, with a “sharp” break, has weak covariance between η_1 and η_2 which “smears” the best-fit η_1. The modified Jaffe fit exhibits very strong covariance between η_1 and η_2, with a wide range of allowed fits for any given simulation profile.[It is important to note that because of the covariance in the fits, models with the modified Jaffe fit with η_1∼ 0 give r_c≪ r_ half, i.e. the “rollover” occurs very slowly down to extremely small radii (often below the simulation resolution) – so the central densities are still large.]
  
  
  Fig. <ref> shows that the best-fit η_1 and η_2 do not depend systematically on either cluster half-mass radius r_ h or mass M_ cl. Likewise our cluster catalogue includes simulations with progenitor clouds of different metallicities (Z/Z_⊙=0.01-1), and we see no dependence on Z. 
  
  Even though the exact values of the “inner slope” η_1 can vary between fits, the fact that these are relatively narrowly constrained (within covariances, i.e. all are “good” representations of the data), and that they do not depend systematically on cloud mass/size, means that we obtain reasonably good estimates for the final central object mass from our full MC calculation, using an idealized mass profile fit with either of the fitting functions or assuming a “universal” mass profile shape across all clusters (Fig. <ref>). 
  The modal fit values for (η_1, η_2) are (η_ 1M, η_ 2M) = (1.2, 3.7) for the modified Jaffe, (1.4, 3.5) for the double-power-law, and (1.9) for the single power-law model. Because of the covariances, however, these are not the same as the values which give the best estimate of M_q compared to our full MC calculation. Instead, we should ask which values (η_1∗, η_2∗), applied to the entire ensemble of clouds (as a “universal” profile shape), most accurately predict M_q from the MC:[Formally we find the (η_1∗, η_2∗) which minimize the variance ∑ |log{M_q,  pred(η_1∗, η_2∗, r_ h, M_ cl)} - log{M_q( MC)}|^2.] these are (η_1∗, η_2∗)=(1.68, 4.95) for modified Jaffe,[The dramatic change in the values for modified Jaffe again indicates the covariance (with the η_2=4.95 value indicating that the outer slope plays a very small role in determining M_q).] (1.81, 3.79) for double power-law, and (1.93) for single power-law, all of which feature a similar, isothermal-like inner slope η_1 ∼ 2.  
  
  
  
  
  
  
  
   §.§ Central Object Growth
  
  
  
  Fig. <ref> shows one example of our full MC simulation, with the ensuing growth of the central object as a function of time M_c(t). The results from the spherical analytic model with the modified Jaffe or double-power-law profiles agree very will with the full MC (the single-power-law works well up to ∼ 1-3Myr, as well, where most of the accretion is from radii ≪ r_ h). The “crossover point” where Ṁ_ acc = M_c/t_0 occurs at t ∼ 1.3Myr, but the growth rate of M_c is slowing down already at this point, so M_q only differs by a factor of ∼ 1.5 if we take M_q to be M_c at t=3Myr, or a factor of 1.9 if we take M_q as t→∞. In any case, this particular cluster, chosen to be relatively extreme (with a total mass ∼ 10^6  M_⊙ and initial central density of ∼ 10^11  M_⊙ pc^-3 at r≲ 0.001  pc) is able to merge most of its massive stars (∼ 10% of the total stellar mass) within <3Myr. 
  
  Fig. <ref> shows M_q from our full MC calculation for each simulated cluster, as a function of cluster mass and half-mass radius. There is a clear trend where more massive clusters M_ cl at the same size give larger M_ q, and a weaker but still evident trend of larger M_ q for more compact clusters at fixed mass. These are expected if cluster profiles are approximately self-similar: to show this we compare the predicted M_ q from analytic models with different M_ cl and r_ h, assuming a universal density profile shape (the double power-law fit with fixed η_1 = η_ 1∗,  η_2 = η_ 2∗).
  
  Fig. <ref> compares the mass fraction which can merge to the center f_q≡ M_q/M_ cl from our full MC calculation to that obtained from the simple spherical analytic models. We compare the results from the modified Jaffe and double-power-law fits, fit individually to each simulation cluster, which predict f_q to within <10% on average – this indicates that deviations from symmetry, “lumpiness” or irregular structure in the potential and density profile, or resolution effects (e.g. numerical flattening or shot noise in the central density profile, as compared to the profile generated by a smooth power-law down to r→ 0) do not strongly influence our results. We also show the results assuming a universal profile with (η_1, η_2) = (η_1∗, η_2∗). This increases the scatter (as expected) but only by a modest amount: we can predict M_q to within an rms <0.15 dex assuming this universal shape at formation.
  
  
  It is unclear exactly at which mass scale dynamical friction ceases to be a good approximation for the “sinking” of massive stars: Fig. <ref> varies the minimum mass m_ min we allow to sink, to show this has only a small effect on our predictions. Varying m_ min from 0.5-8  M_⊙ changes M_q by a factor ∼ 2, because (a) lower-mass stars sink more slowly (even if we allow them to sink), and (b) the Salpeter IMF is not extremely steep, so the total mass of stars “sinking” only changes with m_ min^-0.3. 
  
  Because most of the mass in the IMF is not in the highest-mass stars, it also makes little difference if we vary the high-mass cutoff (e.g. changing the upper-mass cutoff of the IMF from 100 M_⊙ to 200 M_⊙ only produces a ≲10% difference in M_ q). 
  
  Another significant uncertainty in our models is how the actual mergers/coalescence occur in the center: we simply populate stars and merge anything within some large radius in the center (reflecting the envelope size of the quasi-star), implicitly meaning there is some “overlap” between the envelope of the quasi-stars and our populated stars in the models. Properly determining if or how mergers once massive stars sink close to the central quasi-star requires dynamical stellar merger simulations. But even within our simple model, stars can still “stall” near the center. In Fig. <ref>, we consider a model variation where we simply merge any star which reaches the radius where the proto-star would dominate (be more than 1/2 of) the enclosed mass M(r), and ignore the mass of M_q itself in calculating v_c in Eq. <ref>: these changes essentially guarantee that any massive star which approaches small r merges. We see that this systematically increases M_q/M_ cl, as expected, by a factor ∼ 2. This in turn means that our “default” model is predicting an order-unity fraction of massive stars near r∼ 0 “stall” or otherwise fail to merge, a reasonable order-of-magnitude approximation to few-body simulations. It also implies that these correction does not change our qualitative conclusions.
  
  
  
  
  
  
  
  
  § COMPARISON TO PREVIOUS WORK
  
  
  As discussed in  <ref>, <cit.> and <cit.> considered detailed N-body simulations to follow runaway merging, but used adopted very different mass profile shapes (similar to post-relaxation clusters today). If we adopt a similar profile shape to their default (e.g. a Plummer-like (η_1, η_2) = (0 , 5)), and then run our full model to calculate M_q for a wide range of M_ cl and r_ h (sampling the values of our simulation library) then – despite the other  simplifications here – we obtain quite good agreement (within a factor of a few) with both their requirement that the cluster must have t_ rlx≲ 25Myr to produce any appreciable growth of M_ c and the peak mass M_ q or resulting mass fraction f_q=M_q/M_ cl in a central massive object produced when this criterion is met. This is reassuring, and implies our methodology is reasonable. The key difference in our predictions, compared to theirs, arises because our MHD star formation simulations predict quite different values of η_1 compared to those they considered.
  
  Some other recent studies have considered runway stellar mergers in initial conditions closer to those here (but with a more limited or ad-hoc choice of initial profiles). <cit.> adopted a similar approach to that here, using hydrodynamic simulations to select dense (M_ cl∼ 10^5 M_⊙, King-profile core r_ c∼ 0.4pc) proto-galactic “clouds” and then using those to set up initial conditions for N-body simulations: although our survey is intended to match much later-forming star clusters, where the cloud properties overlap we find similar M_q within a factor ∼ 2 for each of the ∼ 8 clusters they simulate (assuming a typical r_ c/r_ h∼ 0.1). 
  
  In parallel, <cit.> considered analytic models for cloud/cluster formation, to estimate typical cloud densities in the early universe, coupled to a simple prescription from <cit.> for the fraction of clusters undergoing runaway, to argue f_q could reach ∼ 0.05 for dense clusters formed in the early Universe because these produce steep central profiles (η_1≳ 5/3), broadly similar to our conclusions. And recently, <cit.> performed semi-analytic calculations qualitatively akin to those here, considering a much more limited range of profiles (but taking steep η_1) but much more detailed models for the (proto)-stellar evolution of the quasi-star and merger criteria, but conclude that effective growth ceases at ∼ 3Myr (as we assume here) with a similar effective radius for merger (versus M_q) as we adopt here. More recently, <cit.> performed a series of N-body simulations of YMCs based on King-type profiles, indicating that massive stars weighing up to ∼ 400 M_⊙ may form within 5–15  Myr and sequentially become IMBHs. The results generally support our semi-analytical model, which is based on a more realistic parameter space of YMC density profiles and limited lifetime for quasi-stars.
  
  
  
  
  § DISCUSSION
  
  
  
  We now consider the implications of our results for real dense stellar systems. Fig. <ref> plots the distribution in M_ cl and r_ h of a wide variety of known dense, stellar-dominated, dispersion-supported systems: globulars and YMCs, super-star clusters (SSCs), nuclear star clusters in different dwarf and late-type galaxies (NSCs), ultra-compact dwarf galaxies (UCDs), nearby and high-redshift compact elliptical galaxies and bulge-dominated galaxies (Es). The sizes and masses are compiled in <cit.>, from observations by <cit.>. We have no way of knowing their properties “at formation,” but because our full simulations can be reasonably approximated by assuming a “universal” profile shape at formation, we compare the contours of M_q predicted by our model with a universal at-formation (η_1∗, η_2∗). This does assume that the total mass and size have not evolved much since formation, an obviously uncertain assumption, but likely plausible since most of these systems have N-body relaxation times longer than the Hubble time.
  
  From the detailed studies of quasi-star evolution noted in  <ref>, we will also assume in what follows that any quasi-star leaves behind a “relic” BH of mass M_ BH∼ 0.1 ϵ_0.1 M_q. This “fudge factor” ϵ_0.1 accounts for processes including inefficiency of final mergers, mass loss/ejection during merges, stellar winds and mass loss during collapse from the quasi-star. 
  
  
  
  
  
  
  
   §.§ Analytic Scalings
  
  
  Assuming the “universal” profile parameters, the contours of constant M_q are approximately power-laws over most of the dynamic range of interest. We can approximate this quite well via a simple purely-analytic estimate for M_q if we assume a single-power law profile with η_1=2 (isothermal), neglect “edge” effects (assume the stars outside r_ h are not sinking efficiently), and approximate the effects of various non-linear terms like varying coulomb logarithms, the finite N-body point mass correction, finite IMF sampling, finite size of the quasi-star, and others as a systematic factor of ∼ 2 normalization correction (reasonably well-motivated by our comparison in Fig. <ref>). With all of these approximations, we obtain the very simple expression:
  
      M_ BH   = 0.1 ϵ_0.1 M_q
         ∼ 250 M_⊙ ϵ_0.1( M_ cl, 5/r_ h, pc)^3/4∼ 250 M_⊙ ϵ_0.1 ( V_ eff/20  km s^-1)^3/2
  
  where V_ eff^2≡ G M_ cl/r_ h, M_ cl, 5=M_ cl/10^5 M_⊙, r_ h, pc=r_ h/ pc.
  
  Despite the many simplifications involved in deriving this expression, it provides a quite reasonable order-of-magnitude approximation to the most important results from our more detailed full model calculations. 
  
  
  
  
   §.§ Globulars and “Typical” Dense Star Clusters
  
  
  Fig. <ref> & Eq. <ref> do suggest IMBHs could form in massive GCs, with typical masses M_ BH∼ 0.0003 ϵ_0.1 M_ cl. In Fig. <ref> we also show the criterion t_ rlx < 100Myr which <cit.> and <cit.> argue is required for a GC with an initially flat (η_1 = 0), Plummer-like density profile to undergo any significant runaway merging. We see, as noted in  <ref>, that almost no known present-day massive clusters meet this criterion. The reason our modeling here predicts they can form central objects is because we argue they likely had steeper slopes at initial formation (allowing some merging near their center at these early times). But, essentially by definition, any interior region which has a steep enough slope to produce runaway merging within <3Myr will, on timescales ∼ Gyr, have undergone relaxation, flattening the central profile seen today. 
  
  However, although our models predict runaway merging could occur in the centers of almost all clusters at formation, the actual mass which we predict successfully merges (for realistic cluster M_ cl and r_ h) is quite modest, giving a rather low M_ BH/M_ cl compared to the clusters which undergo “complete runaway core collapse” (with t_ rlx < 25Myr) in <cit.>. In Table <ref>, we explicitly list a number of individual observed GCs from Fig. <ref> which have claimed detections or upper limits for central IMBHs. Many contradictory observational claims exist for some clusters, a well-known issue in the literature. Using the same models from Fig. <ref>, and the observed cluster properties, we give our best estimate of M_ BH (approximately given by Eq. <ref>), and compare to these observations. We see that the predicted relic mass from our calculations is typically ∼ 10^-4-10^-3 M_ cl, in many cases a factor ∼ 10 or more below the claimed detections/upper limits. There is no case where our predicted M_ BH exceeds even the most stringent upper limits. The most constraining examples we find are M15 and NGC 6388: here our “default” prediction is only a factor ∼ 2.5 below the current upper limits or smallest values among the claimed detections. This implies ϵ_0.1≲ 3 (i.e. M_ BH≲ 0.3 M_q, if our models are to be believed at this level of accuracy). 
  
  If some of the most-massive detections (with claimed IMBH masses up to ∼ 70 times larger than our prediction) are indeed correct (although almost all of these cases are controversial with much lower limits claimed by other studies), then it would most likely imply that the central BHs grew rapidly after formation via some other process such as gas accretion (from e.g. stellar mass-loss in the cluster).
  
  
  
   §.§ Connection to SMBHs and More Massive Stellar Systems
  
  
  The predicted quasi-stars/IMBHs suggested in Fig. <ref> & Eq. <ref> become more massive, on average, in more massive systems, reaching M_ q∼ 3×10^4 M_⊙ in the most massive and dense NSCs and UCDs, and up to M_q∼ 3×10^5 M_⊙ in the centers of the most compact local and high-redshift bulges/Es. 
  
  These are systems which are known to host super-massive BHs, with M_ BH∼ 10^4-10^10 M_⊙ unambiguously detected (with the smallest BHs in dwarf NSCs, the most massive in compact Es), obeying a tight correlation with the velocity dispersion σ of the surrounding stars M_ BH^M-σ∼ 3×10^8 M_⊙ (σ / 200  km s^-1)^4.3 <cit.>. However, noting that σ≈ V_ eff in Eq. <ref> for an isothermal profile, this implies that the present-day SMBHs observed are much more massive than the IMBH we predict from runaway merging at formation, for any σ≳ 7 ϵ_0.1^0.3  km s^-1 (present-day M_ BH^M-σ≳ 100 ϵ_0.1^1.5 M_⊙). In other words, while the masses here are potentially interesting for very first initial seeds of the SMBHs, runaway merging cannot establish most of the mass of any observed BHs on the BH-host galaxy (or BH-NSC) scaling relations. It does not substantially reduce the amount of BH accretion (nor the time required for that accretion, if it occurs at e.g. a fixed Eddington ratio), nor even radically change the initial seed mass relative to commonly-assumed M_ BH^ seed∼ 100 M_⊙ as the “most optimal normal stellar relic” remnant mass. 
  
  An important additional caveat in these massive systems is that our models assume the stars are approximately co-eval. This is reasonable in the centers of dense GCs where the dynamical times are ≪ Myr. However, in e.g. elliptical galaxy centers, the dynamical times can be ≫ 3Myr; since stars cannot form much faster than the dynamical time, it is almost certainly the case that the massive star formation was extended in time relative to the lifetime ∼ 3Myr of the quasi-star. In the center, later-forming stars can still sink, but they will merge with a central IMBH instead of quasi-star, producing a tidal disruption event and building up an accretion disk rather than directly forming a quasi-star. 
  
  
  
  
  § CONCLUSIONS
  
  
  
  Using the outputs of high-resolution numerical hydrodynamic simulations of star cluster/complex formation and destruction which have been shown to reproduce a wide range of GMC and cluster observables, we develop a semi-analytic model for the sinking of massive stars to cluster centers and their merger into a massive quasi-star. We find:
  
  	
    * The mass profile of YMCs “at formation” (centered on local peaks, as there can still be substructure) can be described by a double power-law with steep, near-isothermal inner slopes common (which flatten at later times as the inner regions dynamically relax). This means that some runaway merging can occur early even in clusters with relatively low mean densities and long relaxation times (e.g. M_ cl∼ 10^5-10^6 M_⊙, r_ h∼ 1-10pc, with relaxation times t_ rlx≳ 100Myr). The runaway ceases at ∼ 1-3Myr, regardless of the details of the quasi-star evolution.
  	
    * Over the parameter space of greatest interest (where massive, dense stellar systems are observed), our predictions can be approximated with a simple scaling, with the total mass of massive stars which sink to the center and could potentially merge, M_q, scaling as M_q∝ V_ eff^1.5, where V_ eff is a characteristic circular velocity (Eq. <ref>). 
  	
    * Although some runaway merging is predicted in nearly all clusters (Fig. <ref>), the actual masses of IMBH relics predicted in our model for observed globulars and typical dense star clusters are quite modest, ∼ 100-1000 M_⊙. For relic mass M_ BH≲ 0.3 M_q (expected allowing for mass-loss, imperfect merging, quasi-star evolution, etc.), our predictions are consistent with even the most stringent upper limits (to our knowledge) on central IMBH mass in all clusters for which such constraints exist. The most constraining clusters for our models at present are M15 and NGC6388; the only well-studied cluster where our model predicts M_ BH≳ 1000 M_⊙ is G1. 
  	
    * In more massive systems such as nuclear star clusters, ultra-compact dwarfs, and the centers of compact ellipticals, the central object mass could reach M_q∼ 10^4-10^5 M_⊙, an interesting range for initial seeds of super-massive BHs. However for any system with velocity dispersion ≳ 10  km s^-1, the SMBHs on the various observed SMBH-host scaling relations (e.g. M_ BH-σ) are far more massive than even the most optimistic IMBH masses resulting from runaway merging. Thus runaway merging does not significantly reduce the need for subsequent accretion to super-massive BH masses.
  
  
  Our models are intentionally simplified in order to survey a wide parameter space efficiently and guide intuition and predictions for future models. Using the models here to identify the most interesting parameter space, in future work we hope to consider explicit N-body simulations of the merging process in cluster centers (necessarily limited to a small number of realizations). In the systems where merging occurs most rapidly, it is also possible that mergers occur even as stars are still forming in the cluster, potentially before massive protostellar cores even complete their pre-main sequence evolution. Exploring this will require fully hydrodynamic+N-body simulations of star formation which can resolve the stellar IMF self-consistently and follow mergers as they occur “live.” Considerable uncertainties also still surround the actual dynamics of massive stellar mergers (including complicated effects not followed here, such as the effect of resolved binaries and hierarchical multiples on merger efficiency) and the evolution (especially as it grows via rapid merging) of the (proto) quasi-star. In addition, if such a system forms, a variety of processes may allow for rapid growth even after it collapses to an IMBH, as it could accrete tidally-disrupted lower-mass stars which sink on longer timescales (e.g. m∼ 2-8 M_⊙), or stellar mass-loss products from AGB stars that can remain gravitationally bound in the cluster potential. All of these remain important subjects for future study.
  
  
  Support for YS & PFH was provided by NSF Collaborative Research Grants 1715847 & 1911233, NSF CAREER grant 1455342, NASA grants 80NSSC18K0562, JPL 1589742. Numerical calculations were run on the Caltech compute cluster “Wheeler,” allocations FTA-Hopkins supported by the NSF and TACC, and NASA HEC SMD-16-7592.
  
  
  
  