The Birth of Binary Direct-collapse Black Holes

More than 300 quasars have now been discovered at z > 6 (e.g., Fan et al. 2003), including 7 at z > 7 (Mortlock et al. 2011; Bañados et al. 2018; Matsuoka et al. 2019). The formation of such massive black holes (BHs) less than 1 Gyr after the Big Bang poses serious challenges for paradigms of early structure formation. A number of processes have been proposed for the origins of these quasars: the collapse of Population III stars to BHs at z ∼ 20–25, runaway collisions in dense nuclear clusters at z ∼ 10–20 that build up a single massive star that collapses to a BH, and the formation of supermassive stars (SMSs) in atomically cooling halos at z ∼ 15–20 that die as direct-collapse black holes (DCBHs).

Population III remnant BHs are expected to range in mass from a few tens to hundreds of solar masses at birth (e.g., Hirano et al. 2014) but are born in low ambient densities that preclude their initial rapid growth (Whalen et al. 2004; Johnson & Bromm 2007), and some are ejected from their halos at high velocities by natal kicks (Whalen & Fryer 2012). They in principle could reach 10⁹ ${M}_{\odot }$ by z ∼ 7 with episodes of super-Eddington accretion at duty cycles of just a few percent (Lupi et al. 2016; Pezzulli et al. 2016) but have not encountered high enough densities to trigger such growth in any cosmological simulation to date (Alvarez et al. 2009; Smith et al. 2018). Runaway stellar collisions in marginally enriched dense stellar clusters can create BHs of up to a few thousand solar masses (Devecchi et al. 2012; Latif et al. 2016a; Reinoso et al. 2018), but even these objects may not be massive enough to become quasars by z > 6 (Smidt et al. 2018).

For these reasons DCBHs have become the leading contenders for the seeds of the first quasars. Catastrophic baryon collapse in atomically cooling halos leads to initial infall rates of ∼0.01–1 ${M}_{\odot }$ yr⁻¹. Standalone models of stellar evolution have shown that if such rates persist they would build up cool, red 100,000–300,000 ${M}_{\odot }$ stars before they collapse to DCBHs via general relativistic instability (Hosokawa et al. 2013; Umeda et al. 2016; Haemmerlé et al. 2018a, 2018b). The low ionizing UV fluxes of these stars cannot slow down accretion onto themselves so DCBHs form in the dense environments that create them and can therefore grow much faster at birth. But for these halos to reach masses of 10⁷–10⁸ ${M}_{\odot }$ and virial temperatures of ∼10⁴ K without having first formed a less massive Population III star via H₂ cooling, they must grow either in the presence of Lyman–Werner (LW) UV sources that sterilize them of H₂ (e.g., Sugimura et al. 2014; Latif et al. 2015; Agarwal et al. 2016) or in highly supersonic baryon streaming motions that delay the collapse of the halo even if H₂ is present (Hirano et al. 2017; Schauer et al. 2017). Once collapse begins it proceeds quickly because of high densities, temperatures, and therefore sound speeds, c_s, in the gas ($\dot{M}\sim {c}_{{\rm{s}}}^{3}/G\sim 0.1\,{M}_{\odot }/\mathrm{yr}{\left(T/8000{\rm{K}}\right)}^{3/2}$).

Numerical simulations of the collapse of atomically cooling halos at high redshifts have steadily improved over the past decade but remain a trade-off between resolution and evolution time. The original simulations either resolved sub-au scales that could only follow the formation of the hydrostatic protostar (but not the atomically cooled disk around it; Wise et al. 2008) or 0.01 pc scales that captured the formation of the disk but still could only follow its evolution for a few dynamical times (Regan & Haehnelt 2009). These studies found large infall rates at early times but could not determine how long they lasted. Later work at high resolution and somewhat longer evolution times found that large accretion rates continued down to scales approaching those of the supermassive star itself but did not run for nearly enough times to follow its evolution (Latif et al. 2013a, 2016b; Regan et al. 2014).

The introduction of sink particles and pressure floors (Machacek et al. 2001) at the highest resolutions extended simulation times to a few tens of thousands of years (Latif et al. 2013c; Shlosman et al. 2016; Becerra et al. 2018; Regan & Downes 2018), confirming that accretion rates at the smallest scales remained high. Most recently, Chon et al. (2018) and Regan & Downes (2018) employed sink particles with radiative feedback from the protostar with a simple treatment of its evolution to follow collapse for 100 kyr and 250 kyr, respectively. They confirmed that radiation from the protostar was unable to prevent high accretion rates and found some small-scale fragmentation in the disk at later times. In particular, Chon et al. (2018) found that almost half of the fragments in one of the simulated halos are formed in binaries that they suspected may become supermassive binaries, corroborating idealized simulations by Bromm & Loeb (2003). However, neither study evolved the disks for long enough times to determine if any of the fragments became stars or were simply subsumed into the central object at later times. Previous work also only considered one or a few halos that were not parameterized by either spin parameter or assembly history, so how these factors influence the masses or numbers of supermassive stars in the halos remains unknown.

We follow here the collapse of atomically cooled halos at high redshift for timescales 12 times longer than in any comparable simulation to determine the final fates of any fragments that form in their accretion disks, how many supermassive stars might result, and their masses at collapse to DCBHs. Our halos were chosen to have a range of spin parameters that bracket their likely values at early epochs, and we have tallied accretion rates at the centers of the disks for later use in supermassive stellar evolution models in the cosmological flows that create the stars. In Section 2 we describe our numerical methods and simulations and discuss our results in Section 3. We examine the consequences of our results and conclude in Section 4.

We model the collapse of atomically cooled halos with the Enzo adaptive mesh refinement cosmology code (Bryan et al. 2014). Enzo utilizes an N-body adaptive particle-mesh scheme to evolve dark matter (DM), a third-order piecewise-parabolic method for fluid dynamics, a multigrid Poisson solver for calculating self-gravity, and a nonequilibrium reaction network to evolve primordial gas chemistry (Anninos et al. 1997). Simulations are initialized with Gaussian primordial density fluctuations at z = 150 with MUSIC (Hahn & Abel 2011) with cosmological parameters from the second-year Planck best-fit lowP+lensing+BAO+JLA+H₀: Ω_M = 0.308, Ω_Λ = 0.691, Ω_b = 0.0223, h = 0.677, σ₈ = 0.816, and n = 0.968 (Planck Collaboration et al. 2016). To approximate the presence of a strong LW background we only evolve ${\rm{H}},{{\rm{H}}}^{+},{\rm{H}}{\rm{e}},{{\rm{H}}{\rm{e}}}^{+},{{\rm{H}}{\rm{e}}}^{++},\,{\rm{a}}{\rm{n}}{\rm{d}}\,{{\rm{e}}}^{-}$ and neglect H₂ and HD chemistry.

We use an L = 1 cMpc h⁻¹ simulation box with periodic boundary conditions and run a number of DM-only simulations at low resolution to identify random seeds that produce halos that exceed 10⁷ ${M}_{\odot }$ at z > 10 with a variety of spin parameters $\lambda =J| E{| }^{1/2}/{{GM}}^{5/2}$, where J is the total angular momentum, E is the total energy (kinetic plus gravitational), and M is the total halo mass. The halos we chose have λ = 0.08, 0.005, 0.002, and 0.02 (labeled A, B, C, and D, respectively), which span the spin parameters found in cosmological simulations (e.g., Bullock et al. 2001). We considered this range of λ to investigate its impact on DCBH formation. In each simulation we center three static nested grids on the halo, a top grid and two additional grids enclosing the central 20% of the top grid with the same resolution (256³), and a number of DM particles. This setup yields an initial effective resolution of 1024³.

We allow up to 15 levels of refinement during the simulation for a maximum spatial resolution of ∼2000 au and minimum DM particle mass of ∼67 ${M}_{\odot }$. The grid is refined on baryon overdensity and DM particle mass, and we ensure that the Jeans length is resolved by at least 64 cells during the run, which has been found to be sufficient to resolve turbulent eddies in past cosmological simulations (Latif et al. 2013b). We employ a pressure floor to stabilize collapse on the smallest scales after reaching the maximum refinement level. This approach enables us to follow the evolution of clumps that could become SMSs out to 3 Myr. Further details on our simulation setup and refinement criteria can be found in Latif et al. (2016b) and Latif & Khochfar (2019).

Halos A, B, C, and D begin to atomically cool at masses of 5.16 × 10⁷ ${M}_{\odot }$, 1.8 × 10⁷ ${M}_{\odot }$, 2.1 × 10⁷ ${M}_{\odot }$, and 2.6 × 10⁷ ${M}_{\odot }$ at z = 10.8, 10.3, 12.7, and 13.3, respectively. We evolve each halo for 3 Myr after the onset of collapse (point when the maximum refinement level is reached in simulations with densities of ∼10⁻¹⁸ g cm⁻³) because this is enough time for an SMS to form and collapse to a DCBH at the infall rates we find in the accretion disks in our simulations, as we discuss below (see also Figure 4 of Woods et al. 2017). Spherically averaged profiles of all four halos at 3 Myr are shown in Figure 1, and projections of the accretion disks forming in them are shown at 0.5, 1.5, 2.5, and 3 Myr in Figure 2. They collapse nearly isothermally with central temperatures of ∼5000–8000 K. The halos have similar profiles at larger radii because of the self-similar nature of runaway gravitational collapse. Gas in the halos at radii greater than 10 pc exhibit Keplerian rotation with velocities of ∼15 km s⁻¹. At smaller radii (between 0.01 and 0.1 pc) these velocities rise sharply because of the formation of self-gravitating disks, as seen in Figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Gas densities vary from ∼10⁻²⁴ g cm⁻³ at the virial radii (several hundred pc) to ∼10⁻¹⁴ g cm⁻³ at the center, and have roughly the r⁻² profiles expected for isothermal collapse. The bumps in the density profiles are due either to massive clumps within disks, as in halos A and D, or to binary disks at a few parsecs, as in halos B and C. The enclosed gas mass increases sharply outward from the center and flattens out further out in the disks, which have masses of a few 10⁵ ${M}_{\odot }$. The disks in the high-spin halos A and D are more massive than the those in halos B and C. The enclosed gas mass within the virial radius reaches a few times 10⁶ M_⊙ in all four cases.

As shown in Figure 2, solitary disks form in halos A and D and partially fragment into a few satellite clumps, but halos B and C form binary accretion disks. Fragmentation is more frequent in the solitary disks, with clumps tidally stripping mass from each other at times (as shown, for example, at 2.5 Myr in halo A). At 1 Myr these fragments have typical masses of a few 10⁴ ${M}_{\odot }$. Most of them spiral into the massive clump at the center of the disk well before they could form stars, but a few are ejected into the surrounding medium and survive for 2 Myr. The accretion disks in halos B and C have masses of (2–3) × 10⁵ ${M}_{\odot }$ at 1 Myr. Although they also exert tidal forces on each other, they survive for 2 Myr, with average separations of 2 pc. The numbers and masses of the clumps in all four halos at 3 Myr are shown in Figure 3. The most massive ones reach a few 10⁵ ${M}_{\odot }$, with a few 10³–10⁴ M_⊙ fragments in the high-spin halos.

Zoom In Zoom Out Reset image size

We plot the masses of clumps that survive for more than 1 Myr along with their accretion rates in Figure 4. Not all the clumps forming in the high-spin halos appear because some already merged with the central object, while others only form in the last 0.6 Myr. The fluctuations in the accretion rates are due to fragmentation and highly turbulent flows in the disks. It is clear that the large inflow rates previously suspected (but never confirmed) to build up SMSs indeed persist down to small enough scales for long enough times to create such stars, and later DCBHs. The ratios of the masses of the two clumps in the binary systems is initially M₁/M₂ = 2/5, but the smaller ones grow to nearly the same masses as their partners by ∼1 Myr. Accretion rates in the binary disks average about 0.1 ${M}_{\odot }$ yr⁻¹ for two Myr, more than enough time for SMSs to form and collapse to DCBH binaries that could later merge into a single object.

Zoom In Zoom Out Reset image size

Even more stars will form in halo A, in which three massive clumps have grown for about 2 Myr at average rates of ∼0.1–0.2 ${M}_{\odot }$ yr⁻¹. One reaches a final mass of (4–5) × 10⁵ ${M}_{\odot }$, while the other two grow to (1–2) × 10⁵ ${M}_{\odot }$. In halo D all the clumps migrate inward and merge with the one at the center, likely forming just one SMS. It grows to 5 × 10⁵ ${M}_{\odot }$ by the end of the run. Although multiple clumps appear in this halo at later times, an SMS will already have formed at its center and collapsed by then, so it is unclear if they can become stars when exposed to X-rays from a DCBH. Most of the clumps lose mass at times in their evolution because of tidal stripping by other fragments. At such times the magnitudes of their accretion rates are plotted in the panel on the right since negative rates cannot appear on logarithmic scales.

Our simulations confirm for the first time that SMSs, and therefore DCBHs, can form in binaries or even small multiples in halos with low spins and high spins, respectively. They also demonstrate that flow rates in these halos continue for sufficient long times on small scales to create such objects. Although fragmentation has been reported in some recent studies they could not follow the evolution of the clumps to determine if they later formed stars or were simply subsumed into the center of the disk. Our simulations suggest that binary SMSs preferentially form in low-spin halos because of their lower angular momenta, while high-spin halos favor the formation of small multiples because the accretion disk breaks up more easily. We note that past numerical simulations exhibit fragmentation on much smaller au scales (Becerra et al. 2015) that are not resolved here. However, these fragments later merge with the central object on timescales of ∼10 yr and do not become stars themselves. Our failure to resolve them therefore does not alter the results of our study.

Although we do not include radiative feedback from SMSs in our simulations, it is not expected to have a large effect on the evolution of the clumps over time. Chon et al. (2018) examined the impact of UV feedback from SMSs in the unlikely scenario that they become blue, hot, and luminous in ionizing UV. They found that the stars created bipolar H ii regions in which the temperature of the gas was at most twice that of the surrounding gas. Mass loading by infall halts the expansion of the I-front and confines the ionized gas onto the disk around the star, with no effect on its growth. However, Smith et al. (2017) post-processed highly resolved simulations of atomically cooling halos with Lyα photon transport and found it might exert some mechanical feedback on flows in the vicinity of the star. Radiation hydrodynamical simulations by Luo et al. (2018) and Ardaneh et al. (2018) without resonant Lyα scattering found that radiation from the protostar in its early stages did not significantly alter flows in its vicinity but did suppress fragmentation, thus promoting the rapid growth of a single supermassive object, but they only evolved these systems for a few years and could not evaluate its effects at later times.

Our simulations have important consequences for the detection of DCBHs (and thus the first quasars) at birth. Inspirals of two or more DCBHs could emit gravitational-wave signals that are powerful enough to be detected at high redshifts by LISA. Likewise, if the environments of multiple DCBHs are dense with other fragments or less massive (and longer-lived) SMSs, they could produce tidal disruption events that would be extremely luminous in the near-IR today (Kashiyama & Inayoshi 2016). They could be found by the James Webb Space Telescope or extremely large telescopes in the coming decade.

M.A.L. thanks the UAEU for funding via UPAR grant No. 31S390 and startup grant No. 31S372. D.J.W. was supported by the Ida Pfeiffer Professorship at the Institute of Astrophysics at the University of Vienna and by STFC New Applicant grant ST/P000509/1.