Radiation Hydrodynamical Simulations of the First Quasars

TL;DR

This paper presents advanced radiation hydrodynamical simulations demonstrating that supermassive black holes observed at high redshifts can form from direct collapse black holes fed by cold accretion streams, matching observed properties.

Contribution

It introduces new cosmological simulations with X-ray coupling to primordial chemistry, explaining early SMBH formation consistent with recent high-redshift quasar observations.

Findings

Simulations reproduce observed SMBH mass, luminosity, and host galaxy properties.

Cold, dense accretion streams enable rapid SMBH growth in early universe.

Models align with measurements from ALMA and VISTA surveys.

Abstract

Supermassive black holes (SMBHs) are the central engines of luminous quasars and are found in most massive galaxies today. But the recent discoveries of ULAS J1120+0641, a $2 \times 10^9$ M$_{\odot}$ BH at $z =$ 7.1, and ULAS J1342+0928, a $8.0 \times 10^{8}$ M$_{\odot}$ BH at $z =$ 7.5, now push the era of quasar formation up to just 690 Myr after the Big Bang. Here we report new cosmological simulations of SMBHs with X-rays fully coupled to primordial chemistry and hydrodynamics that show that J1120 and J1342 can form from direct collapse black holes (DCBHs) if their growth is fed by cold, dense accretion streams, like those thought to fuel rapid star formation in some galaxies at later epochs. Our models reproduce all of the observed properties of J1120: its mass, luminosity, and H II region as well as star formation rates and metallicities in its host galaxy. They also reproduce the dynamical mass of the innermost 1.5 kpc of its emission region recently measured by ALMA and J-band magnitudes that are in good agreement with those found by the VISTA Hemisphere Survey.

Figures

Full text

1. Introduction

Nearly 160 quasars have now been discovered at $z>$ 6, including ULAS J1120+0641, a $2\times 10^{9}$ $M_{\odot}$ BH at $z=$ 7.1, and ULAS J1342+0928, an $8.0\times 10^{8}$ $M_{\odot}$ BH at $z=$ 7.5 (Mortlock et al., 2011; Bañados et al., 2018). Their origins and how they grew to such large masses by such early times are not yet clear, but they might have formed from either Population III star black holes (Pop III BHs) at $z\sim$ 25 or direct-collapse black holes (DCBHs) at $z\lesssim$ 20 (Kuhlen & Madau, 2005; Johnson et al., 2013). Pop III BHs are born in low densities that preclude rapid initial growth, so they would later require extended periods of hyper-Eddington accretion to reach 109 $M_{\odot}$ by $z\sim$ 7 to overcome the low duty cycles imposed by radiation from the BH (e.g., Whalen et al., 2004; Milosavljević et al., 2009a, b; Park & Ricotti, 2011, 2012, 2013; Madau et al., 2014; Volonteri et al., 2015). However, even if such modes of accretion could occur it is not clear that the BH would encounter enough gas to become an SMBH by $z\sim$ 7 (Alvarez et al., 2009; Smith et al., 2018). Many low-mass Pop III BH seeds would also be ejected from their birth halos, and thus their fuel supplies, by asymmetries in their core-collapse engines (Whalen & Fryer, 2012).

On the other hand, rapid baryon collapse in atomically-cooled halos (Lodato & Natarajan, 2006; Wise et al., 2008; Regan & Haehnelt, 2009; Shang et al., 2010; Latif et al., 2013) and the formation of supermassive stars that collapse to DCBHs (Hosokawa et al., 2013; Umeda et al., 2016; Woods et al., 2017; Haemmerlé et al., 2018a, b) require unusually strong sources of Lyman-Werner (LW) UV photons (Agarwal et al., 2012; Dijkstra et al., 2014; Yue et al., 2014; Johnson et al., 2014; Regan et al., 2017; Chon & Latif, 2017; Regan & Downes, 2018) or supersonic baryon streaming motions that may or may not occur in any given halo (Greif et al. 2011; Stacy et al. 2011; Latif et al. 2014; Hirano et al. 2017; Schauer et al. 2017 – see also Inayoshi & Omukai 2012; Inayoshi et al. 2015). But DCBHs can grow at much higher rates than Pop III BHs because they are born with much larger masses in much higher densities in host galaxies that can retain their fuel supply, even when heated by X-rays. A DCBH may now have been discovered in the Ly-$\alpha$ emitter CR7 at $z=$ 6.6 (Sobral et al., 2015; Hartwig et al., 2016; Smidt et al., 2016).

To reach 109 $M_{\odot}$ by $z\sim$ 7, SMBH seeds must also be fueled for long times, either by a succession of mergers with other gas-rich halos (e.g., Li et al., 2007) or by unusually strong cold accretion flows, like those thought to drive rapid star formation in some galaxies at later epochs (e.g., Dekel & Birnboim, 2006; Dekel et al., 2009; Bournaud et al., 2011). There may be up to a few dozen regions per Gpc*-3* with cold flows capable of forming quasars by $z\gtrsim$ 7 (Di Matteo et al. 2012; Feng et al. 2014; see also Costa et al. 2014; Hirschmann et al. 2014).

After birth, X-rays from the nascent BH regulate its growth by heating and ionizing flows in its environment and photoevaporating cold streams outside the halo. Winds, ionizing UV and supernovae (SNe) from stars later perturb flows onto the BH, either promoting or suppressing its growth (see, e.g., Dubois et al., 2015; Habouzit et al., 2017). Stars also govern their own rates of formation because their ionizing UV can trigger or quench star formation in nearby clouds (O’Shea et al., 2005; Whalen et al., 2008a, 2010). Metals from SNe also cool gas and create new stars (Kitayama & Yoshida, 2005; Whalen et al., 2008b; Ritter et al., 2012; Safranek-Shrader et al., 2014; Smith et al., 2015; Ritter et al., 2016; Sluder et al., 2016; Chen et al., 2017a, b).

Radiation from the BH can also promote or suppress star formation in its vicinity (Machacek et al., 2003; Aykutalp et al., 2014). X-rays produce energetic photoelectrons that cause secondary ionizations that enhance free electron fractions in primordial gas. Free electrons in turn catalyze the formation of H2 via the H*-* channel, cooling and collapsing gas and creating more stars (e.g., Glover & Abel, 2008). On the other hand, radiation from the BH can also evaporate star-forming clouds and global LW backgrounds can destroy H2, quenching star formation.

To model the formation of primordial quasars, simulations must resolve all these processes deep in the host galaxy of the BH and in cold streams in the intergalactic medium (IGM). We have now bridged these scales with new simulations that include X-rays from the BH and winds, ionizing UV and SNe from star-forming regions, and have evolved a quasar in cold flows from birth at $z=19$ down to $z=6$. We describe these models in Section 2 and examine the evolution of the BH and its host galaxy in Section 3, where we also compare its properties with those observed for J1120. We discuss our results and conclude in Section 4.

2. Numerical Method

We use the Enzo adaptive mesh refinement (AMR) cosmology code (Bryan et al., 2014) with the MORAY radiation transport package (Wise & Abel, 2011) to model the quasar in this study. X-ray and ionizing UV transport in MORAY includes radiation pressure on gas due to photoionizations and is self-consistently coupled to hydrodynamics and nine-species nonequilibrium primordial gas chemistry in Enzo. Secondary ionizations due to energetic photoelectrons and Compton heating by X-rays are taken into account in the chemistry and energy equations along with the usual primordial gas cooling processes: collisional excitational and ionizational cooling by H and He, recombinational cooling, bremsstrahlung cooling, H2 cooling, and inverse Compton cooling by the cosmic microwave background.

The BH is represented by a modified star particle whose X-ray luminosity, $L_{\mathrm{r}}$, is $\epsilon_{\mathrm{r}}\dot{m}_{\mathrm{BH}}c^{2}$, where $\epsilon_{\mathrm{r}}$ is the mean radiative efficiency, taken to be 0.1, and $\dot{m}_{\mathrm{BH}}$ is the accretion rate. This luminosity is all in the form of 1 keV photons to maximize the heating of gas because at higher energies they have smaller ionization cross sections and at lower energies they deposit less energy per ionization (Xu et al., 2014; Hummel et al., 2016). This energy is also consistent with new observations showing 90% of the X-ray flux of J1120 to be at 0.5 - 2 keV (Nanni et al., 2017). Our simulations do not resolve the accretion disk of the BH, which is $\sim$ 1 pc in diameter, so we use an alpha disk model to compute accretion rates to approximate angular momentum transport out of the disk on subgrid scales (DeBuhr et al., 2010). We also approximate a disk wind by depositing 10*-4* $L_{r}$ as thermal energy above and below the midplane of the BH perpendicular to its angular momentum vector (Ciotti et al., 2009).

We use a stochastic prescription for star formation that is based on Cen & Ostriker (1992) (hereafter CO92) but is slightly modified to accommodate a minimum star particle mass to avoid having too many low-mass particles (section 8.2.2 of Bryan et al., 2014). The algorithm tallies the mass accumulated in a region until it exceeds the minimum star particle mass, at which point it creates a star particle there if the conditions for SF in CO92 have also been met. This accumulated mass exceeds the minimum mass at random times so particles are formed sporadically. Furthermore, the particles will have a range of masses because the CO92 criteria might not be satisfied when the accumulated mass in a region exceeds the minimum mass, so it will grow in mass until they are met. We adopt a minimum star particle mass of 107 $M_{\odot}$.

Star particles are assigned a Salpeter initial mass function (IMF) for simplicity and are tagged as sources of ionizing UV photons, which are propagated throughout the simulation box by MORAY. Each particle emits photons at four energies: 12.6 eV (LW photons), 21.62 eV (the average energy of UV photons from massive, low-metallicity stars), 30 eV (He I ionizing photons), and 60 eV (He II ionizing photons). The relative numbers of photons apportioned to these bins are determined from the fits for the $Z=0$, no mass-loss tracks in Table 6 of Schaerer (2002). The particles also deposit momentum and metals from stellar winds into the ISM over their lifetimes. SN feedback is modeled as thermal energy, with 1051 erg deposited per explosion assuming one SN per 200 $M_{\odot}$ of stars formed. Cooling by metals from SNe is included with rates from Glover & Jappsen (2007), assuming solar yields that are consistent with our chosen IMF. Our models also have a uniform Lyman-Werner background due to early stars that evolves with redshift.

Our simulation box is 100 $h^{-1}$ Mpc on a side, with a 2563 root grid and three nested 25 $h^{-1}$ Mpc grids that are centered on the host halo for an effective resolution of 20483. These grids yield initial dark matter and baryon mass resolutions of 8.41 $\times$ 106 $h^{-1}$ $M_{\odot}$ and 1.57 $\times$ 106 $h^{-1}$ $M_{\odot}$, respectively. The grid is initialized with gaussian primordial density fluctuations at $z=$ 200 with MUSIC (Hahn & Abel, 2011) with cosmological parameters from the second-year Planck best fit lowP+lensing+BAO+JLA+H0: $\Omega_{\mathrm{M}}=$ 0.308, $\Omega_{\Lambda}=$ 0.691, $\Omega_{\mathrm{b}}=$ 0.0223, $h=$ 0.677, $\sigma_{8}=$ 0.816, and $n=$ 0.968 (Planck Collaboration et al., 2016). We use a maximum refinement level $l$ = 10 and refine the grid on baryon overdensities of 3 $\times$ 2*-0.2l* to obtain a maximum resolution of 35 pc (comoving). We also refine on a dark matter overdensity of 3 and resolve the local Jeans length with at least 32 zones at all times to avoid artificial fragmentation during collapse.

The box was chosen to be large enough to enclose the cold flows feeding the quasar on cosmological scales while resolving gas flows, photoionization and star formation deep within its host galaxy. However, given that there are only about a dozen regions per Gpc*-3* with reservoirs capable of sustaining rapid quasar growth in this manner, no single 100 Mpc box at random would be expected to enclose one. We therefore tested multiple random seeds at lower resolution until we obtained a halo that exceeded 1012 $M_{\odot}$ by $z\sim$ 7 and was the product of at most a few major mergers, i.e., with other halos similar in mass to the host halo (Trenti et al., 2008). The host halo of the BH is shown at $z=7.1$ in Fig. 1 to illustrate the morphology of the cold streams. We evolve the box from $z=$ 200 to 19.2, when the host halo reaches 3.0 $\times$ 108 $M_{\odot}$ and begins to atomically cool and collapse. At this redshift we initialize a 105 $M_{\odot}$ BH particle at the center of the halo and turn on X-rays and star formation.

3. SMBH / Galaxy Coevolution

BH masses and accretion rates are shown in Fig. 2 and star formation rates (SFRs) in the host galaxy and average accretion rates normalized to the Eddington limit are shown in Fig. 3. Actual accretion rates fluctuate rapidly between the limits shown in blue, with average values in red. Prior to the onset of star formation in the host galaxy, X-ray heating by the nascent BH limits average accretion rates to 0.8 $\dot{m}_{\mathrm{Edd}}$, where $\dot{m}_{\mathrm{Edd}}$ is the Eddington accretion limit, but this fraction soon falls to 0.3 by 400 Myr as the BH grows two decades in mass. X-rays regulate the growth of the BH even at these early times because cold streams carry dense gas deep into its host galaxy, raising its average density to 10 - 100 cm*-1*. At these densities 1 keV photons have mean free paths of 200 - 300 pc, less than the virial radius of the halo, so they efficiently heat the gas.

The host halo merges with a $1.2\times 10^{9}$ $M_{\odot}$ halo at 400 Myr and a $1.0\times 10^{10}$ $M_{\odot}$ halo at 450 Myr (mass ratios of 1:1 and 1:4, respectively). These collisions trigger the small starburst that begins at 420 Myr and peaks at 540 Myr. Disruption of the halo by the mergers and SN feedback from the burst perturb flows onto the BH, inducing brief episodes of super-Eddington accretion at 475 Myr and 510 Myr and causing the slight jump in BH mass at 450 Myr ($z=$ 10.7) in Fig. 2. The accretion ratio again begins to fall after this first starburst subsides.

After the two major mergers, the halo grows primarily by accretion from $5\times 10^{9}$ $M_{\odot}$ at $z\sim 11$ to $1.2\times 10^{12}$ $M_{\odot}$ by $z=7.1$. However, we note that while the halo mass does not exhibit any later jumps indicative of major mergers, smaller halos merge with the host halo over its entire evolution and are also carried into it by the cold streams themselves. After the first starburst, metals injected into the ISM by SNe trigger gas cooling and collapse that lead to a second, much larger starburst at 580 Myr that peaks at 770 $M_{\odot}$ yr*-1* ($z=$ 7.8). The absence of jumps in halo mass over this period indicates that the burst is triggered by metal cooling and primarily fueled by accretion.

The much stronger SN feedback from this burst and larger X-ray fluxes from the now more massive BH cause average accretion efficiencies to fall after 500 Myr, eventually to 0.2 $\dot{m}_{\mathrm{Edd}}$ by 900 Myr. Even so, the upper limit to the fluctuations in these rates is still at nearly the Eddington limit at $z=$ 7.1, consistent with observations of J1120. At no time do we artificially cap accretion rates onto the BH; except for brief episodes, X-rays and SN feedback limit them to below the Eddington rate throughout its life. The second, much larger burst in star formation is quenched by its own UV and SN feedback, falling to $\sim$ 150 $M_{\odot}$ yr*-1* by $z\sim$ 7.3. The SFR then recovers but levels off, never again exceeding $\sim$ 300 $M_{\odot}$ yr*-1* because of rising X-ray fluxes from the BH, even though cold streams continue to flow into the halo as shown in Fig. 1.

At $z=7.1$ the BH mass is $2.15\times 10^{9}$ $M_{\odot}$, within the error bars in mass for J1120. The SFR in the host galaxy is 245 $M_{\odot}$ yr*-1*, consistent with observations (105 - 340 $M_{\odot}$ yr*-1*; Venemans et al., 2017). The dynamical mass of the central 1.5 kpc of the host galaxy is $3.95\times 10^{10}$ $M_{\odot}$, in good agreement with recent measurements of (4.3 $\pm$ 0.9) $\times$ 1010 $M_{\odot}$ by the Atacama Large Millimeter Array (ALMA; Venemans et al., 2017). As shown in Fig. 4, metallicities are mildly supersolar ($\sim$ 2 $Z_{\odot}$) at the center of the host galaxy and approximately solar throughout the rest of its interior, as found in observations (Fig. 2 in Dunlop, 2013). Metals are also visible in nearby halos and in some of the accretion filaments threading the host galaxy. The spherically-averaged temperature profile of the H II region of the quasar at $z=7.1$ is shown in the right panel of Fig. 4. Both X-rays from the BH and SNe heat gas in the host galaxy to a few 106 K and the X-rays heat the surrounding IGM to temperatures of 104 - 105 K. Temperatures of 104 K extend out to $\sim$ 2 Mpc from the quasar, which is the observed radius of the ionized near zone of J1120. Using the CLOUDY code (Ferland et al., 2013) to compute the NIR luminosity of the quasar, we obtain a J band AB mag $=$ 20.6 at $z=$ 7.1, in good agreement with 20.2 from the VISTA High Latitude survey (VHS; Mortlock et al., 2011; Reed et al., 2017).

We show a census of star particles for the host galaxy of the quasar at $z=$ 7.1 in the left panel of Figure 5. The masses of the particles range from the minimum of 107 $M_{\odot}$ up to a few 109 $M_{\odot}$. While there are a large number of particles with masses of 107-7.5 $M_{\odot}$ most of the stellar mass of the galaxy is in the 108-9 $M_{\odot}$ particles, as can be seen in the right panel of Figure 5. The total mass of the stars is 3.1 $\times$ 1010 $M_{\odot}$, or 1.6% of the mass of the halo. This ratio is consistent with Figure 9 of Behroozi & Silk (2015), who predict that the stellar mass of a 1012 $M_{\odot}$ halo at $z=$ 7 should be 1 - 2% of the mass of the halo. There is a total of 2641 star particles in the host galaxy at this redshift.

SN feedback can be unphysically quenched in numerical simulations if the energy of the explosion is deposited as heat in star forming regions with high densities. Rapid cooling can dissipate this energy before it can launch flows that disrupt star formation or the growth of the BH. But as we show in Fig. 6, this never happens in our models: the ratio of the cooling time to the sound crossing time (equation 19 in Dalla Vecchia & Schaye, 2012) in regions of high metallicities and temperatures ($>10^{7}$ K) in the host galaxy is always greater than 1 and usually above 10. These ratios are due to the injection of momentum into the ISM by stellar winds and ionization heating by UV from stars, both of which clear gas away from star particles prior to explosion. Thermal energy from SNe is therefore deposited in lower densities that cannot cool quickly. X-rays from the BH also permeate the host galaxy in our simulations, heating gas and expanding it to lower densities.

4. Discussion and Conclusion

We adopted 1 keV photons for the BH spectrum in our models in part because it is the band in which most X-ray energy is found in high-$z$ quasars (Nanni et al., 2017). This choice produced good fits to all observed properties of J1120 at $z\sim$ 7 in our model. However, actual quasar spectra are expected to evolve from X-rays down to extreme UV over the large range in BH mass in our study. We show theoretical spectra for 105 - 109 $M_{\odot}$ BHs from Pacucci et al. (2015) in Fig. 7. The peak in spectral emission moves steadily towards longer wavelengths as the BH grows in mass, but the spectral energy distribution itself remains flat from 0.5 - 10 keV. These models neglect Compton upscattering by hot gas in the disk, which would redistribute photons at the spectral peak into the keV range (Done et al., 2012). The average photon energy of these spectra ranges from 2.0 keV at 105 $M_{\odot}$ to 0.3 keV at 109 $M_{\odot}$, so our choice of 1 keV is an appropriate average for this interval in BH mass.

That said, lower-energy photons in actual SMBH spectra may have a measurable effect on the evolution of the quasar because they have larger cross sections to ionization in H, He and some metal ions. But this effect appears to be modest at lower BH masses. Smidt et al. (2016) evolved a DCBH over two decades in mass with both monoenergetic X-ray and multigroup spectra that included ionizing UV and found little change in its evolution. Radiation pressure from IR on dust in the host galaxy of the BH may also drive strong but short-lived outflows (Bieri et al., 2017; Costa et al., 2018a, b), although such outflows have only been observed in a fraction of quasars at $z>$ 6. The softer components of SMBH spectra may also quench their images at 21 cm in emission because they ionize rather than heat the IGM. The effects of evolving BH spectra on the growth of the first quasars will be examined in future simulations.

While our models trace metal production due to star formation, they do not include dust formation and cooling or opacity due to dust and metals in the radiation transport. Dust cooling would promote star formation in environments in which radiation does not destroy it first, and both metals and dust would enhance X-ray heating in the galaxy by absorbing more photons. How these processes would offset each other and other factors governing star formation and the growth of the BH will be explored in future studies.

X-ray and SF feedback limit accretion to 0.2 - 0.8 $\dot{m}_{\mathrm{Edd}}$ in our models, rates at which quasars are not observed to have jets, so we do not include them in our models. Steady jets are observed in active galactic nuclei (AGNe) at $L\lesssim 0.01\,L_{\mathrm{Edd}}$ and intermittent jets are seen in quasars at $L\sim L_{\mathrm{Edd}}$ (Merloni & Heinz, 2008). In these cases the disk is geometrically thick and radiatively inefficient. No jets are seen at 0.01 $L_{\mathrm{Edd}}<L<L_{\mathrm{Edd}}$, when the disk is geometrically thin and radiatively efficient.

While our numerical simulations show that J1120 can be explained as a DCBH that formed at the intersection of cold accretion streams at $z\sim$ 20, they also reveal that J1342 could have formed in the same way. The BH reaches $1.18\times 10^{8}$ $M_{\odot}$ by $z=7.5$, slightly above the error bar in observed J1342 mass but well within the uncertainty of a factor of three associated with the use of the Mg II line as a proxy for its luminosity. As shown in the right panel of Fig. 2, we also obtain accretion rates of 27 $M_{\odot}$ yr*-1* at this redshift, in good agreement with observations.

Other scenarios could produce less massive BH seeds at high redshift. For example, if a 108 $M_{\odot}$ halo is marginally enriched by a SN the gas can fragment locally, forming a dense nuclear cluster of stars that build up a single massive star via runaway collisions (Devecchi & Volonteri, 2009; Reinoso et al., 2018; Boekholt et al., 2018). It then collapses to a BH that reaches 103 - 104 $M_{\odot}$ before the first SNe in the cluster begin to disrupt it (Whalen et al., 2008b). Some Pop III stars in less massive halos may reach masses up to 103 $M_{\odot}$ (Hirano et al., 2015). These processes could produce low-luminosity quasars at high redshift whose role in early reionization could be important, and they will be examined in future studies.

We thank Fabio Pacucci, Britton Smith, Marta Volonteri, John Wise and Hao Xu for valuable discussions. J. S., J. L. J. and H. L. were supported by a LANL LDRD Exploratory Research Grant 20170317ER. D. J. W. was supported by STFC New Applicant Grant ST/P000509/1 and the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007 - 2013) via the ERC Advanced Grant ”STARLIGHT: Formation of the First Stars” (project number 339177). Work at LANL was done under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy at Los Alamos National Laboratory under Contract No. DE-AC52-06NA25396. All Enzo models were performed under Institutional Computing (IC) allocations on Turquoise network platforms at LANL (Wolf).