Making the first supermassive black holes Supermassive black holes existed less than a billion years after the Big Bang. Because black holes can grow at a maximum rate that depends on their current mass, it has been difficult to understand how such massive black holes could have formed so quickly. Hirano et al. developed simulations to show that streaming motions—velocity offsets between the gas and dark matter components—could have produced black holes with tens of thousands of solar masses in the early universe. That's big enough to grow into the supermassive black holes that we observe today. Science, this issue p. 1375

Abstract The origin of super-massive black holes in the early universe remains poorly understood. Gravitational collapse of a massive primordial gas cloud is a promising initial process, but theoretical studies have difficulty growing the black hole fast enough. We report numerical simulations of early black hole formation starting from realistic cosmological conditions. Supersonic gas motions left over from the Big Bang prevent early gas cloud formation until rapid gas condensation is triggered in a protogalactic halo. A protostar is formed in the dense, turbulent gas cloud, and it grows by sporadic mass accretion until it acquires 34,000 solar masses. The massive star ends its life with a catastrophic collapse to leave a black hole—a promising seed for the formation of a monstrous black hole.

Recent discoveries of super-massive black holes (SMBHs) at redshift z ~ 7, when the universe was just 5% of its present age, pose a serious challenge to the theory of black hole (BH) formation and evolution (1). The physical mechanisms that form the BH seeds and drive their growth are not yet known but must explain how an initial seed BH can reach a mass of 10 billion times that of the Sun ( ) within 1 billion years after the Big Bang. It is thought that the mass growth is a self-regulating process, limited by the so-called Eddington rate that is proportional to the BH mass; therefore, starting from a massive BH holds a key to the rapid formation of SMBHs.

A previously proposed physical model posits that an early BH forms through direct gravitational collapse of a large primordial gas cloud with a mass of about (2). However, the model must invoke several critical conditions, such as the formation of a bright galaxy in the immediate vicinity (3–6). The overall occurrence rate of such a supposedly rare event in a cosmological volume remains rather uncertain. We present an ab initio simulation of the formation of massive BHs and examine the importance of the supersonic streaming motion between dark matter and ordinary baryonic matter in the early universe (7). In regions with a large streaming velocity, gas condensation—and hence star formation—is suppressed until a deep gravitational potential is generated by a clump of dark matter with mass (8), where similar conditions to that of the conventional direct-collapse model might be realized (9).

We performed cosmological hydrodynamics simulations of early structure formation that incorporate the relative bulk velocity. We selected target dark-matter halos whose dynamical properties were consistent with those of the host halos of SMBHs (10, 11). We then regenerated the initial conditions for three zoom-in simulations with a fixed streaming velocity of three times the root mean square (RMS) value (11). In the following material, we concentrate on describing one case, Run-B, which has , corresponding to 1.5 times the cosmic mean value, where is the RMS density fluctuation in a sphere of radius 8 Mpc. We also performed two other runs: Run-A, with a larger fluctuation amplitude ; and Run-C, starting from a different realization of the density field with (table S1).

At early epochs in Run-B, the streaming motions prevented gas cloud collapse in small-mass dark-matter halos with (12, 13), which would otherwise host the first generation of stars (14). Dark-matter halos grew hierarchically through mergers and accretion, and a massive halo with 2.2 × 107 was assembled when the cosmological redshift was 30.5 (Fig. 1A). The host halo was about 100 times as heavy as that found in previous simulations without the streaming motions (table S1). The strong gravity trapped the streaming gas (Fig. 1B), and the excess momentum generated a wide wedge-shaped structure (Fig. 1C).

Fig. 1 Large-scale density distribution and the structure around an accreting protostar. (A) Projected density distribution of dark-matter component around the star-forming cloud at redshift z = 30.5 in Run-B. The box size is 2500 pc on a side. The virial mass of the main dark-matter halo located at the center is 2.2 × 107 . (B to D) Projected density distribution of the gas component in regions of 2500, 500, 100, and 10 pc on a side, from left to right. The horizontal arrow in (B) shows the direction of the initial supersonic gas stream. The dashed circle in (D) indicates the Jeans length, within which the cloud is gravitationally unstable given its mass of 26,000 . (E to H) Evolution of the temperature and density structure in the protostellar accretion phase after the protostar formation. Colored in white, red, and magenta are the isocontours of gas density (at 106, 105, and 3 × 104 cm–3), photoionized hydrogen abundance with ≥50% (H ii region), and the number fraction of hydrogen molecules with ≥0.2%.

The trapped gas was first shock-heated to a temperature of ~10,000 K, and the temperature soon decreased to ~8000 K by atomic-hydrogen (H) cooling (Fig. 2). The compressed gas further cooled through molecular-hydrogen (H 2 ) emission and formed a dense and cold gas cloud (fig. S1). The cloud was still supported against gravitational collapse by thermal pressure and also by turbulence induced by streaming motions. For such a dynamically supported gas with density ρ, the Jeans mass is given by , where G is the gravitational constant, and , with and denoting the sound speed and relative velocity between baryon and cold dark-matter components, respectively (12, 15). The assembled cloud became gravitationally unstable when its mass exceeded 26,000 at particle number density 4000 cm–3 and gas temperature 400 K (Fig. 1D and fig. S2). The rapidly inflowing gas was accumulated on the contracting cloud and produced a massive and dense envelope (fig. S3A). During the subsequent collapse, a fully molecular cloud formed by efficient H 2 formation via three-body reactions at 108 cm–3 (16), causing the molecular core to collapse further. After this epoch, the temperature evolution of the core was similar to that in the ordinary primordial star formation (17), but the whole cloud contracted more rapidly, and the instantaneous gas mass accretion rate exceeded in the outer part of the cloud (figs. S3 and S4).

Fig. 2 Thermal evolution of the collapsing cloud. The red, blue, green, and black lines show our Run-A, B, C, and Ref, respectively. In each run, the cloud became Jeans-unstable at the points marked by the circles. The background colored regions are distinguished by the major coolant: atomic hydrogen (H) and molecular-hydrogen (H 2 ). In all the cases, H 2 cooling operated at densities greater than ~100 cm–3.

The conventional model of direct BH formation (9, 18) assumes that the gas evolves nearly isothermally at a high temperature of ~8000 K. That model invokes a few necessary conditions so that the gas can evolve on a high-temperature track to maintain a large gas mass accretion rate in the later accretion phase. Often, photodissociation of hydrogen molecules by a nearby galaxy is assumed (4). The massive gas cloud in our simulation cooled first by H cooling and then by H 2 cooling. The temperature drop at a density of 100 cm–3 and the onset of H 2 cooling have been considered as a failure of direct collapse, possibly leading to cloud fragmentation (19), but our simulations show that the cloud collapse continues in a highly dynamical manner.

We stopped our cosmological simulation when a tiny protostar with a mass of ~0.01 was formed at the gas cloud center. We then turned to a three-dimensional gravito-radiation-hydrodynamic simulation coupled with a self-consistent modeling of the stellar evolution (11) in order to follow directly the complex interplay between the gas mass accretion and the stellar radiation feedback that reduces the mass accretion rate (Fig. 1, E to H). When the mass accretion rate is greater than a critical value of , the excess entropy carried by rapidly accreting gas causes the stellar envelope to inflate. A star with an extended envelope has a low effective temperature of ~6000 K and hence does not emit copious amounts of ultraviolet (UV) photons. The resulting radiative feedback on the surrounding gas is weak and cannot halt the gas accretion. The central star in our simulation grew rapidly to a mass of in the first 2000 years (Fig. 3A).

Fig. 3 Time evolution of stellar properties. (A) Stellar masses, (B) mass accretion rates, (C) stellar radii, and (D) extreme ultraviolet (EUV) (with hv ≥ 13.6 eV; solid) and far ultraviolet (FUV) (with 11.2 eV ≤ hv ≤ 13.6 eV; dashed) emissivities, where h is the Planck constant and v is the frequency of a photon. The red, blue, and green lines represent our Run-A, B, and C, respectively. When the mass accretion rate was larger than the critical value of year–1 (Run-A), the star expanded continuously during the accretion phase, and the FUV/EUV radiation could not halt gas accretion. When the gas accretion rate fluctuated around the critical value due to sporadic accretion (Run-B), the stellar radius and the resulting UV emissivity, repeated rises and falls (21–23). The accretion rate stayed well below the critical value for a sufficiently long time in its final phase in Run-C, and the protostar entered the zero-age main sequence.

Gas accretion onto the star occurred sporadically owing to the combined effect of the temperature structure of the infalling gas and mass accretion through the circumstellar disk, where the gravitational torque dominated the angular momentum transport (19, 20). The mass accretion rate fluctuated around the critical value of year–1 for the first 300,000 years (Fig. 3B). During the early phase, the star contracted and its UV emissivity increased temporarily when the accretion rate fell below the critical value (Fig. 3, C and D). However, gravitational fragmentation of the cloud and subsequent merging of the fragments caused accretion bursts, and the star quickly recovered to have an extended envelope (21–23). Because successive accretion bursts occurred at short time intervals in the late accretion phase, the stellar envelope remained extended because there was insufficient time for it to contract.

The hydrogen molecules formed along the dense filament (Fig. 1E) were dissociated by the far-UV radiation emitted by the central star. A bipolar ionized atomic hydrogen (H ii) region appeared temporally (red contour in Fig. 1F) due to the increasing intensity of the stellar UV radiation. The surrounding dense gas filament fragmented to yield a gas clump (Fig. 1G), but it did not cool and contract further because it contained few hydrogen molecules, and gradually approached the central star, avoiding the bipolar H ii region. The clump finally reached the center (Fig. 1H), activated an accretion burst (Fig. 3B) onto the star, and weakened the radiative feedback. After the protostellar mass reached ~10,000 , it evolved as a stable supergiant protostar without radial contraction and continued growing steadily up to 34,000 . Run-A and Run-C also showed the formation of massive primordial stars with 100,000 and 4400 the mass of the Sun at the simulation end. Such very massive accreting stars undergo direct collapse to produce equally massive BHs owing to the general relativistic instability or exhaustion of nuclear fuel, depending on the accretion rate (11, 24).

A 34,000 BH formed at z = 30.5 must grow at about 55% of the canonical Eddington rate until z = 7.1 to acquire a mass of 2 × 109 , matching the estimated mass of the SMBH observed in a luminous quasar (1, 25, 26). The number density of the intermediate mass BHs formed by this mechanism is estimated to be ~1 per cubic comoving Gpc (11), which is similar to the abundance of the observed SMBHs (27).

Unlike previous studies of direct-collapse BH formation, our simulations did not assume additional conditions, such as the existence of strong radiation sources (28) and/or high-speed collisions of gas clouds (29) (see also supplementary text). Other than the primeval density fluctuations, whose statistical properties are well understood from both theory and observation, the only newly introduced element in our study is the relative velocity between the dark matter and baryonic components. The baryonic streaming motions are intrinsically generated in the early universe according to the standard model of structure formation (7). Therefore, our ab initio cosmological simulations show a viable formation path for massive BHs.

Supplementary Materials www.sciencemag.org/content/357/6358/1375/suppl/DC1 Materials and Methods Supplementary Text Figs. S1 to S9 Table S1 References (30–51)