From Seeds to Supermassive Black Holes: Capture, Growth, Migration, and Pairing in Dense Proto-Bulge Environments

The origins and mergers of supermassive black holes (BHs) remain a mystery. We describe a scenario from a novel multi-physics simulation featuring rapid ($\lesssim 1\,$Myr) hyper-Eddington gas capture by a $\sim 1000\,{\rm M}_{\odot}$ ``seed'' BH up to supermassive ($\gtrsim 10^{6}\,M_{\odot}$) masses, in a massive, dense molecular cloud complex typical of high-redshift starbursts. Due to the high cloud density, stellar feedback is inefficient and most of the gas turns into stars in star clusters which rapidly merge hierarchically, creating deep potential wells. Relatively low-mass BH seeds at random positions can be ``captured'' by merging sub-clusters and migrate to the center in $\sim1$ free-fall time (vastly faster than dynamical friction). This also efficiently produces a paired BH binary with $\sim 0.1$\,pc separation. The centrally-concentrated stellar density profile (akin to a ``proto-bulge'') allows the cluster as a whole to capture and retain gas and build up a large (pc-scale) circum-binary accretion disk with gas coherently funnelled to the central BH (even when the BH radius of influence is small). The disk is ``hyper-magnetized'' and ``flux-frozen'': dominated by a toroidal magnetic field with plasma $\beta \sim 10^{-3}$, with the fields amplified by flux-freezing. This drives hyper-Eddington inflow rates $\gtrsim 1\,\rm M_\odot yr^{-1}$, which also drive the two BHs to nearly-equal masses. The late-stage system appears remarkably similar to recently-observed high-redshift ``little red dots.'' This scenario can provide an explanation for rapid SMBH formation, growth and mergers in high-redshift galaxies.


INTRODUCTION
Observations of high-redshift (z ≳ 7, ∼ 0.7 Gyr since the Big Bang) quasars demonstrate the existence of supermassive black holes (SMBHs) of > 10 9 M ⊙ at these very early stages of the Universe (Fan et al. 2001;Yang et al. 2020;Wang et al. 2021).How supermassive BHs form at all -let alone how they grow so quickly -remains a major unsolved problem (Inayoshi et al. 2020;Volonteri et al. 2021).Evidence from low-redshift studies argues such SMBHs grow from "seed" masses (perhaps even as low as stellar 10 − 100 M ⊙ ) via gas accretion (Yu & Tremaine 2002).But assuming accretion is radiatively inefficient (and spherical), the Eddington limit is commonly cited as setting a classical Corresponding author: Yanlong Shi yanlong.astro@outlook.comupper limit to the BH accretion rate, which is defined as ṀEdd = M BH /(ϵ ref t Sal ), where t Sal = κ es c/(4πG) ≈ 45 Myr is the Salpeter time and ϵ ref is an assumed radiative efficiency (typically 0.1; Inayoshi et al. 2020).Thus, to grow to SMBH masses from "seed" masses ≪ 10 6 M ⊙ (Greene et al. 2020) at high redshifts generally requires either accretion well above this naive Eddington limit (Inayoshi et al. 2016;Shi et al. 2023), or exotic IMBH seed formation scenarios like the direct collapse of un-fragmented giant molecular clouds to single hyper-massive stars (with SMBH masses themselves; Bromm & Loeb 2003), exotic dark matter/new physics processes (Xiao et al. 2021), or runaway mergers of stars (Portegies Zwart et al. 2004;Shi et al. 2021;Rantala et al. 2024).
Theoretical work on small-scale (∼ au or horizonscale) accretion physics has shown that super-Eddington accretion is possible (Begelman 1979;Blandford & Begelman 2004;Inayoshi et al. 2016).This is primar-ily enabled by the "photon trapping" effect, in which photons are trapped in the strong accretion flow (∼ 1000 ṀEdd ) and advected into the BH (radiating inefficiently).This is also demonstrated by numerical simulations (Jiang et al. 2014;Sądowski et al. 2015;Jiang et al. 2019) where a sustainable phase of super-Eddington accretion is observed.But these simulations focus only on the accretion disk scales well interior to radii where the BH dominates the potential: another (perhaps more challenging) requirement for super-Eddington accretion is that the BH can gravitationally capture sufficient gas (from its turbulent, star-forming ISM environment) to sustain such accretion for the time needed to grow.This could be especially challenging in high-redshift galaxies where star formation is clumpy and bursty, potentially scattering BH seeds and inhibiting large-scale, coherent accretion flows (Ma et al. 2021;Byrne et al. 2023).Furthermore, classical accretion disk models (Shakura & Sunyaev 1973) predict such a super-Eddington accretion disk should be violently gravitationally unstable and fragment outside of tens of gravitational radii (GM BH /c 2 ) from the horizon (Goodman 2003).
Some of the problems above are explored in Shi et al. (2023) where we embedded BH seeds into star-forming giant molecular clouds (GMCs).We found that significant BH accretion can "stochastically" occur if a dense clump (pc or sub-pc scale) formed via turbulence and stellar feedback shocks (Klessen 2000;Mac Low & Klessen 2004;McKee & Ostriker 2007) is sufficiently gravitationally bound and close (sub-pc) to the BH, and if the gas is sufficiently dense and cold (Inayoshi et al. 2016).In practice, this "clump accretion" only occurs if both the gas and the BH dynamics are dominated by the GMC's self-gravity, and stellar feedback is globally inefficient (unable to rapidly unbind all the gas in the GMC as soon as the first massive stars form; Shi et al. 2023), which is true for GMCs with high initial surface density (Grudić et al. 2018b(Grudić et al. , 2021a;;Chevance et al. 2023).In Shi et al. (2024) we extended these to include physicallymotivated feedback from BH accretion in the forms of radiation (Jiang et al. 2019), accretion-disk winds and jets (Silk & Rees 1998;Blandford & Begelman 1999), and cosmic rays (Guo & Oh 2008;Guo & Mathews 2012;Zweibel 2017).We showed that for reasonable feedback parameters, BHs can still grow rapidly and their growth is primarily "fueling-limited" until they reach supermassive masses and the "feedback-limited" regime.But this alone does not address many of the key questions above.
Meanwhile, in a series of papers using a different numerical method, physics, and technical approach, Hopkins et al. (2024a,b) followed gas accreting onto a preexisting supermassive BH from cosmological initial con-ditions including star formation, and found a novel type of BH accretion disk: one which is "hyper-magnetized" (ratio of thermal-to-magnetic pressure β ≪ 10 −2 ) and "flux-frozen" (dominated by a primarily toroidal magnetic field amplified from typical ISM fields by fluxfreezing).In these papers and a subsequent analytic study (Hopkins et al. 2024c), they showed that these disks could resolve a number of problems posed by classic accretion disk models (like Shakura & Sunyaev 1973) which assume β ≫ 1, most notably that the disk can remain stable (with Toomre Q ≫ 1) out to ≳ pc scales even with accretion rates thousands of times larger than the nominal Eddington limit.
Here, we present a case study of one of the massive, high-density cloud-complex simulations in Shi et al. (2023Shi et al. ( , 2024)), representative of plausible starburst conditions in high-redshift progenitors of massive galaxies.We show that this demonstrates a number of remarkable phenomena that could unify many previously-proposed mechanisms to enable rapid formation of a well-defined dynamical center, "trapping" of BHs, efficient and rapid migration of BH pairs to the center, efficient gas capture, and hyper-Eddington accretion through a hypermagnetized accretion disk.We show this is able to take a BH from "seed" masses to truly supermassive masses in less than one Myr, providing a natural theoretical scenario for high-redshift BH growth.

SIMULATION
The simulation is one of the suite in Shi et al. (2024), where the numerics are described in detail, so we only briefly summarize the salient properties here.We utilize the magnetohydrodynamics (MHD) code gizmo1 with MHD solved using the Meshless Finite Mass (MFM) method (Hopkins 2015;Hopkins & Raives 2016) and self-gravity using a tree-particle mesh method with self-consistent adaptive force softenings (Hopkins et al. 2023).In addition to self-gravity and (ideal) MHD, the code includes radiative cooling and heating, star formation, and stellar feedback, following the treatments in Grudić et al. (2018b); Shi et al. (2021Shi et al. ( , 2023Shi et al. ( , 2024)).In detail stars form from gas which is dense, molecular, rapidly-cooling and Jeans-unstable below the resolution limit and stars subsequently influence the medium via 5-band (FUV through IR) radiation (heating and photon momentum), stellar mass-loss, and supernovae, with rates tabulated assuming a well-sampled universal stellar initial mass function per the FIRE-2 implementation (Hopkins et al. 2018a,b) of the Feedback In Realistic En-Figure 1. Visualization of the simulation.We show the gas morphology (column density in logarithmic scales) and positions of star particles (colored dots; only a subset is displayed for aesthetic purposes) and two selected BHs (black stars) at different stages of the simulation.Panel (a): Early star formation (t ≪ t ff ), where seed BHs wander in the GMC.Panel (b): Hierarchical cluster assembly and mergers (t ∼ t ff ), where bursty star formation happens but the GMC is bound due to strong self-gravity.Massive sub-clusters form and one of the BHs is "captured" by (gravitationally bound to) a cluster (red dots in the inserted zoom-in panel, where the radius of the red circle equals 5 half-mass radii of the cluster).Panel (c): Proto-bulge formation and steady BH accretion (t ≳ 2 t ff ), where subclusters merge into the final massive cluster, and the seed BH captured by the subcluster is migrated to the center.The inserted zoom-in panel shows the gas morphology around the BHs, which suggests convergent inflow and disk structure.
In order to achieve a controlled experiment, we consider a non-cosmological idealized experiment of an initially turbulent, magnetized cloud complex, with some initial BH seed population.The complex has initial mass 10 8 M ⊙ (all gas), spherical radius 50pc, so surface density Σ gas ∼ 10 4 M ⊙ yr −1 , with initial turbulent virial parameter of unity (and initial turbulent/zero mean magnetic fields with magnetic energy equal to 1% of gravitational energy), chosen to be representative of extreme but well-known conditions in starburst galaxy nuclei or high-redshift star-forming gas clumps in massive galaxies (Swinbank et al. 2012;Izumi et al. 2016;Scoville et al. 2017;Tacconi et al. 2018).We randomly distribute BH seeds spatially throughout the cloud, with initial masses uniformly distributed in 2 < log 10 (M seed /M ⊙ ) < 4 and isotropic random velocities with the cloud virial velocity dispersion.Once the simulation begins, seeds can capture gas that is strictly bound to the BH (including thermal, magnetic, and kinetic energies) interior to the BH sink radius (Bate et al. 1995;Shi et al. 2023), which here is ∼ 0.1 pc (sufficient to easily resolve the Bondi radius of the relevant gas for accretion even at warm molecular temperatures, and the tidal radius of clumps/cores to the BHs, though of course much larger than the BH gravitational radii) comparable to state-ofthe-art "hyper-resolution" simulations from ISM/galaxy scales (Anglés-Alcázar et al. 2021).Since this still leaves the inner accretion disk unresolved, we use a sub-grid model (for detailed description see Shi et al. 2024) to approximate the rate of this captured gas accreting onto the BH's inner disk ( ṀBH,0 ≈ M gas, disk /t dep ) with some disk depletion time following a Shakura & Sunyaev (1973)-like scaling.Apart from the mass accretion rate ṀBH that eventually contributes to the BH's mass growth, some mass Ṁw = η w ṀBH = η w /(1 + η w ) ṀBH,0 is assumed to be ejected in winds and bipolar jets following the methods in Torrey et al. (2020); Su et al. (2021) with effective large-scale (emergent from un-resolved sales) velocity v j , while a fraction ϵ r emerges as radiation with the Shen et al. (2020) template spectrum.We set η w = 1, v j = 400 km s −1 , ϵ r = 10 −5 .These are chosen to be representative of properties of super-Eddington "slim disk" models accreting at the hyper-Eddington rates we will find (Watarai et al. 2000;Jiang et al. 2019), but a parameter study of the effects of these choices is the subject of Shi et al. (2024).We run the simulation for 1.8 Myr (a few times the initial global complex free-fall time), at which point we show the system is in quasisteady state, with most star formation and BH growth complete.
To build a star sub-cluster merger history in postprocessing, we first identify gravitationally bound star clusters in each snapshot. 2To check the inheritance relationship between cluster i in one snapshot and cluster j in another one, we use a "similarity measure" which is defined as where n i (n j ) is the number of particles in the cluster i (j), and n ij is the number of common particles between the two clusters.For a cluster in one snapshot, we iterate over clusters in the next snapshot to rank and determine the most possible successor (above a preset threshold).
A "merger tree" is constructed after iterating over all clusters and all snapshots (densely sampling the local dynamical times) in the simulation.

PROTO-BULGE FORMATION AND BLACK HOLE CAPTURE AND MIGRATION
Fig. 1 shows the gas and stellar morphology at different evolution stages in the simulation.The simulation starts with a smooth density profile with initial turbulence seeded.Then gravitational collapse occurs and the medium becomes more turbulent, creating dense clumps, shocks, and filaments.Star formation begins in these dense regions well before one global free-fall time.The high surface density (or equivalently "gravitational pressure" or acceleration) scale is critical: it means the gas is not fully disrupted by "early" (pre-SNe) stellar feedback, despite turning most of its initial gas mass into stars, owing to the strong self-gravity.That in turn means that the bound-cluster star formation efficiency in the various sub-clumps is large (i.e.we form dense, bound subclusters, and not an unbound, low-density open cluster/association), consistent with both analytic expectations and results from a wide variety of simulations (Fall et al. 2010;Grudić et al. 2019Grudić et al. , 2020;;Wada et al. 2009;Hopkins et al. 2016Hopkins et al. , 2022)).
Given the high star formation efficiency, we see these form into hierarchical sub-clusters that merge rapidly, as in local young massive clusters (Grudić et al. 2018a;Li et al. 2019;Guszejnov et al. 2022).In Fig. 2 we show the trajectories of (sub)cluster centers, following their merger history.We can see that clusters merge rapidly after their formation and migrate to the center.We also plot trajectories of the two BH seeds with the most significant mass accretion (same as those in Fig. 1), Figure 2. BH capture, migration, and growth as subclusters merge.Left: Trajectories of the centers-of-mass of the two subclusters (1,2; solid) which capture the two BHs (1,2; dashed) from Fig. 1.The BHs are "captured" (become bound to the subcluster) at the times labeled, well before they reach the center, and are carried to the center by the merging subclusters.Minor subclusters merged into these two massive clusters are tracked with gray lines.Upper right: Mass evolution of the subclusters and BHs.Clusters acquire mass initially through star formation then hierarchical merging, which completes (and the system relaxes) by ∼ 1 Myr.The BHs grow in early "jumps" as they accrete clumps near their subcluster centers, then at a steady rate ∼ 1 − 5M⊙yr −1 after ∼ 0.7 Myr from the nuclear disk.Bottom right: Boundedness/sinking of each BH in its sub-cluster.We show the spatial (|r rel |; solid) and velocity (|v rel |; dashed) separation between BH and its subcluster center-of-mass, normalized by the subcluster half-mass radius r h or escape velocity −2Φc(r rel ).The BHs sink and become more confined as their subclusters grow and merge.
Through evolution, each cluster also changes its mass M c and half-mass radius r h .For the two clusters highlighted in Fig. 2, mass growth, and accompanying deepening of their potential well, enhance capture of their respective BHs.This is illustrated in the bottom right panel of Fig. 2, where we track the relative distance (r rel ≡ r BH − r cluster , normalized to r h , where r cluster is the cluster's gravitational center) and velocity (v rel = dr rel /dt, normalized to the escape velocity at the BH position v esc = −2Φ c (r rel ), where we calculate the gravitational potential Φ c (r rel ) ≡ − i Gm i /|r BH − r i | at the BH position from all stars in the cluster at that time).We see distance and velocity decay in ≲ 1 Myr (the same is true if we measure the specific angular momentum), for both the BHs within their sub-clusters and the subclusters with respect to the parent cloud/cluster.But we see that capture is not just due to a chance encounter between a BH and subcluster at low speeds (which would be rare).Rather Fig. 2 shows that the two BH-cluster pairs are not gravitationally bound initially, and indeed the relative velocities of the BHs to clusters actually increase (from ∼ 50 to ∼ 200 km s −1 ) from early to late times, but the ratios of kinetic to potential energy become < 1 and the BHs are captured at 0.4 and 0.7 Myr, the same times that their clusters grow rapidly in mass, increasing |Φ c | and the subcluster escape velocity rapidly relative to the (relatively weaky-evolving) v rel .
The final merged star cluster has a mass of ∼ 5 × 10 7 M ⊙ , but with a dense central mass concentration (half-mass radius of ∼ 3 pc) and isothermal-sphere like density profile ρ * ∝ r −2 steepening to r −4 at larger radii, after undergoing violent relaxation from multiple hierarchical subcluster mergers.This means the characteristic circular velocity is large, ∼ 200 km s −1 (see §4), and the cluster lies in-between massive typical nearby star clusters and galactic bulges in most structural and kinematic scaling relations (where the two are known to form a continuous family; see Hopkins et al. 2010;Grudić et al. 2019).Indeed, these are very similar stellar masses, sizes, surface densities, velocity dispersions, and density profiles to known massive nuclear star clusters and ultracompact dwarfs observed (Geha et al. 2002;Walcher et al. 2005;Haşegan et al. 2005;Evstigneeva et al. 2007).So given the initial gas and final stellar properties, this suggests a picture where these "proto-bulges" form from sufficiently high-density gas complexes in high-redshift massive galaxy progenitors, which rapidly fragment into sub-clusters that can capture (even relatively modest mass) BH "seeds" and hierarchically merge in a dynamical time (as argued from some cosmological simulations in e.g.Ma et al. 2018Ma et al. , 2020)), forming a deep, centrallyconcentrated potential which traps and brings together the BHs.
It is worth mentioning that hierarchical star cluster assembly is a process with ample dynamics that may produce massive seed BHs.For example, Rantala et al. (2024) found that the process may boost the formation of very massive stars (VMS) through repeated star mergers, which can be progenitors of BH seeds.Still, for gas-rich environments which this work is focused on, BH accretion can play an important role in the "protobulge" system.We expand on these details in the next section.

SUSTAINED HYPER-EDDINGTON ACCRETION THROUGH A HYPER-MAGNETIZED DISK
As the clusters begin to dynamically "settle" after ∼ 0.4 − 0.7 Myr, we see an extended accretion disk appear and grow "inside-out," beginning just outside the sink radius until ultimately extending to ≳ 10 pc scales.Fig. 1 plainly shows the disk in the last inserted panel, while Fig. 3 shows the magnetic field structures.We overplot the magnetic field vectors, which are clearly dominated by a dominant mean toroidal/azimuthal field, though there are also poloidal fields vertically threading through the disk.This results in the nearlyconstant growth shown in Fig. 2 at ∼ 1 − 10 M ⊙ yr −1 for the remaining ≳ 1 Myr.From this and the BH masses, we see this is well above the Eddington limit which ranges from ∼ 0.006 − 0.04 M ⊙ yr −1 over the same time.We also see the disk extends well outside the BH radius of influence We summarize the gas disk properties in Fig. 4. We confirm both that it is indeed a disk (primarily rotation supported with ⟨v ϕ ⟩ ≈ V total c = GM tot (< R)/R out to ≳ 10 pc), and that it clearly appears to be an example of the "hyper-magnetized" or "flux-frozen" disks proposed in Hopkins et al. (2024a,b,c).The defining properties of these disks, all of which we see in Fig. 4, are that magnetic pressure dominates in the midplane (plasma β ∼ P thermal /P B ∼ c 2 s /v 2 A ∼ 10 −3 here), 3 coming primarily from a mean toroidal field (with ⟨B ϕ ⟩ much larger than either the mean radial/poloidal fields B R , B z or turbulent fields), with trans (or weakly sub) Alfvenic turbulence (v turb ∼ v A here).One can verify from Fig. 4 that the accretion rates through the disk are approximately constant (∼ 1 − 10 M ⊙ yr −1 ) with R, with values consistent with the measured mean radial ⟨v R ⟩, and Maxwell+Reynolds stresses directly extracted from the simulation (note the expected ⟨B R ⟩−⟨B ϕ ⟩ anticorrelation, discussed in Hopkins et al. 2024b, which ensures a strong Maxwell stress of the desired sign) as well as that the (thick) vertical height behaves as expected from magnetic+turbulent support (H/R ∼ σ z /V c ∼ v A /V c ∼ 0.3 − 1), and the midplane density scales accordingly (ρ ∼ Σ gas /2H).For a detailed analysis of Upper left: surface/column density, and mass accretion rate, and magnetic Toomre Qmag profiles of the disk (measured in cylindrical annuli).We also plot the predicted surface density Σ(R) for the analytic hyper-magnetized disk models in Hopkins et al. (2024c), assuming an accretion rate Ṁ /ϖν = 30 M⊙/yr.Lower left: mean rotation velocity v ϕ of the gas measured from the simulation (with the gray band representing the standard deviation), and the total circular velocity ≡ G Mtot(< R)/R with the contribution ≡ G Mi(< R)/R the mass of stars, BHs, and gas.Upper right: average radial and azimuthal velocities (dashed if negative) of the gas in the disk, and Alfven speed (vA), thermal sound speed (cs), and root-mean-square turbulent velocities (v turb ).Lower right: strength of the mean radial, toroidal, and poloidal magnetic field (dashed if negative).The overall magnetic field strength is in agreement with the flux-freezing assumption (|B| ∝ ρ 2/3 , dotted ).
these properties in such disks, we refer to Hopkins et al. (2024b).
The magnetic field strengths appear to be consistent with growth via flux-freezing/advection of magnetic flux from the diffuse gas into and throughout the disk: B ϕ rises approximately as r −1 from a couple hundred µG (the value of the turbulent field in the initial conditions, comparable to that observed in gas with a mean density similar to the initial cloud value n ∼ 10 4 cm −3 in the ISM; see Crutcher 2012;Ponnada et al. 2022) to ∼ 0.1 G on sub-pc scales.If we follow the componentwise analysis in Hopkins et al. (2024b), this amplification as well as the B ϕ − B R anticorrelation are consistent with the mean-field Lagrangian induction equation, as the gas is captured and compressed (and the fields "wound up" toroidally) as it accretes. 4As well, we overplot the specific analytic model dependence of ρ(R), B(R) from Hopkins et al. (2024c), after modifying their model to account for the non-Keplerian potential seen here (below), and see reasonable agreement.
As in Hopkins et al. (2024a), the strong magnetic field appears to stabilize the disk against fragmentation and star formation.We do not see fragmentation and the star formation rate after the disk forms is small (Fig. 2) compared to BH growth rates.Indeed, from Fig. 4 we can see that the Toomre Q parameter accounting for magnetic support (Q mag ≡ v A κ/πGΣ) is ≫ 1 at all 4 We can also verify that our assumption of ideal MHD is valid.
The microphysical Spitzer/Braginskii conduction/viscosity coefficients are tiny under these gas conditions (equivalently particle mean free paths are extremely small compared to ∼ H).And comparing non-ideal Ohmic/Hall/ambipolar resistivities (e.g.Kunz & Mouschovias 2009) to the turbulent diffusivities ∼ vt H shows that the microphysical non-ideal terms are smaller by at least ∼ 5 − 10 orders of magnitude, owing to a combination of the relatively low gas densities and high ionization fractions f ion ∼ 0.01 given the extreme stellar irradiation of the lowdensity disk.
radii here (typically ∼ 20 − 50), while if we included only thermal pressure (Q therm ≡ c s κ/πGΣ) we would obtain Q therm ∼ 0.1 ≪ 1 at all radii.The thermal properties themselves, while unimportant for the disk structure since β ≪ 1 (and because we are well outside the near-horizon radii where most of the accretion disk emission originates), are reasonable.At these large (∼ pc) radii the disk self-heating flux ∼ (3/8π) Ṁ Ω 2 is small compared to the flux from stellar irradiation.Assuming a cluster mass profile dominated by a young (≲ 10 Myr) stellar population, with a well-sampled IMF, and standard dust composition as assumed in-code, the equilibrium dust temperature (equating the stellar flux absorbed by grains to their emission) is ∼ 400 K (V c /200 km s −1 ) 2/5 (R/pc) −1/5 , which at the densities here should be in rough equilibrium with the gas temperature.Indeed this appears to reasonably predict the thermal properties in Fig. 4, though we note the disk is somewhat multi-phase with a mix of warm and cool atomic, cold molecular gas, and very compact HII regions around some stars.
There are some notable differences from the simulation presented in Hopkins et al. (2024a).First, the disk here forms under very different circumstances and in a very different regime of parameter space (in terms of BH mass, gas/galaxy properties, accretion rate, initial absolute magnetic field strength, etc.), and numerically the simulations here adopt a different MHD solver, treatments of star formation and stellar feedback (FIRE versus STARFORGE; Grudić et al. 2021a), resolution, BH feedback treatment, and initial conditions.All of which argues that such disks should not be uncommon in highaccretion rate BHs.Second, and perhaps most notably, the disk in those papers formed entirely within the ROI, where we see it extend well past the ROI here at all times.Fig. 4 explicitly shows that stars dominate the potential outside ≳ 0.5 pc even at the final time (when the BHs are most massive), with a somewhat flat, then again weakly declining V c (R).This implies that there is nothing uniquely special about a Keplerian potential required to obtain such a disk, so long as there is a well-defined centrally-peaked potential, and the magnetic fields, gravitational stability, and other properties of the gas permit.Indeed we can still apply the analytic models for the predicted structure of such disks derived in Hopkins et al. (2024c), if we replace Ω in their scalings (which they took to be Keplerian) with the appropriate value here Ω = V c /r.If we take their default model and make this replacement, we predict −∂ ln B ϕ /∂ ln R ∼ 8/9−4/3 and −∂ ln ρ/∂ ln R ∼ 4/3−2 (with the shallower slopes around ∼ 0.5 − 3 pc, where V c ∼ R 0 deviates most from Keplerian), consistent with the behavior in Fig. 4.This also has implications for the disk heating/thermal structure (above).Third, it is worth noting that this is formally a circum-binary disk, given the two nearly equal-mass BHs (Fig. 2).
It is striking that the two BHs, once brought together at ≳ 0.7 Myr, rapidly converge to nearly identical mass growth rates in the interior of the circumbinary disk.This effect has been seen before in circumbinary disk simulations (Lai & Muñoz 2023) and may explain populations of "twin stars" (El-Badry et al. 2019), but will clearly have important implications for BH mergers.

CONCLUSIONS
We present a scenario for extremely rapid (∼ 1 Myr) transition from seed-to-supermassive BHs in a forming proto-bulge or massive nuclear star cluster, with efficient growth and pairing of BH binaries.This unifies many theoretical ideas which have been proposed for different aspects of the problem, demonstrating for the first time all of them operating in concert in an MHD-thermochemical-star formation simulation.Our key conclusions, which define the steps in this scenario (also see Fig. 5), unfold as follows: 1.In a massive gas cloud complex of ∼ 10 8 M ⊙ and surface density Σ ∼ 10 4 M ⊙ pc −2 , typical of starburst nuclei or "clumps" in massive high-redshift galaxies, fragmentation and star formation are very efficient and stellar feedback is weak, so many bound sub-clusters rapidly form.A small fraction of low-mass BH seeds can grow essentially "stochastically" to IMBH masses by encountering dense, small clumps of gas.
2. Those IMBHs are efficiently "captured" by their adjacent subclusters during this process.As the subclusters grow and merge the BHs become even more effectively confined.
3. Subclusters merge hierarchically in essentially one dynamical time of the complex (< 1 Myr).This builds a proto-bulge, with properties similar to observed massive, dense nuclear star clusters and UCDs, with a sharply centrally-concentrated (isothermal) stellar mass profile.
4. Carried by their subclusters, the trapped BHs migrate to the center in the same time, orders of magnitude faster than their dynamical friction time if they were isolated.This also produces efficient "pairing" down to < 1 pc scales of multiple BHs.
5. The proto-bulge provides a deep, centrallyconcentrated potential which allows gas to be globally captured and retained in the system, with Pre-star forming giant molecular cloud complex with seed BHs.
Stellar feedback destroys the gas reservoir and prohibits seed BH growth in the complex.a coherent center and angular momentum, so it forms a large-scale accretion disk extending well beyond the BH radius of influence (to ≳ pc scales).Accretion transitions from stochastic/turbulent clump-clump encounters to disk feeding.
6.As gas is captured with modest pre-existing turbulent magnetic fields, the toroidal field lines are stretched and form a magnetically-dominated, flux-frozen disk.The disk pressure comes from the toroidal field, with plasma β ∼ 10 −3 in the midplane.This stabilizes the disk against catastrophic fragmentation and star formation and allows it to exist at all, while simultaneously ensuring strong Maxwell and Reynolds stresses that support extremely large accretion rates ∼ 1 − 10 M ⊙ yr −1 for ≳ Myr timescales, orders of magnitude larger than the usual Eddington limit.
7. The disk rapidly builds the BH from IMBH to SMBH masses, reaching ≳ 2 × 10 6 M ⊙ in ∼ 1 Myr.Moreover if pairing has occurred, the inner BH binary appears to undergo the expected process whereby the secondary accretes more rapidly until the two BHs are brought to nearly equal mass, and they appear to grow at nearly identical rates, promoting near equal-mass BH-BH mergers.
This scenario therefore provides a natural explanation for the origins and rapid formation of high-redshift quasars in the early Universe.It also has critical implications for mergers of supermassive BHs in the LISA bands, and if the processes above can operate at even higher mass scales (∼ 10 9 M ⊙ ), it could also provide an explanation for the seemingly-high rate of inferred SMBH mergers from pulsar timing (Agazie et al. 2023).We also note that we see a hyper-magnetized, fluxfrozen disk akin to those in Hopkins et al. (2024b) but under very different physical conditions and with distinct numerical methods and resolution and detailed physics treatments (as did e.g.Gaburov et al. 2012) which suggests these are robust and can indeed power super-Eddington mass accretion rates.It also demonstrates that one can quickly reach true supermassive BH masses without invoking exotic channels such as "direct collapse" (Volonteri & Rees 2005;Natarajan & Treister 2009;Mayer et al. 2010) of primordial gas clouds without fragmentation (indeed, here, fragmentation plays a key role promoting rapid BH growth) or other "new physics" ( § 1).
In followup work we hope to make more detailed observational predictions.But observations will be challenging, since we predict this occurs on short time and small spatial scales in high-redshift systems, with comparable stellar and BH accretion bolometric luminosity, and the medium is still dusty and obscured at these times owing to the high gas densities.Still it is immediately striking that the bolometric luminosities (∼ 5 − 50 × 10 44 erg s −1 , depending on the BH radiative efficiency), sizes (∼ 10 − 100 pc), broadly comparable mix of stellar-and-AGN luminosity, and colors (given the predicted dust column densities at large radii from Fig. 4) this simulation would exhibit during its later evolution (≳ 0.7 Myr, after coalescence and disk formation) are all remarkably consistent with the observed properties of the unexpected population of "little red dots" seen by JWST at redshifts z ∼ 4 − 10 (compare e.g.Kokorev et al. 2024;Andika et al. 2024;Kocevski et al. 2024).
Interestingly, the BHs here end up well above the M BH − M * (bulge mass), but well below the M BH − σ * (velocity dispersion) correlations observed at z = 0 (Kormendy & Ho 2013;McConnell & Ma 2013).This may therefore provide an explanation for apparently "overmassive" BHs observed at high redshift (and some lowredshift outliers; Seth et al. 2014;Trakhtenbrot et al. 2015;Walsh et al. 2016;Liepold et al. 2020), but also implies typical systems must grow both in bulge mass and BH mass before becoming "typical" z = 0 massive classical bulges (as indeed predicted for more classical feedback-regulated models of the evolution of these properties; see Hopkins et al. 2008b,a).Of course, there are many caveats and limitations to the simulations here which merit further study.Our resolution is limited: future work with super-resolution approaches like those in Anglés-Alcázar et al. ( 2021) would enable robust mapping of the behavior interior to the radius of influence, and connection to GRMHD simulations like Kaaz et al. (2023) to horizon scales.Improving the star formation model with approaches like Grudić et al. (2021a) would free us from the universal IMF assumption and allow for actual predictions of exotic IMFs (and their effects) in high-redshift systems.In order to consider a "controlled" experiment, we utilize idealized initial conditions with pre-existing gas and magnetic fields and BH seeds, which we could improve by following these systems in fully-cosmological simulations with self-consistent seed BH formation models.And we consider only one, simplified (and relatively "weak") model for "feedback" (in the form of jets, winds, and radiation) from BH accretion (motivated by the study in Shi et al. 2024 and slim disk models in Watarai et al. 2000;Madau et al. 2014).But given that the hyper-magnetized accretion disk here appears to be qualitatively, fundamentally distinct from classical α-disks on which these feedback models are based, it is important to understand the horizon-scale (hence radiation, jet and wind) properties of such disks in order to better inform next-generation sub-grid feedback models.

Figure 3 .
Figure 3. Magnetic field lines (near the midplane of each box) and gas morphology (logarithm of the column density) at the late stage of the simulation.The main panel suggests the large-scale magnetic field is advected to smaller scales due to gravitational collapse throughout the simulation.The inserted zoom-in panels focus on the disk structure, adapted to face-on (left panel) and edge-on (right panel) views.The face-on view shows a dominant toroidal magnetic field spiraling inward anti-clockwise, while the edge-on view shows the (weaker) poloidal magnetic field lines threading through the disk structure.

Figure 4 .
Figure4.Dynamics and magnetic field of the disk.Upper left: surface/column density, and mass accretion rate, and magnetic Toomre Qmag profiles of the disk (measured in cylindrical annuli).We also plot the predicted surface density Σ(R) for the analytic hyper-magnetized disk models inHopkins et al. (2024c), assuming an accretion rate Ṁ /ϖν = 30 M⊙/yr.Lower left: mean rotation velocity v ϕ of the gas measured from the simulation (with the gray band representing the standard deviation), and the total circular velocity ≡ G Mtot(< R)/R with the contribution ≡ G Mi(< R)/R the mass of stars, BHs, and gas.Upper right: average radial and azimuthal velocities (dashed if negative) of the gas in the disk, and Alfven speed (vA), thermal sound speed (cs), and root-mean-square turbulent velocities (v turb ).Lower right: strength of the mean radial, toroidal, and poloidal magnetic field (dashed if negative).The overall magnetic field strength is in agreement with the flux-freezing assumption (|B| ∝ ρ 2/3 , dotted ).

Figure 5 .
Figure 5. Illustration of this scenario of seed BH capture, growth, migration, and pairing in the star-forming GMC complex.The evolution of the star-forming complex is largely determined by the surface density, while low-surface-density clouds are expelled by stellar feedback, halting possible seed BH accretion (upper left → lower left).In high-surface-density clouds, seed BHs may undergo steps (a) → (e), corresponding to the qualitative trend of star formation and BH mass growth (lower middle panel; also see Fig. 2).Image credits -M42 (upper left panel ): NASA, ESA, M. Robberto (Space Telescope Science Institute/ESA) and the Hubble Space Telescope Orion Treasury Project Team; M60-UCD1 (panel e): NASA, ESA, and the Hubble Heritage (STScI/AURA).