Overmassive black holes and little red dots naturally form in simulations

Our cosmological radiation-hydrodynamic simulations naturally produce massive BH seeds that grow into overmassive (≳10⁷ M_⊙) BHs by z ≈ 10. This rapid growth occurs in overdense protocluster regions exposed to intense FUV radiation from nearby star-forming galaxies, in which the radiation suppresses early star formation. We follow the collapse of one such halo using three-dimensional radiation-hydrodynamic simulations performed with the moving-mesh code AREPO (ref. ¹⁶) (Methods). The halo—identified by Ishiyama and Hirano¹⁷ as a promising heavy-seed site—is located about 10 kpc from a luminous neighbour that provides the strong FUV flux needed to form heavy seeds (Extended Data Fig. 2). As a result, the halo accumulates a large gas reservoir before collapsing at z ≃ 14 (Fig. 1a,b). Once collapse sets in, the central protostars grow rapidly, reaching 5–9 × 10⁵ M_⊙, well above the canonical value of approximately 10⁵ M_⊙ predicted by standard direct-collapse models^{18,19,20,21,22}. The unusually deep potential well of the host halo (virial temperature about 4 × 10⁴ K; Extended Data Fig. 1) enables ionized gas to remain gravitationally bound and sustain the high accretion rates needed to form such massive seeds. Several seed-forming sites emerge within the same overdense environment, indicating that heavy-seed formation can naturally occur in protocluster regions.

These massive seeds do not remain static; each heavy seed subsequently experiences a brief but intense phase of super-Eddington accretion (Fig. 1f). Immediately after collapse, a dense and optically thick envelope develops around the nascent BH, efficiently feeding the accretion disk. Radiation trapping becomes important: photons are advected inward faster than they can diffuse outward, enabling accretion rates several to a few tens of times the Eddington limit for less than a million years^23,24,25. The inflowing gas is supplied from the circum-BH disk, inside which the gas loses angular momentum through gravitational torques. This sustains rapid growth until the reservoir is exhausted, when the BHs reach several times 10⁶ M_⊙ shortly after formation. The accretion rate then gradually decreases to 0.1–1 times the Eddington limit as the gas supply dwindles, allowing the BHs to grow to ≳10⁷ M_⊙ by redshift z ≃ 10.

At z ≃ 10, the accretion rate sharply declines to less than 1% of the Eddington limit. As the BHs pass through the central region of the massive neighbouring halo, their envelope gas is stripped by the surrounding hot medium²⁶. Dynamical friction subsequently brings them to the galaxy centre, where they merge with the host system. By redshift z ≃ 8, the BHs have settled in the central kiloparsec of a galaxy with stellar mass approximately 10⁹ M_⊙, establishing a BH-to-stellar mass ratio of about 1%, an order of magnitude above the local relation. After settling, the BHs resume sub-Eddington accretion, maintaining slow but steady growth. For a brief period, the accretion rate reaches the Eddington limit, during which the inferred BH masses and accretion luminosities are consistent with those of observed high-redshift active galactic nuclei (AGN) (ref. ²⁷). Our simulation yields a lower limit of a comoving number density of about 10⁻⁴ Mpc⁻³ for massive BH-hosting galaxies, consistent with that required to explain the observed LRDs^1,3,5,6.

To demonstrate that such rapid growth is unique to heavy seeds, we also follow the evolution of BHs seeded by Population III remnants (roughly 800 M_⊙) formed earlier than the heavy seeds at redshift z ≃ 22, hereafter referred to as a light-seed black hole (LBH). Their growth remains highly inefficient. The initial gas supply is strongly suppressed by radiative feedback from the Population III stars, resulting in final masses orders of magnitude smaller than those of the massive seeds (grey line in Fig. 1e). The gas in the host halo quickly photoevaporates during seed formation because its shallow gravitational potential cannot confine the fuel required for accretion. This contrast highlights the necessity of heavy-seed formation to account for the overmassive BHs already present at z > 7.

Figure 2 shows the time evolution of the most massive BH and the stellar mass of the source galaxy. Points with error bars denote the observed BH and stellar masses of the LRDs and the blue line traces the simulated relation. Before the formation of the massive seed, the BH mass lies far below the observed values. After the emergence of the heavy seed and its brief phase of super-Eddington growth, the simulated BH mass becomes comparable to the observations, reaching the BH-to-stellar mass ratio greater than 0.01. Subsequently, the BH and stellar masses coevolve, whereas the BH mass remains well above the local scaling relation indicated by the dashed lines.

The progenitor cloud does not form a single collapsing object but instead undergoes hierarchical fragmentation into a dense multiple system whose members experience sustained supergiant phases that suppress radiative feedback. To examine the feasibility of forming extremely massive-seed BHs, we performed a further high-resolution radiation-hydrodynamic simulation for the progenitor cloud of MBH1 that resolves gas densities up to 10⁸ cm⁻³. Figure 3a,b shows the gas density distribution during the protostellar evolution. Each protostar accretes mass at a rate of 0.1–1 M_⊙ yr⁻¹ and the most massive one reaches about 4 × 10⁵ M_⊙ about 2 Myr after the formation of the first protostar (Fig. 3c,d). Because of these high accretion rates, the stars remain in an inflated supergiant phase, thereby weakening radiative feedback and allowing continued growth²⁸. Figure 3e shows that the stars mostly stay in this supergiant phase, with radii larger than 100 R_⊙. Owing to the bursty nature of accretion, their radii also oscillate with time, repeating cycles of Kelvin–Helmholtz contraction and reinflation driven by entropy injection. This behaviour greatly reduces the ultraviolet emissivity and prevents photoevaporation, keeping gas available for subsequent BH growth. The stars collapse into BHs roughly 2 Myr after their formation, marking the end of the high-resolution simulation. At the time of collapse, dense gas remains tightly bound around the stars. These complementary calculations confirm that rapid accretion and limited feedback naturally lead to the formation of supermassive stars that ultimately collapse into heavy BH seeds.

Dense circum-BH disks in our simulation naturally reproduce the physical conditions required to explain the V-shaped continua and Balmer series absorption features observed in LRDs^6,7,8,29. Gas at densities n ≳ 10⁸ cm⁻³ enhances the n = 2 population of hydrogen, producing strong Balmer absorption^30,31. Such high densities imply large electron column densities and hence Thomson optical depths sufficient to broaden the Hα emission line to widths exceeding 1,000 km s⁻¹, consistent with observations^8,32,33. Notably, the Hα emission arises directly from the dense gas around the BHs without requiring an externally imposed broad-line region, closely resembling the recently proposed ‘BH-star’ systems^{9,29,31,33,34}.

To compare our heavy-seed BH system with the BH-star interpretation, we carry out a dedicated zoom-in simulation for MBH2, restarting at the moment of seed formation and resolving the circum-BH environment down to 500 au. Figure 4 shows snapshots at 26 and 523 kyr after BH formation. Shortly after formation, the BH is embedded in a dense, geometrically thick disk with hydrogen densities exceeding n_H > 10¹⁰ cm⁻³, sufficient to maintain a large n = 2 population (Fig. 4c). The resulting Hα luminosity within r < 10⁴ au is L_Hα = 1.5 × 10⁴³ erg s⁻¹, which is comparable to values inferred for LRDs^1,6. The Thomson optical depth reaches τ_e = 10.2 at this radius, large enough to generate strong broadening. Although this exceeds observationally inferred values of τ_e ≈ 0.1–1 (ref. ⁸), the effective optical depth for Hα photons is smaller because most of the emission originates near the outer disk (r ≈ 10⁴ au). In this early phase, the heavy-seed BH system therefore satisfies the physical conditions associated with the spectral features of a BH-star.

The obscured phase ends after several hundred kyr. At 523 kyr, the BH remains surrounded by dense gas (n_H ≳ 10⁸ cm⁻³) with a moderate Thomson depth of τ_e = 0.43, but the Hα luminosity declines to 5.3 × 10³⁹ erg s⁻¹, below typical LRD values (Fig. 4f). Producing strong Hα emission in this later stage requires further line emission from the inner (≲ 100 au) region, analogous to a compact broad-line region. The accretion rate remains near the Eddington limit, yielding an inner disk luminosity of about 10⁴⁴ erg s⁻¹, of which 0.1–1% is expected to emerge in Hα (ref. ³⁵). Thomson scattering in the outer disk then redistributes this radiation and produces broadened emission, yielding Hα luminosities of 10⁴¹–10⁴² erg s⁻¹, consistent with AGN. These results indicate that the system can naturally transition from an LRD-like, heavily obscured phase into a more AGN-like, less obscured accretion state on timescales of about 0.1–1 Myr following the bursty accretion episode.

Our results demonstrate that heavy-seed formation is a natural outcome of early structure formation in overdense regions, with profound implications for the observed abundance of early quasars and their gravitational-wave signatures. The overmassive BHs produced in these simulations represent plausible progenitors of the first quasars, bridging the gap between the birth of seed BHs and the luminous active nuclei observed by the JWST.

Our model shows that BHs can rapidly grow to about 10⁷ M_⊙ through the combination of direct collapse and subsequent super-Eddington accretion. This configuration—a massive BH embedded in a radiation-pressure-supported, optically thick envelope and fed at highly super-Eddington rates—closely resembles the ‘quasi-star’ solutions proposed in analytic models (for example, refs. ^34,36,37,38). Our simulations provide a fully cosmological realization of such quasi-star-like systems and follow their formation, migration and subsequent growth within a hierarchically assembling galaxy. This early phase of rapid growth helps overcome the long-standing difficulty that canonical massive BH seeds (about 10⁵ M_⊙) cannot efficiently accrete from hot, feedback-heated gas^26,39,40. Although our present simulation does not yet resolve the full galactic gas inflows required to sustain long-term Eddington accretion, it suggests that, once such large-scale supply is captured, the BHs could continue growing to about 10⁹ M_⊙, accounting for the luminous quasars already observed at z ≈ 7 (refs. ^41,42,43).

Although the super-Eddington accretion phase is short compared with the cosmic age, several heavy seeds form nearly simultaneously in the vicinity of massive halos, each experiencing a brief period of super-Eddington growth. The presence of several such sources increases the likelihood of observing galaxies in an active, LRD-like phase. Because our calculation covers only part of the regions surrounding massive halos, the true occurrence rate of heavy-seed formation may be even higher than inferred here. As these BHs continue to merge during hierarchical galaxy assembly, they are expected to generate strong gravitational-wave signals detectable by the Laser Interferometer Space Antenna (LISA). Extrapolating to cosmological scales, the predicted abundance of heavy seeds can account for the luminous quasar population observed by the JWST at z > 7. Their subsequent coalescences contribute gravitational-wave events with characteristic strains of h_c ≈ 10⁻¹⁷–10⁻¹⁶ at millihertz frequencies⁴⁴. Heavy-seed scenarios similar to ours therefore imply LISA event rates of tens to hundreds over a 3-year mission^45,46. Taken together, these findings suggest that the early Universe hosted numerous short-lived, rapidly growing BHs that profoundly influenced the assembly of the first massive galaxies.

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^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.

1. Harikane, Y. et al. A JWST/NIRSpec first census of broad-line AGNs at z = 4–7: detection of 10 faint AGNs with M_BH ~ 10⁶–10⁸M_⊙ and their host galaxy properties. Astrophys. J. 959, 39 (2023).

   Article  ADS  CAS  Google Scholar 

2. Übler, H. et al. GA-NIFS: a massive black hole in a low-metallicity AGN at z ~ 5.55 revealed by JWST/NIRSpec IFS. Astron. Astrophys. 677, A145 (2023).

   Article  Google Scholar 

3. Maiolino, R. et al. JADES. The diverse population of infant black holes at 4 < z < 11: merging, tiny, poor, but mighty. Astron. Astrophys. 691, A145 (2024).

   Article  CAS  Google Scholar 

4. Maiolino, R. et al. A small and vigorous black hole in the early Universe. Nature 627, 59–63 (2024).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

5. Taylor, A. J. et al. Broad-line AGNs at 3.5 < z < 6: the black hole mass function and a connection with little red dots. Astrophys. J. 986, 165 (2025).

   Article  ADS  CAS  Google Scholar 

6. Matthee, J. et al. Little red dots: an abundant population of faint active galactic nuclei at z ~ 5 revealed by the EIGER and FRESCO JWST surveys. Astrophys. J. 963, 129 (2024).

   Article  ADS  Google Scholar 

7. Hviding, R. E. et al. RUBIES: a spectroscopic census of little red dots: all point sources with v-shaped continua have broad lines. Astron. Astrophys. 702, A57 (2025).

   Article  CAS  Google Scholar 

8. Rusakov, V. Little red dots as young supermassive black holes in dense ionized cocoons. Nature 649, 574–579 (2026).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

9. Naidu, R. P. et al. A gas-enshrouded and gas-reddened black hole at cosmic dawn. Nature 656, 329–333 (2026).

10. de Graaff, A. et al. RUBIES: a complete census of the bright and red distant Universe with JWST/NIRSpec. Astron. Astrophys. 697, A189 (2025).

    Article  Google Scholar 

11. Lupi, A., Quadri, G., Volonteri, M., Colpi, M. & Regan, J. A. Sustained super-Eddington accretion in high-redshift quasars. Astron. Astrophys. 686, A256 (2024).

    Article  ADS  CAS  Google Scholar 

12. Trinca, A. et al. Episodic super-Eddington accretion as a clue to overmassive black holes in the early Universe. Preprint at https://doi.org/10.48550/arXiv.2412.14248 (2024).

13. Jeon, J. et al. Little red dots and their progenitors from direct collapse black holes. Astrophys. J. 998, 148 (2026).

    Article  ADS  Google Scholar 

14. Matthee, J. et al. Environmental evidence for overly massive black holes in low-mass galaxies and a black hole–halo mass relation at z ~ 5. Astrophys. J. 988, 246 (2025).

    Article  ADS  Google Scholar 

15. Bulichi, T.-E. et al. High-redshift AGN population in radiation-hydrodynamics simulations. Mon. Not. R. Astron. Soc. 544, 355–371 (2025).

    Article  ADS  CAS  Google Scholar 

16. Springel, V. E pur si muove: Galilean-invariant cosmological hydrodynamical simulations on a moving mesh. Mon. Not. R. Astron. Soc. 401, 791–851 (2010).

    Article  ADS  Google Scholar 

17. Ishiyama, T. & Hirano, S. A semianalytic framework of Population III and subsequent galaxy formation on cosmological N-body simulations. Astrophys. J. 994, 107 (2025).

    Article  ADS  Google Scholar 

18. Bromm, V. & Loeb, A. Formation of the first supermassive black holes. Astrophys. J. 596, 34–46 (2003).

    Article  ADS  Google Scholar 

19. Chon, S., Hosokawa, T. & Yoshida, N. Radiation hydrodynamics simulations of the formation of direct-collapse supermassive stellar systems. Mon. Not. R. Astron. Soc. 475, 4104–4121 (2018).

    Article  ADS  CAS  Google Scholar 

20. Regan, J. A. et al. The formation of very massive stars in early galaxies and implications for intermediate mass black holes. Open J. Astrophys. 3, 15 (2020).

    Article  ADS  Google Scholar 

21. Latif, M. A., Whalen, D. J., Khochfar, S., Herrington, N. P. & Woods, T. E. Turbulent cold flows gave birth to the first quasars. Nature 607, 48–51 (2022).

    Article  ADS  CAS  PubMed  Google Scholar 

22. Kiyuna, M., Hosokawa, T. & Chon, S. Sequential formation of supermassive stars and heavy seed BHs through the interplay of cosmological cold accretion and stellar radiative feedback. Mon. Not. R. Astron. Soc. 534, 3916–3935 (2024).

    Article  ADS  CAS  Google Scholar 

23. Begelman, M. C. Black holes in radiation-dominated gas: an analogue of the Bondi accretion problem. Mon. Not. R. Astron. Soc. 184, 53–67 (1978).

    Article  ADS  Google Scholar 

24. Ohsuga, K., Mori, M., Nakamoto, T. & Mineshige, S. Supercritical accretion flows around black holes: two-dimensional, radiation pressure-dominated disks with photon trapping. Astrophys. J. 628, 368–381 (2005).

    Article  ADS  CAS  Google Scholar 

25. Sadowski, A. & Narayan, R. Three-dimensional simulations of supercritical black hole accretion discs – luminosities, photon trapping and variability. Mon. Not. R. Astron. Soc. 456, 3929–3947 (2016).

    Article  ADS  Google Scholar 

26. Chon, S., Hosokawa, T. & Omukai, K. Cosmological direct-collapse black hole formation sites hostile for their growth. Mon. Not. R. Astron. Soc. 502, 700–713 (2021).

    Article  ADS  CAS  Google Scholar 

27. Larson, R. L. et al. A CEERS discovery of an accreting supermassive black hole 570 Myr after the Big Bang: identifying a progenitor of massive z > 6 quasars. Astrophys. J. Lett. 953, L29 (2023).

    Article  ADS  Google Scholar 

28. Hosokawa, T., Yorke, H. W., Inayoshi, K., Omukai, K. & Yoshida, N. Formation of primordial supermassive stars by rapid mass accretion. Astrophys. J. 778, 178 (2013).

    Article  ADS  Google Scholar 

29. de Graaff, A. et al. A remarkable ruby: absorption in dense gas, rather than evolved stars, drives the extreme Balmer break of a little red dot at z = 3.5. Astron. Astrophys. 701, A168 (2025).

    Article  Google Scholar 

30. Ferland, G. J. et al. The 2017 release of Cloudy. Rev. Mex. Astron. Astrofis. 53, 385–438 (2017).

    ADS  CAS  Google Scholar 

31. Inayoshi, K. & Maiolino, R. Extremely dense gas around little red dots and high-redshift active galactic nuclei: a nonstellar origin of the Balmer break and absorption features. Astrophys. J. Lett. 980, L27 (2025).

    Article  ADS  CAS  Google Scholar 

32. Chang, S.-J., Gronke, M., Matthee, J. & Mason, C. Impact of resonance, Raman, and Thomson scattering on hydrogen line formation in Little Red Dots. Mon. Not. R. Astron. Soc. 545, staf2131 (2026).

    Article  CAS  Google Scholar 

33. Torralba, A. et al. The warm outer layer of a little red dot as the source of [Fe II] and collisional Balmer lines with scattering wings. Astron. Astrophys. 707, A75 (2026).

    Article  CAS  Google Scholar 

34. Kido, D., Ioka, K., Hotokezaka, K., Inayoshi, K. & Irwin, C. M. Black hole envelopes in little red dots. Mon. Not. R. Astron. Soc. 544, 3407–3416 (2025).

    Article  ADS  CAS  Google Scholar 

35. Greene, J. E. & Ho, L. C. Estimating black hole masses in active galaxies using the Hα emission line. Astrophys. J. 630, 122–129 (2005).

    Article  ADS  CAS  Google Scholar 

36. Begelman, M. C., Rossi, E. M. & Armitage, P. J. Quasi-stars: accreting black holes inside massive envelopes. Mon. Not. R. Astron. Soc. 387, 1649–1659 (2008).

    Article  ADS  Google Scholar 

37. Volonteri, M. Formation of supermassive black holes. Astron. Astrophys. Rev. 18, 279–315 (2010).

    Article  ADS  Google Scholar 

38. Santarelli, A. D. Evolutionary tracks and spectral properties of quasi-stars and their correlation with little red dots. Astrophys. J. Lett. 998, L4 (2026).

    Article  ADS  CAS  Google Scholar 

39. Dubois, Y., Devriendt, J., Slyz, A. & Teyssier, R. Self-regulated growth of supermassive black holes by a dual jet-heating active galactic nucleus feedback mechanism: methods, tests and implications for cosmological simulations. Mon. Not. R. Astron. Soc. 420, 2662–2683 (2012).

    Article  ADS  Google Scholar 

40. Latif, M. A., Volonteri, M. & Wise, J. H. Early growth of typical high-redshift black holes seeded by direct collapse. Mon. Not. R. Astron. Soc. 476, 5016–5025 (2018).

    Article  ADS  CAS  Google Scholar 

41. Mortlock, D. J. et al. A luminous quasar at a redshift of z = 7.085. Nature 474, 616–619 (2011).

    Article  ADS  CAS  PubMed  Google Scholar 

42. Bañados, E. et al. An 800-million-solar-mass black hole in a significantly neutral Universe at a redshift of 7.5. Nature 553, 473–476 (2018).

    Article  ADS  PubMed  Google Scholar 

43. Wang, F. et al. A luminous quasar at redshift 7.642. Astrophys. J. Lett. 907, L1 (2021).

    Article  ADS  CAS  Google Scholar 

44. DeGraf, C. et al. High-redshift supermassive black hole mergers in simulations with dynamical friction modelling. Mon. Not. R. Astron. Soc. 527, 11766–11776 (2024).

    Article  ADS  Google Scholar 

45. Sesana, A., Volonteri, M. & Haardt, F. The imprint of massive black hole formation models on the LISA data stream. Mon. Not. R. Astron. Soc. 377, 1711–1716 (2007).

    Article  ADS  Google Scholar 

46. McCaffrey, J. et al. A heavy seed black hole mass function at high redshift – prospects for LISA. Open J. Astrophys. 8, 11 (2025).

    Article  ADS  Google Scholar 

47. Kocevski, D. D. et al. The rise of faint, red active galactic nuclei at z > 4: a sample of little red dots in the JWST extragalactic legacy fields. Astrophys. J. 986, 126 (2025).

    Article  ADS  CAS  Google Scholar 

48. Izumi, T. et al. Subaru High-z Exploration of Low-luminosity Quasars (SHELLQs). XIII. Large-scale feedback and star formation in a low-luminosity quasar at z = 7.07 on the local black hole to host mass relation. Astrophys. J. 914, 36 (2021).

    Article  ADS  CAS  Google Scholar 

49. Kormendy, J. & Ho, L. C. Coevolution (or not) of supermassive black holes and host galaxies. Annu. Rev. Astron. Astrophys. 51, 511–653 (2013).

    Article  ADS  CAS  Google Scholar 

50. Reines, A. E. & Volonteri, M. Relations between central black hole mass and total galaxy stellar mass in the local universe. Astrophys. J. 813, 82 (2015).

    Article  ADS  Google Scholar 

51. Ishiyama, T. et al. The Uchuu simulations: Data Release 1 and dark matter halo concentrations. Mon. Not. R. Astron. Soc. 506, 4210–4231 (2021).

    Article  ADS  CAS  Google Scholar 

52. Hahn, O. & Abel, T. Multi-scale initial conditions for cosmological simulations. Mon. Not. R. Astron. Soc. 415, 2101–2121 (2011).

    Article  ADS  CAS  Google Scholar 

53. Planck Collaboration et al. Planck 2018 results. VI. Cosmological parameters. Astron. Astrophys. 641, A6 (2020).

    Article  ADS  Google Scholar 

54. Behroozi, P. S., Wechsler, R. H. & Wu, H.-Y. The ROCKSTAR phase-space temporal halo finder and the velocity offsets of cluster cores. Astrophys. J. 762, 109 (2013).

    Article  ADS  Google Scholar 

55. Behroozi, P. S. et al. Gravitationally consistent halo catalogs and merger trees for precision cosmology. Astrophys. J. 763, 18 (2013).

    Article  ADS  Google Scholar 

56. Agarwal, B. et al. Ubiquitous seeding of supermassive black holes by direct collapse. Mon. Not. R. Astron. Soc. 425, 2854–2871 (2012).

    Article  ADS  Google Scholar 

57. Visbal, E., Bryan, G. L. & Haiman, Z. Self-consistent semianalytic modeling of feedback during primordial star formation and reionization. Astrophys. J. 897, 95 (2020).

    Article  ADS  Google Scholar 

58. Kauffmann, G., White, S. D. M. & Guiderdoni, B. The formation and evolution of galaxies within merging dark matter haloes. Mon. Not. R. Astron. Soc. 264, 201–218 (1993).

    Article  ADS  CAS  Google Scholar 

59. Cole, S., Lacey, C. G., Baugh, C. M. & Frenk, C. S. Hierarchical galaxy formation. Mon. Not. R. Astron. Soc. 319, 168–204 (2000).

    Article  ADS  CAS  Google Scholar 

60. Hirano, S., Hosokawa, T., Yoshida, N., Omukai, K. & Yorke, H. W. Primordial star formation under the influence of far ultraviolet radiation: 1540 cosmological haloes and the stellar mass distribution. Mon. Not. R. Astron. Soc. 448, 568–587 (2015).

    Article  ADS  CAS  Google Scholar 

61. Johnson, J. L., Dalla, V. C. & Khochfar, S. The First Billion Years project: the impact of stellar radiation on the co-evolution of Populations II and III. Mon. Not. R. Astron. Soc. 428, 1857–1872 (2013).

    Article  ADS  CAS  Google Scholar 

62. Scalo, J. The IMF revisited: a case for variations. In The Stellar Initial Mass Function (38th Herstmonceux Conference) (eds Gilmore, G. & Howell, D.) 201–236 (Astronomical Society of the Pacific, 1998).

63. Barkana, R. & Loeb, A. In the beginning: the first sources of light and the reionization of the universe. Phys. Rep. 349, 125–238 (2001).

    Article  ADS  CAS  Google Scholar 

64. Greif, T. H. & Bromm, V. Two populations of metal-free stars in the early Universe. Mon. Not. R. Astron. Soc. 373, 128–138 (2006).

    Article  ADS  CAS  Google Scholar 

65. Kulkarni, M., Visbal, E. & Bryan, G. L. The critical dark matter halo mass for Population III star formation: dependence on Lyman–Werner radiation, baryon-dark matter streaming velocity, and redshift. Astrophys. J. 917, 40 (2021).

    Article  ADS  CAS  Google Scholar 

66. Chon, S., Hirano, S., Hosokawa, T. & Yoshida, N. Cosmological simulations of early black hole formation: halo mergers, tidal disruption, and the conditions for direct collapse. Astrophys. J. 832, 134 (2016).

    Article  ADS  Google Scholar 

67. Truelove, J. K. et al. The Jeans condition: a new constraint on spatial resolution in simulations of isothermal self-gravitational hydrodynamics. Astrophys. J. 489, L179 (1997).

    Article  ADS  Google Scholar 

68. Matsukoba, R., Tanaka, K. E. I., Omukai, K., Vorobyov, E. I. & Hosokawa, T. Protostellar-disc fragmentation across all metallicities. Mon. Not. R. Astron. Soc. 515, 5506–5522 (2022).

    Article  ADS  CAS  Google Scholar 

69. Glover, S. C. O. Simulating the formation of massive seed black holes in the early Universe – I. An improved chemical model. Mon. Not. R. Astron. Soc. 451, 2082–2096 (2015).

    Article  ADS  CAS  Google Scholar 

70. Susa, H. Smoothed particle hydrodynamics coupled with radiation transfer. Publ. Astron. Soc. Jpn. 58, 445–460 (2006).

    Article  ADS  CAS  Google Scholar 

71. Chon, S. & Latif, M. A. The impact of ionizing radiation on the formation of a supermassive star in the early Universe. Mon. Not. R. Astron. Soc. 467, 4293–4303 (2017).

    Article  ADS  CAS  Google Scholar 

72. Shang, C., Bryan, G. L. & Haiman, Z. Supermassive black hole formation by direct collapse: keeping protogalactic gas H₂ free in dark matter haloes with virial temperatures T_vir ≳ 10⁴ K. Mon. Not. R. Astron. Soc. 402, 1249–1262 (2010).

    Article  ADS  CAS  Google Scholar 

73. Sugimura, K., Omukai, K. & Inoue, A. K. The critical radiation intensity for direct collapse black hole formation: dependence on the radiation spectral shape. Mon. Not. R. Astron. Soc. 445, 544–553 (2014).

    Article  ADS  CAS  Google Scholar 

74. Latif, M. A. et al. A UV flux constraint on the formation of direct collapse black holes. |Mon. Not. R. Astron. Soc. 443, 1979–1987 (2014).

    Article  ADS  CAS  Google Scholar 

75. Regan, J. A., Johansson, P. H. & Wise, J. H. Forming supermassive black hole seeds under the influence of a nearby anisotropic multifrequency source. Mon. Not. R. Astron. Soc. 459, 3377–3394 (2016).

    Article  ADS  CAS  PubMed  Google Scholar 

76. Wolcott-Green, J., Haiman, Z. & Bryan, G. L. Photodissociation of H₂ in protogalaxies: modelling self-shielding in three-dimensional simulations. Mon. Not. R. Astron. Soc. 418, 838–852 (2011).

    Article  ADS  CAS  Google Scholar 

77. Hosokawa, T. & Omukai, K. Evolution of massive protostars with high accretion rates. Astrophys. J. 691, 823–846 (2009).

    Article  ADS  CAS  Google Scholar 

78. Haemmerlé, L., Woods, T. E., Klessen, R. S., Heger, A. & Whalen, D. J. The evolution of supermassive Population III stars. Mon. Not. R. Astron. Soc. 474, 2757–2773 (2018).

    Article  ADS  Google Scholar 

79. Nandal, D. et al. Critical accretion rates for rapidly growing massive Population III stars. Astron. Astrophys. 677, A155 (2023).

    Article  CAS  Google Scholar 

80. Nandal, D. & Chon, S. Growth of metal-enriched supermassive stars by accretion and collisions. Astrophys. J. 999, 110 (2026).

    Article  ADS  CAS  Google Scholar 

81. Hosokawa, T., Omukai, K. & Yorke, H. W. Rapidly accreting supergiant protostars: embryos of supermassive black holes? Astrophys. J. 756, 93 (2012).

    Article  ADS  Google Scholar 

82. Sakurai, Y., Hosokawa, T., Yoshida, N. & Yorke, H. W. Formation of primordial supermassive stars by burst accretion. Mon. Not. R. Astron. Soc. 452, 755–764 (2015).

    Article  ADS  CAS  Google Scholar 

83. Heger, A. & Woosley, S. E. The nucleosynthetic signature of Population III. Astrophys. J. 567, 532–543 (2002).

    Article  ADS  CAS  Google Scholar 

84. Toyouchi, D., Inayoshi, K., Li, W., Haiman, Z. & Kuiper, R. Radiative feedback on supermassive star formation: the massive end of the Population III initial mass function. Mon. Not. R. Astron. Soc. 518, 1601–1616 (2023).

    Article  ADS  CAS  Google Scholar 

85. Fuller, G. M., Woosley, S. E. & Weaver, T. A. The evolution of radiation-dominated stars. I - Nonrotating supermassive stars. Astrophys. J. 307, 675–686 (1986).

    Article  ADS  CAS  Google Scholar 

86. Shibata, M. & Shapiro, S. L. Collapse of a rotating supermassive star to a supermassive black hole: fully relativistic simulations. Astrophys. J. 572, L39–L43 (2002).

    Article  ADS  CAS  Google Scholar 

87. Woods, T. E., Heger, A., Whalen, D. J., Haemmerlé, L. & Klessen, R. S. On the maximum mass of accreting primordial supermassive stars. Astrophys. J. Lett. 842, L6 (2017).

    Article  ADS  Google Scholar 

88. Uchida, H., Shibata, M., Yoshida, T., Sekiguchi, Y. & Umeda, H. Gravitational collapse of rotating supermassive stars including nuclear burning effects. Phys. Rev. D 96, 083016 (2017).

    Article  ADS  Google Scholar 

89. Nandal, D. et al. The evolution of accreting Population III stars at 10⁻⁶–10³ M_⊙ yr⁻¹. Astron. Astrophys. 689, A351 (2024).

    Article  CAS  Google Scholar 

90. Umeda, H., Hosokawa, T., Omukai, K. & Yoshida, N. The final fates of accreting supermassive stars. Astrophys. J. Lett. 830, L34 (2016).

    Article  ADS  Google Scholar 

91. Sazonov, S. Y., Ostriker, J. P. & Sunyaev, R. A. Quasars: the characteristic spectrum and the induced radiative heating. Mon. Not. R. Astron. Soc. 347, 144–156 (2004).

    Article  ADS  CAS  Google Scholar 

92. McKinney, J. C., Tchekhovskoy, A., Sadowski, A. & Narayan, R. Three-dimensional general relativistic radiation magnetohydrodynamical simulation of super-Eddington accretion, using a new code HARMRAD with M1 closure. Mon. Not. R. Astron. Soc. 441, 3177–3208 (2014).

    Article  ADS  Google Scholar 

93. Inayoshi, K., Haiman, Z. & Ostriker, J. P. Hyper-Eddington accretion flows on to massive black holes. Mon. Not. R. Astron. Soc. 459, 3738–3755 (2016).

    Article  ADS  CAS  Google Scholar 

94. Watarai, K.-y., Mizuno, T. & Mineshige, S. Slim-disk model for ultraluminous X-ray sources. Astrophys. J. Lett. 549, L77–L80 (2001).

    Article  ADS  Google Scholar 

95. Regan, J. A. et al. Super-Eddington accretion and feedback from the first massive seed black holes. Mon. Not. R. Astron. Soc. 486, 3892–3906 (2019).

    Article  ADS  CAS  Google Scholar 

96. Massonneau, W., Volonteri, M., Dubois, Y. & Beckmann, R. S. How the super-Eddington regime regulates black hole growth in high-redshift galaxies. Astron. Astrophys. 670, A180 (2023).

    Article  ADS  Google Scholar 

97. Omukai, K. Primordial star formation under far-ultraviolet radiation. Astrophys. J. 546, 635–651 (2001).

    Article  ADS  CAS  Google Scholar 

98. Becerra, F., Marinacci, F., Bromm, V. & Hernquist, L. E. Assembly of supermassive black hole seeds. Mon. Not. R. Astron. Soc. 480, 5029–5045 (2018).

    ADS  CAS  Google Scholar 

99. Chon, S. & Omukai, K. Formation of supermassive stars and dense star clusters in metal-poor clouds exposed to strong FUV radiation. Mon. Not. R. Astron. Soc. 539, 2561–2582 (2025).

    Article  ADS  Google Scholar 

100. Dijkstra, M., Ferrara, A. & Mesinger, A. Feedback-regulated supermassive black hole seed formation. Mon. Not. R. Astron. Soc. 442, 2036–2047 (2014).

     Article  ADS  Google Scholar 

101. Ventura, E. M., Qin, Y., Balu, S. & Wyithe, J. S. B. Semi-analytic modelling of Pop. III star formation and metallicity evolution – I. Impact on the UV luminosity functions at z = 9–16. Mon. Not. R. Astron. Soc. 529, 628–646 (2024).

     Article  ADS  CAS  Google Scholar 

102. Yoshida, N., Abel, T., Hernquist, L. & Sugiyama, N. Simulations of early structure formation: primordial gas clouds. Astrophys. J. 592, 645–663 (2003).

     Article  ADS  CAS  Google Scholar 

103. Wise, J. H. et al. Formation of massive black holes in rapidly growing pre-galactic gas clouds. Nature 566, 85–88 (2019).

     Article  ADS  CAS  PubMed  Google Scholar 

104. Bonnor, W. B. Boyle’s Law and gravitational instability. Mon. Not. R. Astron. Soc. 116, 351–359 (1956).

     Article  ADS  MathSciNet  Google Scholar 

105. Ebert, R. Über die Verdichtung von H I-Gebieten. Mit 5 Textabbildungen. Z. Astrophys. 37, 217–232 (1955).

     ADS  Google Scholar 

106. Yoshida, N., Oh, S. P., Kitayama, T. & Hernquist, L. Early cosmological H II/He III regions and their impact on second-generation star formation. Astrophys. J. 663, 687–707 (2007).

     Article  ADS  CAS  Google Scholar 

107. Chon, S. & Omukai, K. Supermassive star formation via super competitive accretion in slightly metal-enriched clouds. Mon. Not. R. Astron. Soc. 494, 2851–2860 (2020).

     Article  ADS  CAS  Google Scholar 

108. Inayoshi, K., Onoue, M., Sugahara, Y., Inoue, A. K. & Ho, L. C. The age of discovery with the James Webb Space Telescope: excavating the spectral signatures of the first massive black holes. Astrophys. J. Lett. 931, L25 (2022).

     Article  ADS  Google Scholar 

109. Hviding, R. E. et al. The X-ray dot: exotic dust or a late-stage little red dot? Astrophys. J. Lett. 1000, L18 (2026).

     Article  ADS  CAS  Google Scholar 

110. Ananna, T. T., Bogdán, Á, Kovács, O. E., Natarajan, P. & Hickox, R. C. X-ray view of little red dots: do they host supermassive black holes? Astrophys. J. Lett. 969, L18 (2024).

     Article  ADS  CAS  Google Scholar 

111. Yue, M. et al. Stacking X-ray observations of “Little Red Dots”: implications for their active galactic nucleus properties. Astrophys. J. Lett. 974, L26 (2024).

     Article  ADS  CAS  Google Scholar 

112. Sacchi, A. & Bogdán, Á Chandra rules out super-Eddington accretion models for little red dots. Astrophys. J. Lett. 989, L30 (2025).

     Article  ADS  Google Scholar 

113. Alvarez, M. A., Wise, J. H. & Abel, T. Accretion onto the first stellar-mass black holes. Astrophys. J. 701, L133–L137 (2009).

     Article  ADS  CAS  Google Scholar 

114. Jeon, M. et al. The first galaxies: assembly with black hole feedback. Astrophys. J. 754, 34 (2012).

     Article  ADS  Google Scholar 

115. Jeon, J., Liu, B., Bromm, V. & Finkelstein, S. L. Observability of low-luminosity AGNs in the early Universe with JWST. Mon. Not. R. Astron. Soc. 524, 176–187 (2023).

     Article  ADS  Google Scholar 

116. O’Shea, B. W. & Norman, M. L. Population III star formation in a ΛCDM universe. II. Effects of a photodissociating background. Astrophys. J. 673, 14–33 (2008).

     Article  ADS  Google Scholar 

117. Chon, S., Hirano, S., Ishiyama, T., Chang, S.-J. & Springel, V. Arepo suite for DCBH formation. Zenodo https://doi.org/10.5281/zenodo.21242352 (2026).

118. Chon, S. particle_visualization. Zenodo https://doi.org/10.5281/zenodo.21761582 (2026).

Download references

Authors and Affiliations

1. Max-Planck-Institut für Astrophysik, Garching, Germany

   Sunmyon Chon, Seok-Jun Chang & Volker Springel

2. Department of Applied Physics, Faculty of Engineering, Kanagawa University, Yokohama, Japan

   Shingo Hirano

3. Department of Astronomy, School of Science, University of Tokyo, Tokyo, Japan

   Shingo Hirano

4. Digital Transformation Enhancement Council, Chiba University, Chiba, Japan

   Tomoaki Ishiyama

Contributions

S.C. led the project, performed the hydrodynamical simulations, analysed the simulation data, produced the figures and wrote the manuscript. T.I. provided the initial conditions for black hole formation and contributed to the selection of the target halos. S.H. contributed to the scientific motivation and interpretation of the study, including the connection between the semi-analytic model of T.I. and S.H. and the hydrodynamical simulations. V.S. contributed to the development and implementation of the simulation code, advised on the numerical methods and supervised the overall project. S.-J.C. contributed to the radiation-transfer modelling and to the interpretation of the results in connection with observations. All authors discussed the results and commented on the manuscript.

Corresponding author

Correspondence to Sunmyon Chon.

Competing interests

The authors declare no competing interests.

Peer review information

Nature thanks John Regan and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions