Heavy Seeds and the Multi-scale Fueling of Early Supermassive Black Holes and Little Red Dots

Introduction

How do supermassive black holes (SMBHs) and Little Red Dots (LRDs) emerge in the early Universe [e.g. 134, 153, 65, 108, 87, 2] (Y. Koyanagi et al. 2024)? This question has further risen in prominence with the recent discoveries from the James Webb Space Telescope (JWST). SMBHs/active galactic nuclei (AGN) have been found in larger numbers and at higher redshifts than previously expected [18, 107, 111, 142, 24, 87, 114, 44, 46, 143], along with the revelation of the mysterious LRDs, exhibiting both AGN (e.g., broad spectral lines) and non-AGN like (e.g., lack of X-ray emission and variability) characteristics [2, 87, 99].

Various theoretical models have been proposed to explain high-redshift SMBHs and LRDs, some overlapping. To produce quasars with their large masses ($\gtrsim10^8 M_\odot$) at $z \gtrsim6$ [7, 160, 41] (X.-B. Wu et al. 2015), extreme conditions have to be invoked. One possibility is that the initial seed black holes (BHs) formed with somewhat higher masses, compared to locally observed stellar remnants at few solar masses. The first, or Population III (Pop III) stars, born out of metal-free pri mordial clouds, are expected to have a top-heavy initial mass function (IMF; A. Stacy et al. [139]; S. Hirano & V. Bromm [56]; M. A. Latif et al. [94]), and produce more massive BH seeds around $\sim100\,M_{\odot}$ (P. Madau & M. J. Rees [106]; A. Heger et al. [55]; J. L. Johnson & V. Bromm [77]). After formation, these light seeds will grow through gas accretion and BH mergers (M. Jeon et al. [75]; B. D. Smith et al. [135]; A. K. Bhowmick et al. [14]; A. Trinca et al. [145]; A. J. Porras-Valverde et al. [122]), requiring super-Eddington accretion to reach the observed high-$z$ SMBH masses (J. Jeon et al. [72]). Interestingly, such conditions of super-Eddington accretion could suppress X-ray emission/variability, in line with what is observed for LRDs (P. Madau & R. Maiolino [105]; H. Liu et al. [98]; A. Secunda et al. [132]; F. Pacucci & R. Narayan [116]).

Alternatively, to mitigate the need for sustained episodes of efficient accretion (Z. Haiman & A. Loeb [54]), heavy seed channels have been proposed, where initial seed masses already reach $\sim10^{5}\,M_{\odot}$ (e.g., F. Becerra et al. [9]). Direct Collapse Black Holes (DCBHs) constitute a well-established model of heavy seed formation – forming from the collapse of a massive gas cloud, with a supermassive star (SMS) and/or quasi-star (QS; M. C. Begelman et al. [12]) as intermediate, short-lived stages. The gas cloud must be extremely metal-poor ($Z \lesssim10^{-3}\,Z_{\odot}$, S. Chon & K. Omukai [30]) and requires rare conditions, prohibiting gas cooling to low temperatures to allow cloud-collapse without vigorous fragmentation (e.g., V. Bromm & A. Loeb [20]; G. Lodato & P. Natarajan [100, 101]; J. L. Johnson et al. [79]; J. H. Wise et al. [152]; L. Haemmerlé et al. [52, 51]; Y. Luo et al. [103]). Given the conditions required for their formation, DCBHs are expected to reside in low stellar-mass regions surrounded by a dense gas envelope, possibly producing LRDs (E. Cenci & M. Habouzit [21]; J. Jeon et al. [73]; F. Pacucci et al. [115]; D. J. Whalen et al. [151]). Similarly, the transitory SMS and QS phases have been proposed to be LRDs (J. Chisholm et al. [25]; V. Kokorev et al. [89]; M. C. Begelman & J. Dexter [11]; A. D. Santarelli et al. [130]).¹

Still, even for heavy seeds, efficient gas accretion, possibly even super-Eddington accretion, will be required to produce the observed early SMBHs with masses orders of magnitude above the proposed seed masses (e.g., J. Jeon et al. [70, 73]). BH mergers could facilitate BH mass growth at early times, but gas accretion will still be necessary later on (A. K. Bhowmick et al. [17]). However, self-consistently modeling BH growth into SMBHs is challenging due to the difference in scales covering their evolution. The BH accretion physics operates at sub-pc scales, whereas the cosmological environment providing the baryonic material for accretion evolves on larger kpc to Mpc scales. Cosmological simulations focused on galaxy-scale evolution cannot simultaneously resolve the small-scale BH physics, and thus require approximate subgrid models to represent BH evolution (e.g., J. Jeon et al. [71]; B. D. Smith et al. [135]; B. Liu & V. Bromm [96]; M. Volonteri et al. [147]).

For LRDs, an emerging model proposes that they host an accreting central object, likely an SMBH, surrounded by a dense, dust-free gas cocoon (K. Inayoshi & R. Maiolino [63]). The dense gas produces a red spectrum with a strong Balmer break and exponential profiles, well-matching the observed LRD spectra (A. de Graaff et al. [36]; V. Kokorev et al. [90]; A. J. Taylor et al. [142]; R. P. Naidu et al. [110]). The required hydrogen densities are extreme ($10^{8}-10^{10}\ {\rm cm}^{-3}$) and a massive reservoir of gas will be needed to sustain the dense gas envelope as the SMBH accretes (K. Inayoshi & L. C. Ho [62]; J. Jeon et al. [69]). Again, similar to the question of BH accretion, cosmological simulations cannot resolve the scale of proposed LRD configurations and cannot confirm the existence of the proposed dense gas cocoons around SMBHs. Instead, a simulation that covers the large-scale cosmological gas inflow to the small-scale BH environment will be required.

To address these questions of BH seed growth and the nature of LRDs, we perform high-resolution cosmological simulations in this work. Specifically, we utilize a cutout technique, where we excise a targeted region from a large-volume cosmological simulation, evolving it for short periods at a higher resolution (P. C. Clark et al. [32]; A. Stacy & V. Bromm [138]). We select regions near BH seeds that formed under different formation channels in the simulation, so that both the gas accretion and structure can be followed at extremely high, sub-pc resolution, elucidating the multi-scale nature of early Universe BH evolution and the emergence of LRDs, from the cosmic web to near-BH scales. Other works have also tried to resolve both the cosmological and near-BH scales simultaneously (e.g., P. F. Hopkins et al. [59]), but at a high computational cost. Our numerical approach, involving relatively low particle numbers, is computationally much less expensive, thus enabling a broader exploration of parameter space.

In this work, we use Planck cosmological parameters [121]: $\Omega_m = 0.315$, $\Omega_b = 0.048$, $\sigma_8 = 0.829$, $n_s = 0.966$, and $h = 0.6774$. We explain our simulation models and cutout setup in § 2, present our results in § 3, discuss the implications for SMBH and LRD assembly in § 4, and summarize the work in § 5.

Numerical Methodology

We use the GIZMO code [57], employing the Lagrangian meshless finite mass (MFM) hydro solver, and the gravity solver inherited from GADGET-2 [136]. We specifically use the modified version of GIZMO [96, 71, 70] with distinct sub-grid models for star formation (SF) and feedback; BH formation, accretion, and feedback; as well as primordial chemistry, cooling, and metal enrichment [67, 66].

Simulation setup

The simulation initial conditions are generated with the MUSIC code [53] at $z = 99$. We first run a parent simulation of 2 cMpc/$h$ in linear size with only dark matter (DM) particles with masses $8 \times10^3$ M$_\odot$. Using Rockstar [13] in post-processing to identify DM halos, we select a $\sim7 \times10^8$ M$_\odot$ halo at $z \simeq6$. We then identify a region that contains all the DM particles of the selected halo as the zoom-in region.

The zoom-in simulation has a linear size of 450 ckpc/$h$, includes baryons, and has an increased resolution with DM particle masses of $\sim800$ M$_\odot$ and baryon particle masses of $\sim150$ M$_\odot$. The high-resolution zooms are run to $z \sim4 - 5$ to examine the formation of SMBHs.

In the zoom-in simulation, we calculate the virial parameter $\alpha$ of the gas particles within 100 pc (in physical units) of each BH particle that formed across redshifts,

where $E_{\mathrm{th,pot}}$ is the thermal/gravitational potential energy of the cutout system, $c_s$ the gas sound speed, $R_{\mathrm{eff}}$ the weighted average of the gas and DM particle distance from the cutout center (BH location), $M_{\mathrm{gas}}$ the mass of the gas, and $M_{\mathrm{DM}}$ the mass of DM particles in the region. Lower values of $\alpha$ suggest a more bound gas system, and we choose, for a DCBH or a stellar BH, systems with high mass but also low $\alpha$, to consider strongly bound configurations so that optimistic BH growth can be followed.

The chosen system is then converted into a cutout containing the gas, DM, and central BH particle. The gas particles within the 100 pc region are split by a factor of 100 (see also [83]), so that they have masses of $\sim1$ M$_\odot$ and $\sim0.1$ pc minimum softening length resolutions. The split gas particles are placed randomly within the smoothing kernel of the original gas particle [32, 138]. We do not split the DM particles, to not introduce unphysical gravitational structure in the cutout not inherited from the original cosmological simulation. We further perform a run where the gas particles are split by a factor of 500, resulting in mass and length resolutions of $\sim0.2$ M$_\odot$ and $\sim0.02$ pc, respectively. In the cutouts, we disable SF and stellar feedback, as the gas particle masses are now too low for our respective models that have been designed for cosmological simulations with higher gas particle masses. For BH models, specifically accretion and feedback (see § 2.5 and 2.5), we perform runs where they are disabled, such that the BH only interacts gravitationally with the gas, and where they are enabled, to assess the effects of high resolution on BH growth and energy output.

The cutouts (see § 2.2 for details) are run for $\lesssim10$ Myr. We thus avoid unphysical edge effects, resulting from the artificial outer (vacuum) boundary conditions of the cutout region, where all external material has been removed. We can estimate the period over which this approximation is valid with the sound crossing time to traverse the cutout region. Given an average gas temperature of $\sim10^3 - 10^4$ K, corresponding to a sound speed of $\sim10^{-5}$ pc yr$^{-1}$, the distance corresponding to the 10 Myr period is $\sim100$ pc, half the size of the cutout region. We accordingly limit our analysis to $< 8$ Myr periods. Table 1 summarizes the simulation setup and Figure 1 and Figure 2 illustrate the wide range of length scales involved. We further verify the validity of the cutout methodology with an idealized test case of a collapsing spherical gas cloud, as described in Appendix § 7.

Simulation	Linear Size	DM particle mass	Gas particle mass	Gas minimum gravitational softening length
Parent	2 cMpc/$h$	$8 \times10^3$ $M_\odot$	N/A	N/A
Zoom	0.45 cMpc/$h$	800	150	10 cpc/$h$
Cutout	200 pc	800	1.5	0.1 pc
Cutout_Highres	200 pc	800	0.3	0.02 pc

Cutout selection

We run multiple cutout simulations for two select BH seeds while varying the included BH subgrid models (accretion and feedback). We specifically consider one stellar remnant BH and one DCBH seed, selected from the zoom-in simulation to represent cases where the nearby gas is closest to being self-gravitating, as determined by Eq. (1). The DCBH cutout had $\alpha\simeq0.3$ and the Stellar cutout $\alpha\simeq0.2$. We note that these values for the virial parameter are lower limits, not including the modest contribution of supersonic, turbulent motions in the velocity dispersion (which would result in $\alpha\sim1$). For both seeds, we run cutouts with BH accretion turned off and on. For the DCBH case only, we further perform runs with increased resolution by a factor of 5 (DCBH_Highres run), and with BH feedback enabled as well (DCBH_Feedback run). We do not perform corresponding runs for the stellar remnant seed, as we find minimal changes in the stellar BH evolution and mass growth (see § 3). In Table 2, we summarize the properties of the cutout simulations that have been carried out.

Cutout	Initial Redshift	Gas Particle Number	Total Gas Mass	Initial BH Mass	BH Accretion	BH Feedback
Stellar	10.1	$2.56 \times10^5$	$3.75 \times10^5$ $M_\odot$	140	No	No
DCBH	9.1	$7.3 \times10^5$	$1.1 \times10^6$ $M_\odot$	$1.5 \times10^5$ $M_\odot$	No	No
DCBH_Accrete	9.1	$7.3 \times10^5$	$1.1 \times10^6$ $M_\odot$	$1.5 \times10^5$ $M_\odot$	Yes	No
Stellar_Accrete	10.1	$2.56 \times10^5$	$3.75 \times10^5$ $M_\odot$	140	Yes	No
DCBH_Feedback	9.1	$7.3 \times10^5$	$1.1 \times10^6$ $M_\odot$	$1.5 \times10^5$ $M_\odot$	Yes	Yes
DCBH_Highres	9.1	$3.7 \times10^6$	$1.1 \times10^6$ $M_\odot$	$1.5 \times10^5$ $M_\odot$	No	No

Stellar Models

For the zoom-in runs, we employ the SF and feedback modules developed in J. Jaacks et al. [67, 66]; B. Liu & V. Bromm [95, 96]. We include models for both the first, metal-free Pop III, as well as metal-enriched Population II (Pop II) stars. The Pop III model has been tested against extremely high-resolution simulations (e.g. J. L. Johnson et al. [78]; J. S. Ritter et al. [127, 128]), and the Pop II model calibrated with observed cosmic star formation histories across $z \sim5\text{--}15$. Individual stars cannot be resolved with the resolution reached in the zoom-in simulations, such that each stellar particle represents a stellar population, characterized by Pop III and Pop II IMFs from J. Jaacks et al. [67, 66].

We further include stellar feedback in the form of Lyman-Werner (LW) radiation, photo-ionization heating, mechanical winds, and supernova (SN) explosions (B. Ciardi & A. Ferrara [31]). Sufficiently strong LW flux enables DCBH formation in the simulation by destroying molecular hydrogen to limit gas cooling at low temperatures. If gas collapse to form stars is delayed until the DM halo grows to be more massive than the atomic-cooling threshold ($> 10^7 M_\odot$), vigorous fragmentation into less massive stars is prevented, and the high rates of gas inflow may instead trigger the formation of massive BH seeds (see also B. Liu et al. [97]). For LW and photo-ionization effects, we include both global and local contributions. For full details of the SF and feedback models, we refer the reader to J. Jaacks et al. [67, 66]; B. Liu & V. Bromm [96].

Black hole seeding

As BH seeds, we consider Pop III stellar remnants and DCBHs. For Pop III remnants, we use the massive seed model of B. Liu & V. Bromm [96], where all of the stellar particle mass is converted into the BH seed mass. When a star particle reaches the end of its life after 3 Myr, we assume that the stars of the particle will merge to form a single BH, effectively converting the stellar particle to a BH sink particle.

For DCBH seeding, we impose multiple conditions to represent the dense, hot, and metal-poor environments necessary for formation (e.g., [20, 112, 5, 152, 29]). When a gas particle’s (hydrogen number) density reaches $n_{\mathrm{DCBH}} \ge2 \times10^{3}\,\mathrm{cm}^{-3}$, has a temperature of $7000 < T/\mathrm{K} < 10^{4}$, a metallicity of $Z < Z_{\mathrm{DCBH}} = 2 \times10^{-4}\,Z_\odot$, and H$_2$ abundance of $x_{\mathrm{H}_2} < 10^{-6}$, it is considered as a DCBH candidate. Furthermore, the gas particles in the neighborhood of the DCBH candidate (within the BH kernel) must also satisfy the temperature, metallicity, and H$_2$ abundance criteria, with the DCBH candidate being the densest gas particle among its neighbors. When all these conditions are met, the DCBH candidate is converted to a BH sink particle.

However, as the gas particle masses in the zoom-in simulation ($\sim150\,M_\odot$; § 2.1 and Table 1) are much lower than the fiducial DCBH mass of $\sim10^{5}\,M_\odot$ (F. Pacucci et al. 2017; F. Becerra et al. 2018b,a; S. Basu & A. Das 2019), the initial sink particle mass formed is much lower than the predicted DCBH mass. Thus, soon after formation, we force the DCBH sink particle to stochastically swallow nearby gas particles ([137, 96]), around 1,000 of them, to reach $\sim10^{5}\,M_\odot$ while maintaining mass conservation.

To enable such rapid initial growth through aggressive gas swallowing, we use two different mass values to characterize the sink particle. The first is the internal BH mass $M_{\mathrm{BH}}$ set at formation to 1,000 times that of the original gas particle mass. This mass is used to determine BH accretion and feedback parameters in the simulation. The second is the dynamical mass $M_{\mathrm{dyn}}$ of the sink particle, set at formation to the progenitor gas particle mass. This mass is used to determine the dynamical interactions of the sink particle with other particles in the simulation. $M_{\mathrm{dyn}}$ of the sink particle is only updated when the sink particle can swallow a gas particle. Thus, when $M_{\mathrm{BH}}$ of the particle is set to be larger than $M_{\mathrm{dyn}}$ at formation, the sink particle will try to swallow nearby gas particles to match the two masses, while maintaining mass conservation (see § 2.5). The criterion that the environment around the DCBH candidate also satisfies DCBH formation conditions establishes that such sum of gas particles is not unphysical, as it results in the formation of DCBH seeds in high gas density environments, such that the seed can immediately swallow the nearby gas particles to reach DCBH masses.

Black hole accretion and growth

For BH accretion, we use the Bondi-Hoyle equation ([19]):

where $M_{\mathrm{BH}}$ is the mass of the BH and $\rho_g$ is the density of gas particles in the BH kernel. The kernel size is adaptive to include 64 gas particles nearest to the BH.

When BH accretion is enabled for the cutouts, we instead use the modified Bondi-Hoyle equation accounting for the angular momentum of the gas ([144, 96]), as the gas structure around the BH is well resolved. We first calculate the characteristic gas rotational velocity around the BH as $v_\theta= j/\epsilon_g$, where $j$ is the specific angular momentum of the nearby gas and $\epsilon_g$ is the gravitational softening length of the BH particle. When compared against the gas bulk velocity, $v_{\mathrm{bulk}}$, if $v_\theta\leq v_{\mathrm{bulk}}$, angular momentum is negligible and we use

where $v_g$ is the mass-weighted average velocity of the gas particles relative to the BH particle. If instead, $v_\theta> v_{\mathrm{bulk}}$, we use

We further limit the accretion rate to 10 times the Eddington accretion rate, defined as

where the radiative efficiency $\epsilon_0 = 0.125$ [113]. $M_{\mathrm{BH}}$ is updated every timestep $\delta t$ as $\delta M_{\mathrm{BH}} = \dot{M}_{\mathrm{Bondi}}\delta t$. $M_{\mathrm{dyn}}$ tries to match $M_{\mathrm{BH}}$ growth by stochastically swallowing nearby gas particles in the BH kernel.

We adopt a probabilistic procedure for BH sink particle swallowing of gas particles [137]. The probability that a gas particle within the BH sink particle kernel will be swallowed (accreted) is

where $w$ is the kernel weight of the gas particle and $\dot{M}_{\mathrm{BH}}$ the BH accretion rate. A random number $p$ is drawn, and if $p < p_{\mathrm{acc}}$ the gas particle will be swallowed.

Black hole feedback

We model thermal feedback from BHs [96, 71, 70], injecting energy into the gas particles in the BH kernel. The distributed energy for $\delta t$ is $\delta E = \epsilon_rL_{\mathrm{BH}}\delta t$, with the radiation-thermal coupling efficiency $\epsilon_r = 0.02$ [144]. $L_{\mathrm{BH}}$ is the BH luminosity defined as

with $c$ being the speed of light. The radiative efficiency $\epsilon_{\mathrm{EM}}$ is defined as

with $A = 100$ [113], asymptotically approaching $\epsilon_0$ for high accretion rates.

Results

After running the zoom-in and cutout simulations, we post-process the simulation output with the yt [146] and Meshoid [48] packages. We compare the gas environment at high resolution for BHs with different origins, see how BH accretion and feedback behaves close to the BH, and predict their spectra from the accretion and nearby gas environment.

Gas distribution near the Black Hole

We first examine the gas density and temperature evolution around the BHs in the cutouts after a few Myr. The gas density distribution is shown in Figure 3, where the DCBH resides within a high-density region, while the stellar remnant BH is located in a low-density, ambient gas environment. To further understand this behavior, in Figure 4 we display the hydrogen gas density distribution versus the distance from the BH sink particle after the cutouts have been run for $\sim4$ Myr. We find that the massive DCBH can affect the gas in the cutout, to reach densities up to $\sim10^7\ \mathrm{cm}^{-3}$ when approaching the resolution limit at $\sim0.1$ pc. In contrast, the light stellar seed cannot influence the gas significantly, so that the bulk of the gas is not influenced by the gravitational potential of the low-mass BH with no buildup of high-density gas as in the DCBH case.

Comparing the resulting gas temperature distribution in Figure 5, we find no significant differences in the temperature range between the DCBH and stellar remnant cutouts, both exhibiting average temperatures of $\sim10^3$ K. We note that for the DCBH case, in difference from the substantially higher densities reached, the temperatures of the gas coalescing close to the BH are similar to those of the more widely distributed gas. In fact, the DCBH cutout region contains higher temperature gas ($>10^4$ K) from SNe in the original Zoom simulation. Still, despite the presence of such hotter gas, the DCBH is able to trigger a dense gas buildup.

We note that the DCBH and Stellar cutouts are two distinct environments, drawn from the Zoom simulation at different redshifts and locations. Thus, some environmental dependence may exist in the gas temperature and density distribution. However, the two cutout systems are similar, having comparable masses with only a factor of $\sim3$ difference (see Table 2), and with similar values of the virial parameter ($\alpha$), implying that the overall thermal energy content of the two systems is also similar. Thus, the environmental differences should not be dominant in determining how the gas responds to the central BH.

To assess whether numerical resolution has any effect on the results, we examine the gas distribution in the

DCBH_Highres cutout. Figure 6 shows the gas density and temperature distributions versus the distance to the BH in the DCBH_Highres run. The gas behaves similarly to the lower resolution case, except that the gas density near the BH reaches even higher values, up to $\sim10^{8}\ \mathrm{cm}^{-3}$. The higher resolution is able to resolve the buildup of dense gas even closer to the central BH seed, but overall, the difference in resolution does not significantly affect the resulting gas distribution.

Black hole accretion and feedback

The cutouts discussed above disabled BH accretion and feedback to first examine how increased resolution affected gas particle interaction with the BH. Now, we enable BH accretion and/or feedback to see how the higher resolution affects BH evolution. Figure 7 shows the BH mass evolution and mass accretion rates for the DCBH_Accrete, Stellar_Accrete, and DCBH_Feedback cutouts.

As expected from the gas distribution near the BHs (§ 3.1), the DCBH_Accrete cutout grows significantly in BH mass, able to swallow the high density gas nearby in brief periods of super-Eddington accretion. Such bursty super-Eddington accretion has previously been seen in other high-resolution simulations of SMBH growth [e.g. 109, 124, 28]. Moreover, large-scale simulations that cannot fully resolve SMBH accretion have introduced artificial burstiness to match observational constraints on AGN luminosity functions [133]. In contrast, the stellar remnant experiences sub-Eddington accretion rates throughout, resulting in in-significant mass growth. Therefore, we find that only heavy BH seeds could facilitate such bursty growth, whereas the light stellar remnant BHs cannot draw in the nearby gas, thus never experiencing high accretion episodes during the cutout run period, with maximal Eddington ratios of $\dot{M}_{\mathrm{BH}}/\dot{M}_{\mathrm{Edd}} \sim10^{-2}$ (see also J. Jeon et al. [71]).

When we include BH feedback, the growth is suppressed even in the DCBH case. After a brief sub-Eddington accretion episode at $\sim2$ Myr, the nearby gas is evacuated from the vicinity of the BH and growth stops. However, the effects of BH feedback are likely overestimated in the cutout setting, as we do not properly model the radiation-hydrodynamical evolution of the initially compact ‘(thermal) feedback bubble’ into the highly-anisotropic surrounding medium, but simply inject feedback energies to the BH kernel (§ 2.6). In reality, photon trapping and advection may act to blunt the impact of feedback (see § 4.4 for a discussion of numerical limitations). In contrast, previous studies have suggested that in thick disks enabling super-Eddington accretion, the effects of feedback may be mitigated (K. Inayoshi et al. [61]; Y.-F. Jiang et al. [76]; K. Park et al. [118]). Thus, our DCBH_Accrete cutout might be more representative of the cases where

super-Eddington accretion is possible with less efficient feedback.

AGN Spectra

To infer how the DCBH_Accrete cutout will appear observationally, we use the Cloudy photoionization code, version C23.01 (M. Chatzikos et al. [23]). In previous work, we did not have information on the close-in gas distribution around a SMBH, given the limited resolution in the larger-scale cosmological simulation considered there. We, therefore, had to make simplifying assumptions in extrapolating to the unresolved scales, such as constant gas density or pressure (J. Jeon et al. [70]). In the cutouts, on the other hand, the nearby gas distribution is well resolved to sub-pc scales. We specifically derive Cloudy spectra for the gas density profile over the $\sim0.1$ to $\sim100\,\mathrm{pc}$ range of distances from the central SMBH (averaging at each radial bin), as found in the DCBH_Accrete cutout run at $\sim5.7\,\mathrm{Myr}$, when the highest accretion rate is encountered. We thus avoid arbitrarily assuming a constant gas density for the previously unresolved scales.

For the AGN part of the broad-band spectrum, we use the Cloudy model implementation of an AGN with a a multicolor disk UV component and a non-thermal X-ray component²:

where $h$ is the Planck constant, $k_{\mathrm{B}}$ the Boltzmann constant, and $T_{\mathrm{disk}}$ the SMBH accretion disk temperature. The $a$-parameter is chosen to fix the X-ray to UV flux ratio, where the non-thermal term $a\nu^{-1}$ is only included when $h\nu\ge1.36$ eV. For the disk temperature, we use the slim disk model of K.-y. Watarai [149] (see their Eq. 20) for super-Eddington accretion, determining $T_{\mathrm{disk}} \simeq10^6$ K. For the X-ray to UV ratio, we use the model of Y. Avni & H. Tananbaum [6] (see their Fig. 8) and determine a ratio of $\sim-1.4$ in log scale. We employ the SMBH accretion rate at 5.7 Myr and Eq. (7) to determine the AGN luminosity as $\sim10^{44.4}\,\mathrm{erg}\,\mathrm{s}^{-1}$. Lastly, we use the average gas phase metallicity of the

gas particles in the cutout, assuming that no significant metal enrichment has occurred since the cutout was initialized, with a value of $\sim0.02\,Z_\odot$.

The resulting spectra are shown in Figure 8. We further show in comparison the model spectrum from a zoom-in cosmological simulation with lower resolution ([70]), as well as normalized observed LRD spectra, including the BH* object ([110]), CAPERS-LRD-z9 ([142]), and RUBIES-EGS-49140, with spectral data from the CANDELS-Area Prism Epoch of Reionization Survey (CAPERS; [39, 37]). Whereas both simulation spectra show a red AGN continuum, the strong Balmer-break feature, similar to that of observed LRDs, can only be seen when the high-density gas in the vicinity of the central BH is properly resolved. Here, the required high resolution is achieved through the cutout methodology.

Discussion

On the basis of the BH evolution resulting from the cutout simulations, we consider the implications for the overall SMBH mass evolution and the possible origin of the compact LRDs.

Supermassive Black Hole Formation

To form the observed SMBHs with masses $\sim10^{8}\,\mathrm{M}_{\odot}$ at $z \gtrsim6$, the seed BHs must grow rapidly, at near-or super-Eddington rates (e.g., [134, 15, 70, 109]). We find that such growth is only possible for heavy seeds. As shown in Figure 7, only the DCBH cutout experiences significant bursty growth and super-Eddington accretion periods. Without having an initially large mass and deep gravitational potential well, the light BH seed could not gather the nearby gas (Figure 4), and the resulting growth was negligible.

This conclusion would imply that most of the observed high-redshift SMBHs originated from heavy seeds. While other channels such as BH mergers exist for BH mass growth besides gas accretion, previous studies have found that accretion significantly contributes to SMBH growth even when considering high merger rates ([16, 17]). Furthermore, there are multiple pathways for heavy seed formation in addition to the DCBH scenario, which is our main focus here, such as hierarchical mergers (see § 1 and § 4.3). If the various heavy seed pathways place the BH into a largely similar environment, allowing for a period of time to lose memory of the specific seeding process, our cutout simulations would then be expected to exhibit similar behavior for the subsequent evolution, provided that the main factor determining the gas distribution and the resulting BH accretion is the BH mass, not the seed formation mechanism. However, the DCBH pathway naturally places a massive seed in a pristine and cold gas environment lacking stars, and may be more effective in inducing dense gas buildup and super-Eddington accretion. Scenarios like hierarchical mergers, on the other hand, will form heavy seeds within dense stellar environments, such that the gas may be heated and/or expelled by stellar feedback even in the presence of a deeper BH gravitational potential well.

To derive a rough estimate for the heavy BH seed mass needed to allow for efficient gas accretion, we consider the Bondi radius, giving the BH’s sphere of gravitational influence:

To calculate $c_{\mathrm{s}}$, we assume an average gas temperature of $10^{3}\,\mathrm{K}$ near the BH (Figure 5). We obtain the mass enclosed within the BH’s region of gravitational influence by assuming a sphere of radius $R_{\mathrm{Bondi}}$, containing hydrogen gas with average number density $100\,\mathrm{cm}^{-3}$ (Figure 4). We find that for a $10^{5}\,\mathrm{M}_{\odot}$ BH, the mass within the Bondi radius is $\sim5 \times10^{6}\,\mathrm{M}_{\odot}$, while for a $100\,\mathrm{M}_{\odot}$ BH the mass enclosed is only $\sim5 \times10^{-3}\,\mathrm{M}_{\odot}$. For the intermediate case of a $10^{4}\,\mathrm{M}_{\odot}$ BH, the enclosed mass is $\sim5 \times10^{3}\,\mathrm{M}_{\odot}$. We recognize that the enclosed mass estimate for the $10^{5}\,\mathrm{M}_{\odot}$ BH is larger than the mass contained in the cutouts (Table 2), so that the single gas density and temperature approximations will break down at larger scales. Still, this argument implies that to enable efficient gas accretion onto a BH, heavy BH seeds of $\sim10^{5}\,\mathrm{M}_{\odot}$ are initially needed to influence a sufficiently large gas mass, comparable to or larger than its current mass, and build up a high-density region as shown in Figure 4. As long as the gas reservoir can be replenished from the cosmic environment, the heavy seed will be able to accrete and grow significantly in mass. In a previous study, we suggested that such replenishment is possible for the biased halos, where DCBHs originate ([70]).

As massive quasars and high-redshift AGN are relatively rare, it would not be surprising that they originate from heavy BH seeds that also require rare conditions for formation. To fully characterize the number density of heavy BH seeds, large volume simulations or semi-analytic models are needed (e.g., [50, 72, 122, 15]), enabling the comparison with the observed SMBH population statistics ([143, 42]). Previous work has indicated that rare and biased halos are able to support efficient assembly of SMBHs at high-redshifts ([38, 40, 85, 69, 82, 84]), which implies that clustering measurements of SMBH could provide further clues to their origin.

Little Red Dots

Regarding LRDs, systems proposed to be AGN with subdominant stellar populations, the spectra derived from the cutout runs, only harboring an AGN component, could explain the red part of the LRD spectrum, as shown in Figure 8. The blue part cannot be fully explained without a stellar population, which we have left out in our cutout runs.

The dense gas region triggered by the massive BH produces the Balmer breaks in the spectrum at $\sim0.35\,\mu\mathrm{m}$ comparable to the observed LRDs. However, for the most extreme LRDs like the BH* object ([110]), their extreme Balmer breaks cannot be reproduced. Even higher gas densities than those achieved in the cutouts may be needed. Higher resolution or a more massive SMBH could accomplish this, as higher resolution in the cutout, up to $0.02\,\mathrm{pc}$, led to higher gas densities near the BH, up to $10^{8}\,\mathrm{cm}^{-3}$ (Figure 6). In addition, the dense gas regions, proposed to be the source of the strong Balmer break in LRDs, are thought to have sizes of $\sim1000$ AU ([62, 90]), corresponding to $5 \times10^{-3}$ pc. Thus, even in our cutouts, we have not yet reached the scale of the dense-gas regions in LRDs (see Figure 1). In future work with even higher resolution, we may be able to directly capture the expected scale of the dense LRD gaseous cocoon, thus reproducing the most extreme cases.

We note that because of stochastic BH accretion, corresponding to intermittent periods when the BH is luminous with high $\dot{M}_{\mathrm{BH}}$, the LRD phase may also be variable ([92]). However, if the SMBH is enshrouded in a dense gas envelope as some LRD models suggest, full radiative transfer calculations will have to be performed to measure how long the gas envelope remains bright, even after BH accretion ends, and to determine its variable nature.

In any case, we find that similar to SMBH growth, heavy BH seeds are necessary to induce the dense gas envelope, and to produce LRD-like spectra for a central AGN engine. Other works have suggested that SMSs or early DCBH phases soon after formation could be the origin of LRDs, as they naturally produce a system with a gas envelope surrounding an accreting compact central object ([25, 89, 151, 115]). We complement such studies to show that to sustain an LRD-like phase after BH seed formation, a heavy BH seed is required. Light BH seeds, on the other hand, corresponding to stellar remnants, will be unable to grow in mass, as indicated by our cutout simulations, and thus cannot form a dense gas region or LRD-like systems.

Comparison with Previous works

We compare our results to select studies with a similar focus on early BH growth in a cosmological context, without aiming for completeness, to gauge the results’ robustness in terms of numerical methodology and different assumptions for the sub-grid physics.

Given that they also use the gizmo code, it is instructive to consider the FORGE’d in FIRE simulations ([59 c,a]). This simulation suite resolves the SMBH accretion disk with even higher (AU-scale) resolution within a larger cosmological simulation, albeit at high computational cost. We generally find consistent results with FORGE’d in FIRE. They resolve a region around a $\sim10^{7}\ M_{\odot}$ SMBH, and find a similar build-up of gas density and temperature ($\sim10^{8}\ \mathrm{cm}^{-3}$ and $\sim10^{3}$ K at $\sim0.1$ pc; [59]). We note that FORGE’d in FIRE finds that SF occurs down to $\sim10^{-2}$ pc away from the BH ([59 a]). However, they further measure SF suppression within $10 \sim100$ pc, inside the region where the BH gravitational potential dominates. With very low expected SF efficiencies, the possibly unresolved SF at the smallest scales should not significantly impact BH accretion. In addition, FORGE’d in FIRE does not find a significant impact from magnetic fields at larger scales. At sub-pc scales, magnetic fields stabilize the accretion disk against fragmentation and enhance the gas inflow near the BH, but these effects are present at $\lesssim10^{-1} - 10^{-2}$ pc scales ([59 c]), below our resolution limit. The absence of magnetic fields in our cutout simulations thus should not affect our conclusions for the dense gas buildup and BH mass evolution.

The BRAHMA simulations ([14]) found that hierarchical mergers of lower mass ($\sim10^{3}\ M_{\odot}$) seeds could produce BHs of $\sim10^{5}\ M_{\odot}$, providing an alternative pathway for heavy seed formation. However, they find that the $10^{5}\ M_{\odot}$ BHs also grow mainly by mergers until $\sim10^{6}\ M_{\odot}$, then through gas accretion to higher masses ([15, 17]). The subdominant growth through accretion may arise from unresolved dense gas near the massive seeds due to the lower resolution in cosmological simulations (Figure 1), and the Bondi-Hoyle formalism underestimating accretion rates. Still, they reach qualitatively similar conclusions - that a massive BH is initially required to enable efficient accretion and produce the observed high-redshift SMBHs and quasars.

We note that other recent studies have found contrary results. The SEEDZ simulations ([109]), employing extremely high resolution (in a separate, smaller volume run than shown in Figure 1), reach $\sim0.1$ pc at the smallest scales, similar to this work. They conclude that light stellar-remnant seeds with initial masses of $\sim100\ M_{\odot}$ could rapidly grow to $\sim10^{3} - 10^{4}\ M_{\odot}$ in bursts of super-Eddington accretion, even without mergers. Below, we consider possible reasons for the conflicting results.

SEEDZ does not include radiative (non-SN) feedback from stellar particles before they collapse into a light BH seed. This choice is based on previous work finding that the dense gas around Pop III stars traps the ionizing radiation, so that stellar feedback would not impact the subsequent evolution of the larger-scale environment ([125, 68]). However, these studies considered relatively limited time evolution ($\sim20$ kyr), or had to post-process the ionizing radiation due to high computational costs. Other studies that include radiation-hydrodynamical effects have found contrasting results, where the (initially D-type) ionization front surrounding massive Pop III stars eventually breaks out (becoming R-type) to reach beyond the host halo, thus severely limiting the accretion onto the BH remnants [150, 3, 1, 74]. Even if the progenitor star is initially located inside a dense gas region, the R-type ionization front after break-out (at $\sim60$ kyr for a $\sim100\,M_\odot$ Pop III star) enables a champagne flow that significantly drives down the density. This low-density environment also makes the subsequent SN feedback more disruptive (e.g., [47, 127]). In contrast to SEEDZ, our simulations do include photoionization feedback from stellar particles before BH formation, which will significantly hinder subsequent BH growth [77, 71]. Therefore, the rapid stellar remnant mass growth may arise from the optimistic assumption of not including stellar radiative feedback.

Another possibility is that the difference arises from the short timescales and limited number of cases simulated in our work (see § 4.4). The cutouts can only be run for short periods ($< 10$ Myr) due to the artificial vacuum boundary conditions. Given that the SEEDZ simulation finds bursts of super-Eddington accretion that are separated by $\sim50 - 100$ Myr, our cutouts could represent a dormant period, and we might not have caught the rare conditions enabling a high-accretion phase. The SEEDZ study further concludes that only a small fraction of stellar remnants experience rapid growth (of order a few percent), whereas we only simulate a small number of light-seed cases. Our choice of a light seed could be one of the common dormant cases. Further work is needed to more fully explore this statistics.

Caveats

A number of caveats needs to be considered in this work. The cutout simulations, by their nature, can only be run for short periods of time. The larger cosmic environment will eventually affect the cutout region, after which the cutouts will be physically unrealistic. Thus, our analysis have been limited to $\lesssim8$ Myr. While we have examined the BH evolution over short periods, this behavior may change over cosmological timescales. One example concerns the gas reservoir replenishment. Figure 7 shows short periods of super-Eddington accretion leading to BH growth, but for such behavior to continue, the accreted gas will have to be continuously replaced. More substantial gas replenishment will be needed as the BH mass grows to maintain super-Eddington bursts. Therefore, the BH evolution predicted in this work might not be long-lasting.

Furthermore, we have only examined one case of DCBH and stellar remnant BH in this work. We have chosen the most optimistic cases for BH growth based on gas environment (Sec. 2.1.1), but a diversity of BH evolution pathways may exist, reflecting environmental factors beyond the BH formation mechanism. However, both the limited time coverage and number of cases probed, do not change our conclusions that to produce the observed high-redshift AGN and LRDs, heavy BH seeds may be favored. How the gas reacts to BH mass/gravitational potential should remain consistent over long periods, and if the light BH seeds are unable to affect the nearby gas initially to grow in mass, they will not be able to grow through accretion over cosmological timescales either. Moreover, previous works have found that heavy DCBH seeds usually form in massive and biased halos, whose gas environment support efficient accretion to be able to produce the observed high-redshift SMBHs and overmassive systems [30, 70, 28], so the short-term efficient mass growth we see for heavy seeds in cutouts could be sustained for longer periods, producing the observed high-redshift SMBHs.

In addition, we are not able to fully encapsulate BH physics even with the cutouts. Our sub-parsec resolution still does not reach the inner accretion disk scales, as shown in Figure 1, such that the approximate Bondi-Hoyle formalism has to be used. While artificial boost factors are not needed, unlike in cosmological simulations with much lower resolution that cannot resolve the dense gas buildup near the BH, the BH accretion modeled in this work is still an approximation. The actual BH growth at the event horizon may be lower or higher than the growth inferred in this work. General relativistic magnetohydrodynamics simulations that directly model the accretion disk have found that the Bondi-Hoyle model actually overestimates BH accretion [26] (K.-Y. Su et al. 2026). Another key aspect of BH physics we cannot fully model is the resulting feedback. As stated above (Sec. 3.2), we currently do not accurately model the expansion of the (initial) compact thermal feedback bubble into the anisotropic surrounding medium. The feedback model built for lower resolution cosmological simulations imparts BH feedback energy to all particles within the BH kernel, while the realistic interaction between the (now resolved) dense anisotropic medium and the expanding BH feedback bubble will be more complex. Processes such as photon trapping and propagation of the feedback front are not modeled. Thus, the gas response in the cutouts to BH feedback may be overestimated. Future work will be needed to reach even higher resolution, explore a diversity of BH seed environments, and resolve this issue of consistently modeling BH feedback.

Conclusions

To address the origin of high-redshift SMBHs and LRDs, we have run targeted hydrodynamical simulations within cosmological initial conditions. Specifically, we have used the cutout technique, where we excise a small region of the gas environment and resimulate it at a higher resolution for short periods to connect the large-scale cosmological environment to small-scale BH evolution. We explore two BH evolution scenarios, a DCBH heavy seed, and a stellar remnant BH light seed. The massive heavy seed with its deeper gravitational potential well is able to influence the nearby gas and build up a dense gas region around it. Such dense gas around the BH enables intermittent super-Eddington accretion, and may produce LRD-like spectra. The light seed, on the other hand, given its shallow potential well, cannot significantly affect the nearby gas. Consequently, there is no dense gas buildup around the BH, resulting in weak accretion and mass growth during the cutout run. Light seeds are thus expected to stay dormant in general, if they only grow through gas accretion.

Our results suggest that rare heavy BH seeds are necessary and most likely the progenitors of high-redshift SMBHs and LRDs, whether they originate from a DCBH or a different heavy BH seed scenario. Follow-up work will be needed to further test this proposal, as our current study is exploratory and only considers individual cases. Other simulations have found contrary results where light BH seeds could grow significantly through gas accretion (up to $\sim10^{4}\ \mathrm{M}_{\odot}$; [109]), although still not to the heavy DCBH seed regime ($\sim10^{5}\ \mathrm{M}_{\odot}$), and invoking idealized assumptions such as neglecting stellar radiative (non-SN) feedback. As different simulations adopt different formalisms and subgrid models for BH evolution, future investigations will be needed to determine what causes such differences, whether they arise from true physical effects or numerical variations. The number densities and formation times of heavy BH seeds are also not well constrained, due to the uncertainty in their formation criteria ([141, 64, 27, 28]). Large volume simulations testing the different heavy BH formation scenarios will be necessary to compare the number of heavy seeds to the observed high-redshift SMBHs and LRDs, thus assessing whether sufficient heavy seeds could be produced to match the observations.

In the future, to further elucidate the emergence of high-redshift SMBHs, observations at higher redshifts or reaching lower BH masses will be needed ([50, 117, 112, 152, 148, 72]). Previous studies have shown that by the time they are observed, SMBHs have lost their seeding signature ([143, 72]). To discriminate SMBH origins, gravitational lensing can enable observations of $\sim10^{5}\ \mathrm{M}_{\odot}$ BHs at the predicted heavy seeding scale for DCBHs ([42]), thus providing direct confirmation of the DCBH scenario in the near-future. Upcoming wide-field surveys with the Roman Space Telescope, predicted to discover AGN at $z \sim6 - 10$ ([155]), will be able to perform population level studies of high-redshift SMBHs and constrain the various theories of BH formation. Farther ahead, multi-messenger searches with gravitational wave observatories will constrain SMBH populations in novel and complementary ways (e.g., [4]). The LISA observatory will detect binary systems to measure SMBH populations and the contribution of mergers to their formation ([129]), a channel not explored in this work.

For LRDs, searches are now being extended to lower redshifts ([102, 81]). As LRD observations have been focused at high-redshift, constraining the evolution of their physical properties and population number densities at later times will be crucial (D. D. Kocveski et al. 2023; Y. Ma et al. 2026). Identifying the LRD population at low redshifts will also provide strong tests for the heavy seed scenario for LRDs ([73, 22, 115]), which have been invoked in part to address their apparent high-redshift abundance.

Overall, the discoveries of AGN and LRDs at high-redshift enabled by JWST, and in close synergy with theoretical modeling, have invigorated questions on the formation of the first SMBHs. These remarkable beacons of the early Universe may well hold the key to unlocking the physics of structure formation at Cosmic Dawn.

Appendix

Numerical test of the cutout method

To confirm the validity of the cutout method, we ran a simulation of an idealized gas sphere. The sphere is initialized with a uniform gas density and total mass of $1.5 \times10^{8}\ M_{\odot}$, radius 1,000 pc, and gas particle mass $150\ M_{\odot}$. We evolve the gas sphere without SF effects, similar to the cutout runs, thus focusing on the gas collapse.

We run the gas sphere for $\sim$ 13 Myr and apply a similar method as described in § 2.1 to prepare a cutout of a 100 pc region inside the gas sphere. The particles within 100 pc of the sphere center are split by a factor of 100, resulting in 1.5 $M_\odot$ particles. Both the cutout and the original sphere are then run until $\sim$ 25.5 Myr to compare the gas collapse for the two resolutions. We find that, as expected, the cutout is able to resolve the gas collapse to smaller scales than the original sphere, as shown in Figure 9. The gas expands at the outer edge of the simulations due to the vacuum boundary conditions, but the collapse closer to the center of the sphere is better resolved in the cutout run. We further fit a functional form of the analytical Larson-Penston model for gravitational collapse of an isothermal sphere ([93, 119]) at two different times, finding good agreement in the inner regions, with the expected behavior of increasing central densities over time. The cutout approach can thus be used to better resolve the small-scale hydrodynamic physics at a limited computational cost, offering an efficient technique to the broader numerical toolkit.