Supermassive black holes from runaway mergers and accretion in nuclear star clusters

Rapid formation of supermassive black holes occurs in dense nuclear star clusters that are initially gas-dominated. Stellar-mass black hole remnants of the most massive cluster sink into the core, where a massive runaway black hole forms as a consequence of combined effects of repeated mergers and Eddington-limited gas accretion. The associated gravitational wave signals of high-redshift extreme mass-ratio inspirals are a unique signature of the nuclear star cluster scenario.

Introduction.The emergence of supermassive black holes (SMBHs) in the early Universe is an outstanding problem in modern cosmology and astrophysics.Many active galactic nuclei (AGN) powered by SMBHs have been identified at redshift  ≳ 6, with spectroscopically inferred black hole (BH) masses in the range 10 6 -10 8  ⊙ (Greene et al. 2024;Juodžbalis et al. 2024;Maiolino et al. 2024).As suggested by recent spectroscopic surveys with the James Webb Space Telescope (JWST), the abundance of red-selected AGN is a hundred times that expected in this redshift range (Matthee et al. 2024;Maiolino et al. 2023).The consensus is that SMBHs have grown out of massive seed black holes (Volonteri 2012), and the problem is then reduced to explaining the formation of those intermediate-mass BH seeds (Greene et al. 2020).One possible formation channel is the direct collapse of pristine gas clouds at  > 10 ( Begelman et al. 2006;Latif et al. 2022).Alternative seeding mechanisms involve repeated BH mergers in dense star clusters (Quinlan & Shapiro 1987, 1989;Davies et al. 2011;Lupi et al. 2014;Antonini et al. 2019;Kroupa et al. 2020;Fragione & Silk 2020;Atallah et al. 2023;Kritos et al. 2023) and feeding by debris of successive tidal disruption events (TDEs; Stone et al. 2017;Inayoshi et al. 2020;Rizzuto et al. 2023).Moreover, several works have studied the rapid formation of supermassive stellar objects with masses in the range 10 3 -10 5  ⊙ within metal-poor massive star clusters (Boekholt et al. 2018;Tagawa et al. 2020;Das et al. 2021;Schleicher et al. 2022;Reinoso et al. 2023) and nuclear star clusters (Escala 2021;Vergara et al. 2023) through runaway collisions of stars and gas accretion.These supermassive stellar objects can then collapse into intermediate-mass BH seeds.
Nuclear star clusters (NSCs) are extremely dense and massive stellar systems that have been observed to reside in the centers of most galaxies in the local Universe (Neumayer et al. 2020).In the nearby Universe, NSCs span the mass range from 10 5  ⊙ to 10 9  ⊙ , with evidence for in situ star formation in the more massive systems (Fahrion et al. 2022).At high redshift, a number of lensed star ★ E-mail: kkritos1@jhu.educlusters with effective radii as small as a parsec (pc) and with masses of order 10 6  ⊙ have been identified (Vanzella et al. 2023;Adamo et al. 2024).There is some theoretical support for the fragmentation of the central region of protogalaxies into millions of metal-poor stars, forming clusters with half-mass radii ≲ 0.5 pc (Devecchi et al. 2010).
In this paper, we use a semianalytic NSC model and demonstrate the possibility of rapidly forming SMBHs through the combined processes of repeated BH mergers and Eddington-limited gas accretion, assuming that a substantial amount of gas can be retained within the cluster for a few tens of millions of years (Myr).Nuclear star cluster model.Our numerical model for the NSC implements a two-mass system consisting of  ★ stars, each of mass  ★ , and  BH stellar BHs, each of mass  BH and nonspinning, with an initial mass  g of residual gas not turned into stars.We find that typically one of the centrally located and initially stellar mass BHs grows dramatically via gas, star, and BH accretion, and we denote its mass and spin by  BH and  BH , respectively.The total cluster mass is therefore  cl =  ★  ★ +  BH  BH +  g +  BH .Since the growing BH is one of the stellar BHs, initially  BH,0 =  BH and  BH,0 = 0.This "oligarchic" growth scenario can occur in dynamical environments (see e.g.Matsubayashi et al. 2004;Kovetz et al. 2018), where the formation of multiple massive BHs within the same system is significantly suppressed by their preferential accretion by the most massive object (Kritos et al. 2023).
The gas, composed primarily of ionized hydrogen with sound speed  s ≃ 10 km s −1 (Krumholz & Matzner 2009), is removed exponentially over a characteristic gas expulsion timescale  ge (Baumgardt & Kroupa 2007).The growing BH may accrete some of this gas.We assume an Eddington-limited Bondi accretion rate onto the BH seed  acc with Salpeter timescale of 50 Myr, and a time-evolving Plummer gas density  g = 3 g /[4( h /1.3) 3 ], where  h is the halfmass radius.As a consequence of gradual gas removal,  h expands adiabatically (Hills 1980).
The half-mass relaxation time  rh , modified for two-mass systems following Eq.( 10) from Antonini & Gieles (2020), is the character-istic time over which a fraction  ≃ 0.1 of the total cluster energy ) can be transferred throughout the cluster via two-body relaxation (Breen & Heggie 2013).The core is dominated by BHs, which partially decouple from stars due to the impossibility of energy equipartition ("Spitzer instability").The core collapses on a time-scale  cc given by Eq. ( 5) from Portegies Zwart & McMillan (2002).
At this point, three-body hard binaries (3bb) can form and heat the system.The rate of energy production  cl = −  cl / rh depends on  h , and is controlled by the efficiency of two-body relaxation ("Hénon's principle").We compute the energy generated by a single 3bb  3bb , and use this condition to calculate the 3bb rate  3bb =  cl / 3bb (as in Heggie & Hut 2003, page 244).If all BHs but one have been ejected from the system, we assume that energy generation proceeds via tidal disruption of stars accreting onto the runaway BH.The TDE rate is  TDE =  cl / TDE , where  TDE =   BH  ★ /(2 T ) is the heat released from a single TDE and  T is the tidal radius (Rees 1988).The main advantage of our method, as in Antonini & Gieles (2020), is that it is not necessary to calculate the properties of the core, while we can incorporate the effect of a cusp around the SMBH.
The differential equations that govern the evolution of the system are given by Equations ( 1a) and (1b) can be compared with Alexander et al. (2014).The main difference is a loss term due to TDEs,  TDE , on the right-hand side of Eq. (1b).The first three terms on the righthand side of Eq. (1f) follow Antonini & Gieles ( 2020), but we also include an adiabatic expansion term due to gas expulsion.The quantity  ejs in Eq. ( 1c) is the number of single BHs ejected per unit binary.We set  = 0.07,  se = 2 Myr (Antonini & Gieles 2020), and  e = 0.0074 (Alexander et al. 2014).
At each iteration, we sample  3bb , the number of 3bbs formed, from a Poisson distribution with parameter  =  3bb , where  is the time-step of the simulation.Every 3bb is formed with components  BH and  BH , with its initial semimajor axis  0 determined by sampling the initial hardness parameter  0 ≥ 1 (Morscher et al. 2015).Binaries harden until they reach their ejection ( ejb ) or gravitationalwave (GW) radius ( GW ), assuming thermal eccentricity.To compute the merger remnant mass  r , spin  r , and relativistic recoil-or GW kick  k -we implement the equations in Sec.V.A-C of Gerosa & Kesden (2016) and references therein, which are based on numerical relativity fits and assume isotropic spin directions.If  GW >  ejb the binary merges within the cluster through the emission of GWs, and if  k <  esc , the merger remnant is retained in the cluster.In this case, we discretely update the runaway BH mass, spin, and number by the substitutions  BH →  r ,  BH →  r , and  BH →  BH −1, while in any case of ejection,  BH →  BH ,  BH → 0, and  BH →  BH −2.Black hole growth.We implement a fourth-order Runge-Kutta scheme to numerically solve the system of equations ( 1) using a time step of  = 0.01 Myr (chosen to be much smaller than  rh ).We set  ★,0 = 10 7 ,  h,0 = 0.1 pc (where the "0" refers to initial values) and  ge = 100 Myr.Using a Kroupa initial mass function in the range [0.08, 150]  ⊙ , the average initial stellar mass is  ★,0 = 0.59 ⊙ , while the fraction of BH progenitors with mass > 20 ⊙ is ≃ 0.18%.We assume  BH = 10 ⊙ , but to quantify uncertainties, we also repeat the analysis with  BH = 50 ⊙ .In Fig. 1, we show the growth of  BH for four selected values of  g,0 = 0, 10 6 , 10 7 , 10 8  ⊙ .
All simulated systems have  esc,0 > 300 km s −1 , allowing repeated mergers to form a runaway BH (Antonini et al. 2019).The BH rapidly grows by hundreds of repeated BH mergers (Fig. 2, left panel) with a merger rate that peaks at ∼ 100 Myr −1 right after core collapse (see inset in the left panel of Fig. 2) and steadily drops as the system expands.
Significant amounts (20%-30%) of gas can be accreted within 200 Myr, rapidly bringing the BH mass into the SMBH regime.The residual gas does not fragment because of heating by UV radiation from cluster stars.Our accretion channel is similar to the direct collapse scenario but occurs on the Salpeter timescale, rather than on the dynamical timescale.The combination of mergers and accretion helps accelerate the growth because Eddington-limited accretion cannot itself grow a stellar-mass BH seed above 1000 ⊙ on a similar timescale (compare the gray lines in Fig. 1, which ignore BH mergers).
Starting from the initial radius  h,0 = 0.1 pc, the cluster inflates by a factor of 10-100 (depending on  g,0 ), so that  h,0 is on the parsec scale after a Hubble time (see Fig. 2, right panel).When  BH = 50 ⊙ , a larger amount of energy is produced, causing  h to expand to above 100 pc (a feature that resembles the bulges of present-day dwarf galaxies).The rapid mass loss in the first few hundred Myr is due to gas expulsion (see the inset of the right panel of Fig. 2).After core collapse, the gas-richest system (dotted red line) ) as a function of detector-frame GW frequency for five extreme mass-ratio inspirals (EMRIs) at redshift  = 10 (colored lines).We use cosmological parameters from Planck 18 (Aghanim et al. 2020).The initial orbital frequency and eccentricity are indicated in the legend; labels on the colored curves refer to the mass of the primary.undergoes a phase of contraction due to balanced evolution that starts at ≃ 22 Myr.Once the gas has been expelled, the energy production term in Eq. (1f) dominates, and the system starts to expand again at  ≃ 200 Myr.
We have demonstrated with an example that a combination of gas accretion and BH mergers is capable of growing stellar-mass BH seeds into the SMBH regime within only 200 Myr.Most notably, if such a cluster forms at  = 15, then (depending on the initial gas content) an SMBH can assemble by  ≃ 10.These SMBH seeds can then further grow by subsequent gas accretion episodes during gas inflows into the NSC following a galaxy merger.This hierarchical scenario can be studied by combining our formalism with galaxy merger trees.
We recall that very massive NSCs are observed locally in many galaxies.Many of them are inferred to have gas because we observe young stars.We find that NSC evolution inevitably results in expansion by an order of magnitude or more (see the right panel of Fig. 2).In contrast, at formation, many NSCs were massive and compact.This supports our contention that NSCs are a promising environment for generating SMBHs via an enhanced merger rate density in the high-redshift Universe.
The ultimate test of our model would be the observation of individual GW signals from SMBH-BH mergers in NSC cores, which may be possible with future GW observatories (Fig. 3).The GW background produced by these sources is a distinctive signature of the NSC scenario that deserves further study.Our predictions motivate the development of the science case for next-generation GW observatories spanning the Hz frequency range between LISA and Pulsar Timing Arrays.Our proposed GW background should be taken into account when estimating astrophysical confusion noise sources that can affect the sensitivity of proposed Hz detectors, such as Ares.For comparison, the black-dotted line in Fig. 3 shows an estimate of the confusion noise from SMBH and white dwarf binaries in the Ares band (Sesana et al. 2021).

Figure 1 .
Figure1.Time evolution of the growing BH mass BH for different initial amounts of gas in the system.The gas expulsion time-scale is set to  ge = 100 Myr, the initial stellar number to  ★,0 = 10 7 , and the initial half-mass radius to  h,0 = 0.1 pc.Thick (thin) lines correspond to simulations with  BH = 10 ⊙ ( BH = 50 ⊙ ).For comparison, we also show how a 10 ⊙ seed (thick gray) and a 50 ⊙ seed (thin gray) would grow under the sole effect of Eddington accretion.

Figure 2 .Figure 3 .
Figure 2. Time evolution of the cumulative number of BH mergers  me (left panel), merger rate  me (left panel, inset), half-mass radius  h (right panel), and cluster mass  cl (right panel, inset) under the same initial conditions used in Fig. 1.