Overmassive black holes and little red dots naturally form in simulations

Our cosmological initial conditions are based on the simulation Phi-4096 of ref. ⁵¹. In this work, a dark-matter-only N-body simulation was performed in a comoving box of side length 16 h⁻¹ Mpc. The initial condition is generated by MUSIC⁵² at a redshift of 127. The cosmological parameters used follow the latest measurement by Planck⁵³: Ω_m = 0.31, Ω_Λ = 0.69, Ω_b = 0.048, h = 0.68, n_s = 0.96 and σ₈ = 0.83. The simulation used 4,096³ dark-matter particles, corresponding to a particle mass of 5.13 × 10³ h⁻¹ M_⊙. This high resolution allows the identification of mini-halos with masses down to 10⁵ h⁻¹ M_⊙, sufficient to resolve the sites of Population III (hereafter Pop III) star formation and to track chemical enrichment across cosmic time. Halo merger trees were constructed from simulation snapshots between z= 35 and z = 7.5 using the ROCKSTAR phase space halo/subhalo finder⁵⁴ and the CONSISTENT-TREES merger tree code⁵⁵. On top of these, a semi-analytic galaxy-formation model was implemented to follow gas cooling, Pop II and Pop III star formation, supernova feedback, metal enrichment and local Lyman–Werner (LW) radiation fields¹⁷. We use the result of their model with no baryon streaming motion (0σ_vbc).

The details of the semi-analytic model largely follow those described in ref. ⁵⁶ with recent updates of the treatment of early star formation⁵⁷. Star formation in Pop II halos follows the standard prescription used in galaxy-formation models, in which the cold-gas component is converted into stars on a timescale t_SF = t_dyn/α_* with efficiency α_* = 0.03 (ref. ⁵⁸). Along with star formation, the model assumes that the rate at which cold gas is reheated into the hot-gas phase by associated supernova feedback is proportional to the star-formation rate, (gamma {dot{M}}_{* }), in which γ = (V_halo/110 km s⁻¹)^−1.74 (refs. ^56,59). When the progenitor halos never formed any stars, the stellar population follows a Pop III initial stellar mass function (IMF). The characteristic Pop III stellar mass is determined by the halo mass growth rate, which correlates with the final stellar mass⁶⁰. Stellar populations emit LW radiation that photodissociates molecular hydrogen and delays star formation in nearby halos. The LW radiation intensity is calculated separately for Pop II and Pop III stars⁶¹:

in which r_i is the distance to halo i and M_PopIII,i and M_PopII,i are the masses of Pop III and Pop II stars formed within the past 5 Myr in halo i, respectively. The sums run over all halos that host Pop III or Pop II star formation. The coefficient in the Pop II expression is derived assuming a Scalo IMF (ref. ⁶²) and a stellar metallicity of Z = 0.001 (refs. ^63,64).

The LW background is primarily contributed by Pop II stars and its amplitude depends on assumptions about the IMF and the stellar models used. We note that uncertainties in the LW intensity modelling have only a small impact on our results, as we later test explicitly. A halo begins to form Pop III stars once its mass exceeds the critical threshold determined by the local LW intensity⁶⁵. Under strong LW irradiation, star formation is delayed until the halo virial temperature reaches about 8,000 K, at which point Lyα cooling triggers rapid collapse. Such halos typically experience high mass accretion rates and host more massive Pop III stars, extending to supermassive stars with M_* ≈ 10⁵ M_⊙, which subsequently collapse into heavy BH seeds. Indeed, Ishiyama and Hirano¹⁷ identified more than about 10⁴ supermassive stars with masses larger than 10⁵ M_⊙ as heavy-seed candidates in their simulation volume.

From this dataset, we select one representative halo predicted to form a massive seed, which serves as the target of our follow-up radiation-hydrodynamic simulations described in this work. The halo is chosen according to two criteria: (1) the local LW intensity at the time of seed formation exceeds the critical value of J_21,crit = 1,000 and (2) the halo is not tidally disrupted by nearby massive galaxies. Here J₂₁ denotes the FUV intensity normalized to 10⁻²¹ erg s⁻¹ Hz⁻¹ cm⁻² sr⁻¹. The second condition is necessary because candidate halos are often located near luminous galaxies that provide strong LW irradiation, but their gravitational collapse may be suppressed by external tidal fields⁶⁶. To exclude such cases, we select halos whose distances from the nearest massive galaxy satisfy r_dist > 10d_tidal, in which d_tidal is defined as the separation at which the tidal radius imposed by the neighbouring galaxy equals the virial radius of the candidate halo. Within the simulated volume, we identify 60 halos that meet both criteria. Among them, we focus on the one that exhibits the most rapid mass growth, reaching a virial temperature of T_vir ≃ 8 × 10³ K, which marks the onset of atomic cooling and subsequent collapse.

We trace back the initial position of the selected candidate halo to the cosmological initial conditions at redshift z = 127. Extended Data Fig. 1 shows the time evolution of the halo mass, in which MBH1 forms. The white star marks the redshift at which the virial temperature of the halo reaches 8,000 K, when it is identified as a DCBH halo. A zoom-in region is defined to be 40 times larger than the Lagrangian radius of the target halo, corresponding to a comoving size of about 400 kpc. We confirm that this region is sufficiently large to encompass all material relevant to the formation of the heavy-seed BH in the selected halo. We measured the overdensity of the specified zoom-in region by computing the mean density within its bounding box of about 400 comoving kpc and found it to correspond to a 3.8σ fluctuation.

Radiation-hydrodynamic calculation

The radiation-hydrodynamic calculation is performed within this zoom-in region using the moving-mesh code AREPO (ref. ¹⁶), extended to include prescriptions for star formation and BH accretion physics. The resolution of the zoom-in region at the initial snapshot is identical to that used in ref. ⁵¹, in which the dark-matter particle and baryonic cell masses are 4.33 × 10³ and 7.94 × 10² h⁻¹ M_⊙, respectively. At this resolution, star formation within mini-halos can be reliably resolved. Adaptive mesh refinement is applied whenever the local cell size falls below 16 times the local Jeans length, ensuring that gravitational collapse is properly captured and that artificial fragmentation is avoided⁶⁷. Unlike the semi-analytic model in ref. ¹⁷, the formation and properties of primordial stars are directly followed by the high-resolution radiation-hydrodynamic simulations described here. Uncertainties in the stellar mass and the resulting FUV background intensity in the semi-analytic model may affect whether a heavy-seed BH forms or not, but we later show that these uncertainties do not substantially alter our main conclusion about the rapid emergence of overmassive BHs.

We solve a non-equilibrium primordial chemical network consisting of eight species: e⁻, H, H⁺, H₂, H⁻, D, D⁺ and HD. The network and associated cooling processes follow the implementation in ref. ⁶⁸ and include H₂ rovibrational line cooling, atomic hydrogen line cooling (including Lyα) and free–free and free–bound emission of H and H⁻. The ionization of H, the dissociation of H₂ by LW radiation through the Solomon process and the photodetachment of H⁻ are also taken into account. The chemical network and reaction rates used are similar to those in the minimal chemistry model presented in ref. ⁶⁹, except that we neglect ({{rm{H}}}_{2}^{+}) and He. ({{rm{H}}}_{2}^{+}) can catalyse H₂ formation in primordial gas but it is less important in environments exposed to the background radiation field assumed in this study, namely a 10⁵ K blackbody spectrum⁶⁹. Heating and chemical reactions induced by FUV, extreme ultraviolet and X-ray radiation from nearby stars and BHs are also included, as described later.

Helium chemistry becomes important at gas temperatures of several 10⁴–10⁵ K, at which it provides further cooling channels and increases the free-electron fraction. Although we do not solve the full non-equilibrium helium chemistry, we estimate the ionization state of helium by assuming chemical equilibrium and include the line-cooling rates of He I and He II. This effect is relevant for gas irradiated by both Pop III stars and accreting BHs. Radiation from the accreting BH raises the gas temperature well above 10⁴ K. The elevated temperature already increases the electron fraction and thereby promotes H₂ formation; including non-equilibrium helium chemistry would further enhance the H₂ abundance to some extent. However, this does not affect our main conclusion. Even if the enhanced H₂ abundance lowers the cloud temperature, the large-scale inflow can be sustained. This is demonstrated by the formation of MBH2, in which a heavy-seed BH forms under a massive inflow despite efficient H₂ cooling that lowers the gas temperature below 1,000 K (Extended Data Fig. 7).

To avoid numerical contamination from the coarse outer region, gas cells outside the zoom-in volume are kept at fixed resolution and both radiative cooling and refinement are disabled. In this calculation, we treat the FUV radiation from the source halos separately from the radiation emitted by stars and BHs formed in the target region. Here we define source halos as halos that have already been labelled as star-forming when the target halo satisfies the DCBH formation criteria in the semi-analytic model. These halos host star-forming galaxies and provide the strong FUV radiation required for DCBH formation. Instead of explicitly following the star formation within the galaxies, we use a time-dependent LW intensity derived from the semi-analytic model in ref. ¹⁷, which self-consistently evolves the star formation and feedback in the same cosmological volume. The resulting LW flux is applied uniformly to all gas cells within the zoom-in region as a uniform background field. We separately solve the ionizing radiation and X-rays emitted by stars and BHs formed inside the zoom-in region, but outside the source halos, using the ray-tracing scheme implemented in regularized smoothed particle hydrodynamic^70,71. We include photoionization of H and the associated photoheating, photodissociation of H₂ and HD and photodetachment of H⁻. The photodetachment rate of H⁻ is calculated under the optically thin approximation, whereas the other photoreaction rates are calculated using regularized smoothed particle hydrodynamics. We do not include the secondary ionization and heating. The luminosities and spectra of these sources are described later in the section ‘Modelling Pop III stars and BHs’.

Extended Data Fig. 2 shows the time evolution of the LW intensity contributed by nearby source galaxies. The intensity increases rapidly and exceeds J₂₁ = 1,000 around redshift z ≃ 18, the critical level required to suppress H₂ cooling and enable heavy-seed BH formation (for example, refs. ^72,73,74,75). The spectrum of the background radiation is modelled as a blackbody with an effective temperature of T_BB = 10⁵ K, which reproduces the relative rates of H₂ photodissociation and H⁻ photodetachment expected from young star-forming galaxies⁷³. We consider the self-shielding of the background LW radiation using the prescription described in ref. ⁷⁶. We assume that the cloud extends over the local Jeans length, l_Jeans, and estimate the H₂ column density as ({N}_{{{rm{H}}}_{2}}={n}_{{{rm{H}}}_{2}}{l}_{{rm{Jeans}}}), in which ({n}_{{{rm{H}}}_{2}}) is the number density of H₂ in the cell. We neglect the contribution of ionizing photons from neighbouring galaxies because their mean free path is short and they are strongly attenuated by the surrounding intergalactic medium⁷¹.

During the radiation-hydrodynamic simulation, we do not assume any specific IMF but instead resolve individual star-formation events directly using a sink-particle method. A sink particle is introduced once the local gas density exceeds 2 × 10⁶ cm⁻³, with an accretion radius set to ten times the local mesh size. The sink accretes surrounding gas following the local gravitational potential and its properties are updated accordingly at each time step. We allow mergers of the sink particles once the distance of two sink particles becomes smaller than the sum of the sink radii.

Modelling Pop III stars and BHs

The luminosity and effective temperature of each protostar are determined from its instantaneous mass and accretion rate by interpolating results from one-dimensional stellar evolution models assuming constant accretion rates⁷⁷. If the stellar accretion rate exceeds the critical accretion rate of ({dot{M}}_{{rm{c}}{rm{r}}{rm{i}}{rm{t}}}equiv 0.02,{M}_{odot },{{rm{y}}{rm{r}}}^{-1}), we assume that the stellar radius expands owing to the injection of the large amount of entropy into the stellar envelope^28,78,79,80. During this phase, we assume the surface temperature of the star to be 6,000 K and the luminosity to be the Eddington value following the results of detailed stellar evolution calculations⁸¹. After the accretion rate becomes smaller than ({dot{M}}_{{rm{crit}}}), we gradually shrink the stellar radius to the radius expected for main-sequence stars on the timescale of the surface Kelvin–Helmholtz time, which is given by ten times the stellar Kelvin–Helmholtz time⁸². Once the stellar age exceeds the lifetime predicted by these models, the star is converted into a BH particle if its final mass exceeds 260 M_⊙ (ref. ⁸³). We note here that all of the stars formed during the simulation (the progenitor stars of LBH, MBH1 and MBH2) have masses greater than 260 M_⊙ and collapse into the BHs without any supernova feedback. The mass accretion rate is measured by the sink method, in which the gas inside the sink radius is assimilated and added to the mass of the star particle. The sink radius is set to be 100 times the cell size, which will be converted into the sink particle, which ranges from 1 to 20 pc. Because the sink radius is much larger than the physical stellar radius, we estimate the final stellar mass from the sink accretion rate using the empirical relation derived from previous high-resolution simulations of Pop III star formation^60,84,

in which ({dot{M}}_{ast }) denotes the time-averaged accretion rate onto the growing protostar. When the stellar age reaches its lifetime, any excess gas mass accreted by the sink that exceeds the estimated stellar mass is returned to the surrounding gas cells to conserve mass. Throughout the simulation, we assign the sink-particle mass as the instantaneous stellar mass when evaluating radiative feedback. This approach overestimates the stellar luminosity and gives stronger feedback effects compared with the actual stellar mass, so the final stellar mass represents a lower limit to the actual value.

We assume a stellar lifetime of 2 Myr, typical for very massive stars with masses greater than 100 M_⊙ (ref. ²⁸). After the stellar age exceeds 2 Myr, we convert the stars into BHs. However, the lifetime of such stars—particularly supermassive stars—remains uncertain. Stars with masses greater than several 10⁵ M_⊙ become unstable owing to general-relativistic effects and collapse into BHs during the hydrogen-burning phase^{28,85,86,87,88,89}. This implies that collapse may occur earlier than the canonical lifetime assumed for massive stars. Recent calculations⁹⁰ show that the threshold mass for general relativistic instability depends on the mass accretion rate, spanning 2–8 × 10⁵ M_⊙. If the stars collapse into BHs earlier, the resulting BHs would remain embedded in dense gas for a longer period, thereby extending the phase of super-Eddington growth.

After the stars collapse into BHs, we convert the corresponding sink particles into BH particles (LBH, MBH1 and MBH2). We keep their original accretion radii if the accretion radius is larger than the Bondi radius of the BHs for the ionized gas, otherwise the accretion radii are set to be the Bondi radii. Because the typical accretion radius is around a parsec, the simulations resolve the Bondi radius for most of the massive-seed BHs, allowing us to directly follow the gravitational capture of gas onto the BHs. The spectral energy distribution of accreting BHs is modelled as a double power-law continuum, F_ν ∝ ν^−0.6, extending from 1 eV to 10 eV, and F_ν ∝ ν^−1.5, extending from 10 eV to 1 keV (ref. ⁹¹). This represents the integrated emission from the accretion disk and its associated corona. We allow the accretion rate to exceed the classical Eddington limit, consistent with theoretical models of super-Eddington accretion flows (for example, refs. ^23,25,92,93). We assume the slim disk solution to give a bolometric luminosity when the accretion rate is higher than the Eddington rate,

in which ϵ_r = 0.1 is the radiative efficiency, ({dot{M}}_{{rm{acc}}}) is the instantaneous gas accretion rate onto the BH, ({dot{M}}_{{rm{Edd}}}) is the Eddington accretion rate and ({L}_{{rm{Edd}}}equiv {dot{M}}_{{rm{Edd}}}/{{epsilon }}_{{rm{r}}}{c}^{2}) is the Eddington luminosity⁹⁴. The emitted radiation is coupled to the surrounding gas through photoionization heating, which regulates the accretion flow and drives intermittent outflows around the BHs. We also allow mergers of BHs once the distance between two BHs becomes smaller than ten physical parsec.

We do not include kinetic feedback from accreting BHs, such as jets or winds. Recent general relativistic magnetohydrodynamic simulations have shown that BHs accreting at super-Eddington rates can launch magnetically arrested jets and disk winds. If such jets break out and transport mass and momentum to larger scales in the surrounding medium, the accretion efficiency may decrease to less than the Eddington value^95,96. The efficiency of jet launching should depend on the BH spin and magnetic field strength, both of which are highly uncertain. This uncertainty arises because modelling them would require following the stellar rotation at the birth of DCBHs and the detailed magnetic dynamo processes during accretion. The jet opening angle and its interaction with the accreting gas are also uncertain. We therefore consider our simulation as a fiducial model that provides an upper limit on BH growth under the assumption that kinetic feedback is absent.

Zoom-in calculation of seed BH formation in the isolated cloud

To assess the feasibility of heavy-seed BH formation, we perform a higher-resolution follow-up simulation that resolves the individual protostars expected to be the progenitors of the seed BHs. We focus on the formation site of the first heavy-seed BH, MBH1, identified in the cosmological radiation-hydrodynamic run. This seed, with a final mass of 6 × 10⁵ M_⊙, forms at redshift z ≃ 14. We extract a cubic region of one comoving kpc on a side centred on the host halo, shortly before the onset of collapse when the central gas density reaches 10³ cm⁻³. At this point, the central region of the cloud, of size a few tens of pc, becomes Jeans-unstable (Extended Data Fig. 6). The re-simulation is performed using the same AREPO framework but with enhanced spatial and mass resolution to capture the internal fragmentation of the collapsing gas cloud. Sink particles are introduced once the local gas density exceeds 2 × 10⁸ cm⁻³, representing the formation of individual protostars. The sink radius is set to be four times the cell size, which will be converted into the sink particle, which ranges from 3,000 to 8,000 au. We allow mergers of the sink particles once the distance of two sink particles becomes smaller than the sum of the sink radii.

The numerical resolution used here is insufficient to resolve the opacity limit at n ≈ 10¹⁶ cm⁻³, above which the gas becomes optically thick and protostars form⁹⁷. Indeed, the physical sizes of growing protostars are expected to be about 0.05–100 au, depending on the stellar mass and accretion rate⁸¹, which is below our numerical resolution. We can, nevertheless, compare our results with higher-resolution studies. Becerra et al.⁹⁸ studied the protostellar evolution changing n_adib, the density above which the gas becomes adiabatic, from 10⁸ cm⁻³ to 10¹² cm⁻³. They have found that the mass evolution of supermassive protostars is not greatly affected by the numerical resolution, because the growth rate is mainly regulated by the large-scale accretion flow. They followed the early stellar evolution for 10⁴–10⁵ years and found that the accretion rate remains at about 1 M_⊙ yr⁻¹, comparable to our results. Similarly, Chon and Omukai⁹⁹ followed protostellar evolution in a DCBH halo for 2 Myr with n_adib = 10¹¹ cm⁻³. They found that the initial mass accretion rate remains around 1 M_⊙ yr⁻¹. In their simulation, however, the accretion rate drops below 0.01 M_⊙ yr⁻¹ at 10⁴–10⁵ years after protostar formation because of the limited gas supply. This contrasts our results, in which efficient accretion continues until around 2 Myr after protostar formation.

Zoom-in calculation of the BH accretion phase

To reproduce the density structure around the massive-seed BH, we performed a high-resolution simulation that resolves gas densities up to 10¹² cm⁻³ and spatial scales down to 500 au around the central BH. We take the snapshot at the moment when MBH2 forms at z = 13.1 and follow the evolution for 30 kyr after its emergence. To capture the later evolution up to 0.5 Myr after MBH2 formation, we also run a complementary simulation with a larger sink radius of 5,000 au. Starting from the snapshot at 0.48 Myr, we then reduce the sink radius back to 500 au and continue the calculation to follow the detailed structure of the circum-BH disk for a further 0.55 kyr. This procedure allows us to track both the long-term disk evolution and the small-scale accretion flow onto the BH.

Host halo evolution of massive BHs

Extended Data Fig. 3 shows the redshift evolution of the distance between the source galaxy and the host halos of MBH1, MBH2 and LBH. All BHs in our simulation form at distances of about ten physical kpc from the source galaxy, outside the virial radius of the source halo (dashed line). After their formation, the BHs migrate towards the source galaxy and enter the virial radius of the source halo at z ≈ 12. Initially, the BHs follow eccentric orbits, as indicated by the oscillatory behaviour of their distances. Over time, however, their orbits gradually decay and the BHs settle towards the galaxy centre by around z ≈ 8. We also find that, around this epoch, their accretion rates increase to about 0.1–1 times the Eddington value (Fig. 1f).

The host halos of the BHs may avoid metal enrichment because they are well separated, by about ten physical kpc, from the source galaxy that drives winds and pollutes the surrounding intergalactic medium with metals. Analytic estimates of galactic-wind expansion in ref. ¹⁰⁰ show that such winds expand only to around 1–2 physical kpc within 300 Myr. Similarly, Ventura et al.¹⁰¹ statistically studied the distribution of superbubbles around star-forming halos in a 10 h⁻¹ Mpc cosmological region and found that the bubble size is at most about 5 kpc. These scales are much smaller than the separation found in our simulation, suggesting that the BH host halos can remain metal-poor despite their proximity to the source galaxy.

We have shown that the key to forming an overmassive BH population is the formation of DCBHs in massive halos, with virial temperatures reaching about 40,000 K. We attribute this to the delayed onset of protostar formation after Lyα cooling becomes effective, caused by the finite collapse timescale of the halo and by dynamical heating. Extended Data Fig. 4 shows the time evolution of the halo growth rate and the peak gas density in the host halo of MBH1. The gas density increases steadily but protostar formation (at 300 Myr) occurs only around 90 Myr after the halo virial temperature first reaches 8,000 K (at 210 Myr), which is often assumed to mark the onset of DCBH formation. This delay arises for two reasons. First, the collapse timescale is on the order of the free-fall time,

in which n is the gas density of the gravitationally unstable region. This implies that DCBH formation should occur only after a timescale on the order of approximately 10 Myr once the halo virial temperature exceeds 8,000 K. The halo mass can grow during the initial free-fall phase to exceed T_vir = 8,000 K. Second, dynamical heating caused by halo mergers further delays the collapse. During the first roughly 30 Myr after the virial temperature reaches 8,000 K, the peak gas density increases only weakly. Around this epoch, the host halo merges with nearby halos or clumps, which deposits further gravitational and kinetic energy into the collapsing gas and suppresses runaway collapse^102,103. After these clump mergers subside and the halo growth rate decreases, the peak gas density begins to rise rapidly, eventually leading to the formation of the protostar that later collapses into MBH1. During this delay, the host halo grows sufficiently massive to enable subsequent super-Eddington accretion.

Formation of extremely massive-seed BHs

Here we show how extremely massive-seed BHs form—through the formation of heavy seeds followed by a brief phase of super-Eddington accretion. Extended Data Fig. 1 shows the time evolution of the mass of the halo that later hosts the heavy-seed BH (MBH1). At the time of protostar formation, the halo mass is about 2 × 10⁸ M_⊙, corresponding to a virial temperature of 4 × 10⁴ K. This is much higher than the canonical virial temperature of 8,000 K at which halo collapse is typically expected. Indeed, the semi-analytic model predicts that this halo would have collapsed at z ≈ 18.1, when its mass was about 10⁷ M_⊙, nearly an order of magnitude smaller than the host halo of MBH1. During the delayed collapse, more gas accumulates within the halo, leading to the formation of more massive stars and enhanced accretion onto the resulting BHs.

Once the gas within a halo becomes gravitationally unstable, it collapses to form stars. The final stellar and BH masses are determined by how much gas is accreted onto the central protostar during its growth phase. Extended Data Fig. 5 shows the radial profiles of gas density (Extended Data Fig. 5a), temperature (Extended Data Fig. 5b) and escape velocity (Extended Data Fig. 5c) at the time of protostar formation. The grey line represents the light-seed case (LBH) and the purple and green lines correspond to the two massive seeds, MBH1 and MBH2, respectively. The density structure differs markedly between light and heavy seeds. In the massive-seed cases, the high-density core extends to radii an order of magnitude larger than in the light-seed case, a consequence of the higher gas temperature in the collapsing cloud. Extended Data Fig. 5b shows that the temperature in heavy-seed formation remains near 10⁴ K—an order of magnitude higher than in the light-seed case—owing to strong ultraviolet irradiation that suppresses molecular cooling. The collapse of this larger, hotter core also raises the escape velocity, which exceeds 10 km s⁻¹. This high binding energy facilitates continued mass accretion, as the ionized gas remains gravitationally bound and can feed the central protostar efficiently.

To visualize how gravitational collapse proceeds, we compare the enclosed mass with the critical Bonnor–Ebert (BE) mass, the threshold above which a cloud becomes gravitationally unstable to collapse^104,105. The top panels in Extended Data Fig. 6 show the BE mass (green) and the enclosed mass (purple) for LBH, MBH1 and MBH2, from left to right. In the light-seed case, only the innermost 0.1–1 pc region is gravitationally unstable, whereas in the heavy-seed cases, the unstable region extends to 10–100 pc. The bottom panels show the ratio of enclosed to BE mass as a function of enclosed mass; regions in which this ratio exceeds unity are gravitationally unstable. Only gas within about 10³ M_⊙ becomes unstable in the light-seed case, whereas up to about 10⁶ M_⊙, it is unstable in the heavy-seed cases. These characteristic mass scales correspond roughly to the resulting BH masses: about 800 M_⊙ for LBH and 3 × 10⁵ M_⊙ and 6 × 10⁵ M_⊙ for MBH1 and MBH2, respectively. The dynamical time at the edge of the gravitationally unstable region is

which is longer than the lifetime of massive stars (about 2 Myr). This implies that the dense envelope remains bound and should continue to fall onto the BHs after the central protostars collapse, unless feedback from the protostars clears it away.

The temperature decrease at r ≲ 1 pc in MBH2 is caused by the enhanced electron fraction and the subsequent increase in molecular hydrogen formation, similar to the Pop III star formation in a fossil H II region¹⁰⁶. Extended Data Fig. 7 shows two-dimensional histograms of temperature (left), molecular hydrogen fraction (middle) and electron fraction (right) as functions of gas density, for MBH1 (top) and MBH2 (bottom). Gas at temperature T ≈ 100–1,000 K in MBH2 shows an enhanced H₂ fraction, indicating that molecular formation and the associated cooling are responsible for the temperature drop. The electron fraction in the outer envelope is almost fully ionized but decreases towards higher densities as recombination proceeds. The increased electron fraction during collapse promotes H₂ formation through the H⁻ channel, the dominant pathway in primordial gas, which is catalysed by free electrons. The ionization source is the strong radiation from MBH1, which fully ionizes and heats the surrounding gas. Indeed, the temperature increase at r ≳ 10³ pc in MBH2 indicates that intense radiation from MBH1 photoionizes the intergalactic medium. Despite the enhanced H₂ formation, a massive-seed BH still forms in MBH2, as large-scale gravitational instability drives a strong gas inflow that overwhelms the cooling. This behaviour is analogous to the heavy-seed formation in the super-competitive accretion scenario, in which fine-structure lines and dust cooling reduce the gas temperature, whereas large-scale collapse collects a substantial amount of gas to form massive-seed BHs^99,107.

Observability of growing massive BHs

To estimate the optical depth and Hα luminosity, we restrict the analysis to gas cells located within 30° of the disk plane. We estimate the electron column density, N_e, by summing the total number of electrons contained within each spherical shell and dividing by the surface area of the shell, yielding an average column density at radius r. The Thomson optical depth is then calculated as τ = N_eσ_T, in which σ_T = 6.65 × 10⁻²⁵ cm² is the Thomson scattering cross-section. In estimating the Hα luminosity, we have assumed case B recombination. This may underestimate Hα luminosity because resonance scattering will increase the Hα emissivity at densities of n ≈ 10¹⁰–10¹¹ cm⁻³ (ref. ¹⁰⁸). The detailed radiative-transfer calculations for similar conditions indicate that the emergent Hα/Hβ is enhanced to values of about 6–10, consistent with the observed values for LRDs³².

Present Chandra constraints indicate that LRDs are generally X-ray weak. Only three LRDs have secure individual X-ray detections so far ^{4: a sample of little red dots in the JWST extragalactic legacy fields. Astrophys. J. 986, 126 (2025).” href=”http://www.nature.com/articles/s41586-026-10985-8#ref-CR47″ id=”ref-link-section-d20440416e4122″>47,109}, whereas the remaining LRD population is largely undetected on an individual basis. Stacking analyses have also yielded mostly tentative signals or non-detections, including stacks of X-ray-undetected sources^110,111,112. To assess the detectability of X-rays emitted by the central BH in our simulation, we evaluated the hydrogen column density, N_H, at two different snapshots, t = 26 and 523 kyr. Extended Data Fig. 8 shows the distribution of N_H measured along 64 different lines of sight evenly spaced around the accreting MBH2. During the early super-Eddington accretion phase at t = 26 kyr, the column density exceeds 10²⁶ cm⁻² in all directions. This implies that X-rays emitted from the BH are strongly Compton scattered and therefore undetectable from any viewing angle, consistent with the non-detections obtained from stacking Chandra data¹¹². At the later epoch, nearly Eddington phase at t = 523 kyr, the BH remains embedded in dense gas but the column density decreases to N_H = 10²⁵–10²⁶ cm⁻², lower than in the early super-Eddington phase. Even in this phase, however, the X-ray emission is still expected to be difficult to observe. This analysis indicates that the LRD-like system found in our study remains X-ray dark for at least the first approximately 0.5 Myr, in agreement with present observational constraints.

Inefficient growth of Pop III remnant BH

We find that LBH undergoes almost negligible growth during the simulation. Its mass remains close to 800 M_⊙ until z ≈ 12, followed only by modest growth thereafter. We mainly attribute this inefficient growth to the shallow gravitational potential of the host halo at the time of Pop III star formation. Extended Data Fig. 5c shows that the escape velocity is only about 2–5 km s⁻¹, smaller than the sound speed of ionized gas, which is about 10 km s⁻¹. This indicates that the halo potential is too shallow to retain ionized gas, so the gas inside the halo is rapidly expelled by photoevaporation^{40,113,114,115}. Indeed, the gas density around LBH remains ≲1 cm⁻³, with a temperature of approximately 10⁴ K (Extended Data Fig. 5a,b). Under these conditions, the Bondi accretion rate is only about 10⁻⁸–10⁻⁷ M_⊙ yr⁻¹, much lower than the Eddington accretion rate of LBH.

Around z ≈ 12, the halo potential becomes sufficiently deep for gas to begin accumulating at the halo centre. However, the inflowing gas is predominantly supplied to MBH1 rather than to LBH. MBH1 formed in a halo located close to the original host halo of LBH and the two host halos subsequently merged. After the merger, MBH1 becomes the primary accretor because of its larger mass, whereas LBH remains less efficiently fed. This is partly because LBH has a smaller mass and therefore a longer dynamical-friction timescale than the massive BHs, further suppressing its orbital decay and growth^26,40. Although stochastic interactions with nearby gas increase the mass of LBH to several thousand solar masses, its accretion rate remains much lower than those of the massive BHs.

DCBH formation in lower-z Universe

To test how ubiquitous DCBH formation and LRD-like systems are, and to examine whether this phenomenon can occur in the later Universe, we performed radiation hydrodynamics calculations for another candidate region selected from ref. ¹⁷. To identify a lower-redshift sample, we chose a candidate halo that satisfies the DCBH formation criteria at z ≈ 14, the latest formation redshift among our 62 samples.

We find eight seed BHs forming in the region around the target halo. Extended Data Fig. 9a shows the spatial distribution of the massive BHs formed in this region at z = 7.011. At this snapshot, the lowest-redshift DCBH, labelled MBH8, has already formed. Some of the BHs have migrated to the centres of massive halos, whereas other clumps around the target halo are still forming BHs. Extended Data Fig. 9b,c shows the redshift evolution of the BH mass and the accretion rate normalized by the Eddington rate, respectively. All seed BHs grow beyond 10⁶ M_⊙ after experiencing brief episodes of super-Eddington accretion. Some BHs, such as MBH1 and MBH3, repeatedly increase their accretion rates and undergo several super-Eddington phases. These results demonstrate that DCBH formation followed by later super-Eddington accretion is not confined to the earliest cosmic epochs but can also occur at later times, overlapping with the redshift range in which LRDs have been detected.

Robustness tests

To assess numerical robustness and model sensitivity, we performed one higher-resolution run (lv13) and three simulations with reduced FUV background intensities (w2, w5 and w10). In the last three runs, we kept the parameters of the semi-analytic model fixed but reduced the LW intensity from its fiducial value obtained from the model. In the higher-resolution run, we use an effective resolution of 8,192³ in the zoom-in region, corresponding to dark-matter and baryonic particle masses of 542 and 99.2 h⁻¹ M_⊙, respectively. This resolution is sufficient to resolve mini-halos in which the first stars form^102,116. The run yields results that are nearly identical to the fiducial case: one LBH forms at z ≈ 22, followed by the formation of two massive seeds. The blue line in Extended Data Fig. 10 shows the evolution of the most massive BH in the higher-resolution run, demonstrating that the mass growth is largely insensitive to the numerical resolution achieved in our default simulation.

To examine uncertainties in the semi-analytic predictions, we also ran simulations with lower FUV background intensities. In our semi-analytic model, the stellar component of the source galaxy is populated using standard galaxy-formation prescriptions, but the star-formation rate and efficiency are subject to substantial uncertainty. A lower star-formation efficiency would result in a smaller luminosity and therefore a reduced FUV intensity incident on the massive-BH-forming halo. The reduced-FUV runs also capture possible spatial variations in the radiation field, which are not explicitly modelled in our semi-analytic treatment. If heavy seeds still form under weaker FUV irradiation, this supports both the robustness of the semi-analytic predictions and the insensitivity to spatial variations.

The green, yellow and red lines in Extended Data Fig. 10 show the BH mass evolution when the background intensity is reduced by factors of two, five and ten, respectively. The dependence on the seed-formation epoch and subsequent growth is weak. As the FUV intensity decreases, the most massive BH forms slightly earlier, but the difference in formation redshift between the fiducial run and the lowest-intensity case is less than Δz ≃ 0.1. This behaviour is consistent with the expectation that weaker FUV irradiation allows more efficient H₂ cooling, leading to earlier collapse of the target halo. The later mass growth is also only mildly affected: lower FUV intensity results in a smaller final BH mass but the variation remains within a factor of three. Even in the weakest-FUV case, the final BH reaches 1.9 × 10⁷ M_⊙ by z = 7, highlighting the robustness of overmassive BH formation despite model uncertainties.