THE FIRST GALAXIES: ASSEMBLY WITH BLACK HOLE FEEDBACK

To survey the relevant parameter space, we have carried out three cosmological simulations. As a reference, simulation "BHN," includes stellar radiative feedback from Pop III stars, whereas the subsequent feedback due to BH accretion is not taken into account. The simulation "BHS" includes both the feedback from Pop III stars and a single isolated BH. Finally, in simulation "BHB," we assume that the Pop III BH remnant has a stellar binary companion, giving rise to an HMXB. We provide a summary of the simulations in Table 1.

For the simulations presented here, we use a parallel version of the combined hydrodynamics and tree N-body code GADGET2 (Springel 2005). The code evaluates hydrodynamical forces using the SPH technique. Our simulations are initialized using a snapshot from the earlier simulation of Greif et al. (2010). We start running the simulations at z ∼ 30, corresponding to the point just before the first Pop III star in the computational box is formed, and terminate them at z ∼ 10, when the assembly of the first galaxy is expected to be complete. The original simulation was initialized at z = 100 in a periodic box of linear size of 1 Mpc (comoving), using ΛCDM cosmological parameters with matter density Ω_m = 1 − Ω_Λ = 0.3, baryon density Ω_b = 0.04, present-day Hubble expansion rate H₀ = 70 km s⁻¹ Mpc⁻¹, spectral index n_s = 1.0, and normalization σ₈ = 0.9, consistent with the WMAP5 measurements (Komatsu et al. 2009). Given that we use a higher value compared to the current measurement of σ₈ = 0.8 (Komatsu et al. 2011), structure formation may be accelerated in our simulations. However, we compensate to some extent for the lack of large-scale power in our simulations due to the limited box size by taking a value larger than the cosmological one. We also stress that in this paper we focus on individual halo properties which are not sensitive to the variation in the σ₈ parameter.

Using a standard zoom-in technique (see Greif et al. 2010), a preliminary run with 64³ particles was hierarchically refined to generate initial conditions with high–mass resolution inside the region destined to collapse into the first galaxy. Employing three consecutive levels of refinement, the mass of DM and gas particles in the highest resolution region is m_DM ∼33 M_☉ and m_sph  ∼  5 M_☉, respectively. The corresponding baryonic mass resolution is M_res ∼ 1.5N_neighm_sph ∼ 400 M_☉, where N_neigh ∼ 50 is the number of neighboring particles within the SPH smoothing kernel (Bate & Burkert 1997). The Jeans mass in primordial gas, where molecular hydrogen cooling imprints a characteristic density of n_H = 10⁴ cm⁻³ and a temperature of 200 K (Bromm et al. 2002), is thus marginally resolved.

We use the same chemistry and cooling network as in Greif et al. (2010), where all relevant cooling mechanisms, such as H and He collisional ionization, excitation and recombination cooling, bremsstrahlung, inverse Compton cooling, and collisional excitation cooling via H₂ and HD, are taken into account. For H₂ cooling, collisions with protons and electrons are explicitly included, in addition to the usually dominant neutral hydrogen atoms. The code self-consistently solves the rate equations for the abundances of H, H⁺, H⁻, H₂, H⁺₂, He, He⁺, He⁺⁺, and e⁻, as well as the three deuterium species D, D⁺, and HD. We can thus accommodate the non-equilibrium chemical evolution which is ubiquitous in early universe structure formation.

Employing the sink particle algorithm of Johnson & Bromm (2007), we convert an SPH particle into a collisionless sink particle if its hydrogen number density exceeds a threshold value of n_max = 10⁴ cm⁻³. When a sink particle forms, gas particles within an accretion radius, r_acc, are immediately accreted onto the sink. The position and velocity of a sink particle are estimated every time step based on a new mass-weighted position and velocity of the accreted gas particles within r_acc = L_res. Here, the resolution length of the simulation is defined as

$$\begin{equation} L_{\rm res} = 0.5 \left(\frac{M_{\rm res}}{\rho _{\rm max}} \right)^{1/3} \simeq 1 \, \mathrm{pc}, \end{equation} \tag{ 1 }$$

where ρ_max ≃ n_maxm_H.

A sink particle can grow in mass through the further accretion of surrounding gas. The criterion for subsequent accretion is that a neighboring SPH particle approaches the sink to within r_acc, which we hold constant throughout the simulation. Note that our simple prescription for sink growth is sufficient for the purpose of providing a marker for the position of a Pop III star and its BH remnant. When determining the mass growth of the BH which is responsible for the X-ray feedback, we use a more sophisticated methodology to determine its accretion rate (see Section 2.5).

After a Pop III star has formed inside its host minihalo, we turn on the radiation field emitted by this star, treated as a point source. We propagate the ionization front (I-front) around it to build up a primordial H ii region, using a well-tested ray-tracing algorithm (Greif et al. 2009). This scheme solves the I-front equation in a spherical grid by tracking 10⁵ rays with 500 logarithmically spaced radial bins around Pop III stars. The hydrodynamical effect is taken into account by self-consistently coupling the ray-tracing module to the chemical and thermal evolution of the gas.

We assume that a 100 M_☉ Pop III star emits blackbody radiation with an effective temperature T_eff = 10^4.9 K, and luminosity L_* = 10^6.1L_☉ (Bromm et al. 2001b; Schaerer 2002). The corresponding production rates for ionizing photons are: $\dot{N}_{\rm ion, {\rm H\,\mathsc{i}/He\,\mathsc{i}}}=9.1 \times 10^{49} \, \mathrm{s^{-1}}$ and $\dot{N}_{\rm ion, {\rm He\,\mathsc{ii}}}=4.1 \times 10^{48} \, \mathrm{s^{-1}}$. The growth of the I-front continues until it reaches its maximum size at the end of the Pop III star's life after t_* = 2.7 Myr. However, not all Pop III stars are endowed with an H ii region because of the computational expense incurred by the ray tracing. We trace the photons only from those stars formed within the Lagrangian volume destined to become the first galaxy at z ∼ 10, roughly corresponding to a 10 kpc radius from the center of this region.

Finally, we include the transfer of the H₂-dissociating Lyman–Werner (LW) photons in the range of 11.2–13.6 eV, isotropically emitted by a Pop III star as a 1/r² field without attenuation of the flux with radius, which propagate far beyond the H ii region. The photo-dissociating rate is $k_{\rm H_2} = 1.1 \times 10^8 F_{\rm LW}$, where F_LW is the radiation flux integrated over the LW bands (Abel et al. 1997). We compute F_LW in the optically thin limit and do not take H₂ self-shielding into account (see Greif et al. 2009, for more details). The neglect of self-shielding is justified due to the low H₂ column densities encountered in our simulations, of the order of ∼10¹³ cm⁻² inside the BH host galaxy (Draine & Bertoldi 1996; Wolcott-Green & Haiman 2011).

In an important difference from Greif et al. (2010), we here assume that all Pop III stars die via direct collapse into BHs, without any preceding SN explosion. For simplicity, however, we allow only one BH to produce radiative feedback due to the accretion onto it, taken as the first BH formed during the simulation. The heat input due to this BH feedback into the pre-galactic region is active for ∼350 Myr, while the radiation from individual Pop III stars is turned on for only 2.7 Myr. We next describe our treatment for estimating the accretion rate onto a BH and the model for the radiation emitted by an isolated BH and an HMXB.

2.5.1. Accretion Rate

We assume that the BH is embedded in a pressure-supported primordial gas cloud, steadily accreting from it. The corresponding accretion rate is given by the Bondi & Hoyle (1944) model, where a homogeneous medium which is at rest at infinity accretes onto a point mass. The Bondi–Hoyle rate can be written as

$$\begin{equation} \dot{M}_{\rm BH} = \frac{4\pi ({\it GM}_{\rm BH})^2 \rho _{\rm gas}}{\big(c_s^2+v_{\rm BH}^2\big)^{3/2}}, \end{equation} \tag{ 2 }$$

where the gas sound speed, c_s, and gas density, ρ_gas, are determined by averaging over the N_neigh ≃ 50 SPH particles closest to the BH. The relative speed of the BH with respect to the surrounding gas, v_BH, is negligibly small throughout the simulation. Varying N_neigh by an order of magnitude does not affect the accretion rates when they are estimated with Equation (2).

By assuming that 10% of the rest-mass energy of the accreted matter is released as radiation, we normalize the total luminosity according to $L = \epsilon \dot{M}_{\rm BH} c^2$, where = 0.1 is the radiative efficiency and c the speed of light. We estimate that the accretion rate $\dot{M}_{\rm BH}$ lies between 10⁻¹² M_☉ yr⁻¹ and 10⁻⁶ M_☉ yr⁻¹. The former corresponds to the initial situation, where the BH is located at the center of an H ii region with a temperature of ∼10⁴ K and hydrogen number density of ∼1  cm⁻³. The upper limit is derived from comparison simulations where stellar and BH radiative feedback is neglected.

It is important to note that the accretion rate used in this work is an upper limit for three reasons. In reality, when the feedback from radiation pressure is also taken into account, the true rates are likely much smaller than the nominal Bondi–Hoyle value (Milosavljević et al. 2009a, 2009b; Park & Ricotti 2012). We estimate, assuming a typical accretion rate of 10⁻⁸ M_☉ yr⁻¹, that the acceleration due to radiation pressure from Thomson scattering and photoionization in the vicinity of the BH, at r = 1 pc, is 30 times larger than the gravitational acceleration, a_grav = GM_BH/r² (e.g., Equation (10) in Johnson et al. 2011). This indicates that radiation pressure, mainly due to photoionization, is important near the BH, acting to limit growth. The other reason is that we have only considered the case of radiatively efficient accretion. If the cooling timescale for the liberation of viscously generated energy is longer than the accretion timescale, then most of the energy is advected inward with the accreting gas instead of being radiated away (Narayan & Yi 1994, 1995; Blandford & Begelman 1999, 2004). A fraction of the energy carried by the outflow could ultimately be converted to radiation, but the resulting spectrum would certainly be very different than that of the thin disk assumed here. A third reason is that the vorticity of the turbulent gas has not been properly resolved here. Studying its effect on Bondi accretion, Krumholz et al. (2005) found that even a small amount of vorticity is able to significantly reduce the accretion onto the BH. This suppression could be important in an atomic cooling halo, where the radial, cold flow along the filaments is converted into turbulent motion, thus generating vorticity in the central region of the galaxy.

In addition to the Bondi–Hoyle prescription, there is another way to estimate BH growth, based on $\dot{M}_{\rm sink}$ calculated with the method in Section 2.3, where the BH mass grows by accreting all gas particles within r_acc. The resulting rate is much higher than the Bondi–Hoyle value, roughly by a factor of 10. We adopt $\dot{M}_{\rm BH}$ rather than $\dot{M}_{\rm sink}$ as an estimator for BH growth, given that the sink accretion radius is about one order of magnitude larger than the Bondi–Hoyle radius, r_B = (μm_H GM_BH)/k_BT, where μ is the mean molecular weight and k_B is the Boltzmann constant, or, r_B ∼ 0.3 pc(M_BH/10² M_☉)(T/10² K)⁻¹. Moreover, only a fraction of the infalling gas will likely reach the vicinity of the BH. Furthermore, we can thus track the change in BH mass in a smooth manner and do not have to contend with discreteness effects owing to the limited mass resolution afforded by SPH.

2.5.2. Non-thermal Radiation

In the local universe, the emergent spectra from accreting BHs are typically modeled as a combination of a power-law component, F_ν∝ν^β, describing non-thermal (NT) synchrotron radiation, and a multi-color disk (MCD), resulting in a soft, thermal continuum (e.g., Mitsuda et al. 1984). In Figure 1, we show typical spectra to illustrate this model where a soft thermal component is fitted by an MCD blackbody with kT_in ∼ 150 eV (Miller et al. 2003), normalized assuming two representative accretion rates. Here, T_in is the temperature of the inner disk region which is related to the BH mass according to T_in∝M^−1/4_BH. It is currently not known whether accreting BHs at high redshift would behave in a similar manner.

Zoom In Zoom Out Reset image size

Conservatively, we here assume that there is no significant difference in the BH emission physics over cosmic time and use the same spectra. However, we ignore the MCD component of the emitted radiation in our simulations. Such a thermal contribution would determine the radiative feedback in the immediate "near zone" from the central BH, but is unimportant on larger scales, where the NT component dominates. We estimate the size of the near zone with a standard Strömgren analysis. Assuming that the ionizing photons originate from the MCD component, and employing a typical accretion luminosity of $\dot{M}_{\rm BH}=10^{-7}{\,M_\odot }$ yr⁻¹ and a typical density of n_H = 10  cm⁻³, we compute an equilibrium I-front radius of ≲ 5 pc, which is comparable to the spatial resolution of the simulations.

We follow Kuhlen & Madau (2005) to model the propagation of high-energy photons where the isotropic radiation field, ∝1/r², is a function of distance from the BH only. The coefficient for the ionization rate, k_ion, can be written as

$$\begin{equation} k_{\rm ion} = \int _{\nu _{\rm NT}}^{\nu _{\rm max}} \frac{F_{\nu } \sigma _{\nu }}{h\nu } d\nu, \end{equation} \tag{ 3 }$$

where ν_NT = 0.2 keV/h is the frequency where the NT component starts to dominate and ν_max = 10 keV/h is the upper cutoff frequency. Our choice of energy range, 0.2–10 keV, is motivated by the study of Miller et al. (2003), where the spectra of select ultraluminous X-ray sources (ULXs) associated with 100 M_☉ BHs were fitted with a power law in this same range. For the hydrogen and helium photoionization cross sections, σ_ν, we use the standard expressions (e.g., Barkana & Loeb 2001; Osterbrock & Ferland 2006).

To consider the secondary ionization effects from energetic electrons released through the absorption of an X-ray photon, we adopt the fitting formulae given by Shull & van Steenberg (1985; see also Valdés & Ferrara 2008; Furlanetto & Stoever 2010). They computed the fractions of the initial electron energy going into secondary ionizations of H i, secondary ionizations of He i, and into heating the surrounding gas. Secondary ionizations of He ii are not important (Shull & van Steenberg 1985) and are thus not included in the overall energy budget. The effective ionization rates can thus be written as

$$\begin{equation} k_{\rm eff}[\rm H\,\mathsc{ i}] = k_{\rm ion}[\rm H\,\mathsc{i}] + k_{\rm sec}[\rm H\,\mathsc{i}], \end{equation} \tag{ 4 }$$

where

$$\begin{equation} k_{\rm sec}[\rm H\,\mathsc{i}] = f_{\rm H}\left(\Gamma _{\rm H\,\mathsc{i}} + \frac{n_{\rm He\,\mathsc{i}}}{n_{\rm H\,\mathsc{i}}} \Gamma _{\rm He\,\mathsc{i}}\right)\frac{1}{13.6 \, \mathrm{eV}} \end{equation} \tag{ 5 }$$

and

$$\begin{equation} k_{\rm eff}[\rm He\,\mathsc{i}] = k_{\rm ion}[\rm He\,\mathsc{i}] + k_{\rm sec}[\rm He\,\mathsc{i}], \end{equation} \tag{ 6 }$$

where

$$\begin{equation} k_{\rm sec}[\rm He\,\mathsc{ i}] = f_{\rm He} \left(\Gamma _{\rm He\,\mathsc{i}} + \frac{n_{\rm H\,\mathsc{i}}}{n_{\rm He\,\mathsc{i}}} \Gamma _{\rm H\,\mathsc{i}}\right)\frac{1}{24.6 \, \mathrm{eV}}. \end{equation} \tag{ 7 }$$

Here, Γ, the rate at which the excess energy is released due to the first photoionization of the gas by X-rays, is given by

$$\begin{equation} \Gamma = n_{\rm n} \int _{\nu _{\rm NT}}^{\nu _{\rm max}} F_{\nu } \sigma _{\nu } \left(1 - \frac{\nu _{\rm min}}{\nu }\right) d\nu, \end{equation} \tag{ 8 }$$

where n_n is the number density of the respective unionized species, ν_min the ionization threshold frequency, and hν_min = 13.6 eV, hν_min = 24.6 eV, and hν_min = 54.4 eV for hydrogen, neutral helium, and singly ionized helium, respectively.

The fractions going into secondary ionizations, f_H and f_He, are a function of the hydrogen ionization fraction, x_ion = n_{H +}/n_H, and can be expressed as

$$\begin{equation} f_{\rm H} = 0.3908 \big(1-x_{\rm ion}^{0.4092} \big)^{1.7592} \hspace{2pt} \end{equation} \tag{ 9 }$$

$$\begin{equation} f_{\rm He} = 0.0554 \big(1-x_{\rm ion}^{0.4614} \big)^{1.6660}\mbox{\,} \end{equation} \tag{ 10 }$$

and the photoheating rates are then $\mathcal {H}=\Gamma (1-f_{\rm H/He})$.

The flux from the NT component is given by

$$\begin{equation} F_{\nu } = F_{\nu _0} \left(\frac{\nu }{\nu _{0}}\right)^{-\beta }, \quad \beta = 1, \end{equation} \tag{ 11 }$$

in the range of 0.2–10 keV, implying a mean hydrogen-ionizing photon energy of 390 eV. Here, $F_{\nu _0}$, the flux at ν₀ = ν_NT, is determined at every time step, based on the normalization

$$\begin{equation} \epsilon \dot{M}_{\rm BH} c^2 = \int _{\nu _{\rm NT}}^{\nu _{\rm max}} L_{\nu } d\nu, \end{equation} \tag{ 12 }$$

where L_ν = 4πr²F_ν is the specific luminosity. We note that since extinction along the line of sight is not taken into account, the estimated ionization and heating rates are likely upper limits. Recall that the accretion rate, $\dot{M}_{\rm BH}$, is updated at every time step, thus allowing our algorithm to reflect the changing conditions in the vicinity of the BH.

The resulting photoionization and photoheating rates are a function of distance from the BH, the instantaneous accretion rate, a function of the hydrogen ionization fraction, and can be conveniently written as

$$\begin{equation} k_{\rm ion} = \dot{K} \left(\frac{r}{\, \mathrm{pc}}\right)^{-2} \left(\frac{\dot{M}_{\rm BH}}{10^{-6} {\,M_\odot }\,{\rm yr^{-1}}}\right), \end{equation} \tag{ 13 }$$

where

$$\begin{equation} \dot{K} = [ 1.96, 2.48, 0.49 ] \times 10^{-11} \; \mathrm{s^{-1}}, \end{equation} \tag{ 14 }$$

and

$$\begin{equation} \mathcal {H} = n_{j} \hspace{2.84544pt} \dot{H} \left(\frac{r}{\, \mathrm{pc}}\right)^{-2} \left(\frac{\dot{M}_{\rm BH}}{10^{-6}{\,M_\odot }\,{\rm yr^{-1}}}\right) (1-f_{\rm H/He})\;, \end{equation} \tag{ 15 }$$

where

$$\begin{equation} \dot{H} = [ 7.81, 9.43, 1.63 ] \times 10^{-21} \; \mathrm{erg \hspace{2.27626pt} s^{-1}} \end{equation} \tag{ 16 }$$

for neutral hydrogen, neutral helium, and singly ionized helium, respectively. The corresponding number densities of each species j are given by n_j.

2.5.3. HMXB Emission

Recent improved simulations of Pop III star formation suggest that the primordial gas could fragment into two or more distinct cores, possibly resulting in a massive binary, or higher-multiple stellar system (Turk et al. 2009; Stacy et al. 2010, 2012; Clark et al. 2011a, 2011b; Greif et al. 2011; Smith et al. 2011; Prieto et al. 2011). This picture provides another possible X-ray source connected to the first stars, where the primary in a binary system leaves a BH behind at the end of its short life (Mirabel et al. 2011; Haiman 2011). The remnant could remain bound to the slightly less massive secondary, where nuclear burning is still going on. The material lost by the secondary in a wind or by Roche lobe overflow will be dumped onto the companion BH, releasing copious X-ray emission. The duration of such an HMXB is limited by the main-sequence lifetime of the donor and is of the order of 10⁷ yr.

It has been suggested that there is a correlation between the number of HMXBs ULXs, which can be explained as a subpopulation of the former, and low-metallicity environment (Majid et al. 2004; Dray 2006; Soria 2007; Mapelli et al. 2009; Linden et al. 2010). We can thus expect that the number densities of HMXBs and ULXs at high redshifts, when the universe was chemically pristine, were higher than what is observed locally (Mirabel et al. 2011; Haiman 2011). Such local surveys show that HMXBs and ULXs have luminosities of L_X ∼ 10³⁸ erg s⁻¹ and L_X ≳ 10⁴⁰ erg s⁻¹, respectively, which are close to the Eddington value, or mildly super-Eddington (Grimm et al. 2003). Consequently, we assume that the luminosity of an HMXB at high redshifts is close to the Eddington limit, L_Edd = 4π GM_BHm_Hc/σ_T, where m_H, σ_T, and c denote the proton mass, Thomson cross section, and the speed of light, respectively. We further assume that the spectrum exhibits the same power-law behavior as in Equation (11). The X-ray luminosity in the HMXB scenario is then determined from

$$\begin{equation} L_{\rm Edd} = 1.38 \times 10^{40} \, \mathrm{erg \hspace{2.84544pt} s^{-1}} \left(\frac{{M}_{\rm BH}}{100{\,M_\odot }}\right) = \int _{\nu _{\rm NT}}^{\nu _{\rm max}} L_{\nu } d\nu, \end{equation} \tag{ 17 }$$

and we assume M_BH = 100 M_☉.

The first Pop III star to appear in our cosmological box forms at a redshift z = 28 inside a 5 × 10⁵ M_☉ minihalo. Since the Pop III star is assumed to be very massive, with M_* ≃ 100 M_☉, it produces copious amounts of ionizing photons which carve out an extended H ii region into the surrounding IGM. The gas within the H ii region is photoheated to temperatures, T ∼ 3 × 10⁴ K, in excess of the virial temperature of the minihalo, thus triggering a hydrodynamic outflow. In Figure 2, we show the projected temperature approximately 1 Myr after the radiation from the central star was switched off. The terminal size of the H ii region is about ∼5 kpc,⁴ such that the I-front extends much farther out than the virial radius of the minihalo, r_vir ∼100 pc. Note that the anisotropic morphology of the H ii region reflects the inhomogeneous density distribution in the neighboring IGM. We describe the physical properties of the photoionized gas in Figure 3, again ∼1 Myr after the Pop III star died. At this time, the gas toward the center has already begun to cool and recombine, leading to a relic H ii region.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The expansion of the I-front begins with a very short lived R-type (supersonic) phase, followed by a D-type (subsonic) one, where the I-front is trapped behind a hydrodynamical shock wave (e.g., Whalen et al. 2004; Alvarez et al. 2006; Johnson et al. 2007). After ∼10⁵ yr, the I-front is able to break out, now again supersonically racing ahead of the shock into the low-density IGM. The resulting density profile, together with the location of the shock at ∼200 pc, and a corresponding velocity of v_sh ∼ 25–35 km s⁻¹ are in good agreement with the analytical champagne-flow solution (Shu et al. 2002). We note that contrary to the minihalos of ∼10⁶ M_☉ common at high redshifts, the extent of the shock radius is likely to decrease in rare higher-mass halos, with ≳ 10⁷ M_☉, due to their deeper potential wells and the effect of infalling gas, possibly trapping the I-front within the host halo, resulting in a compact H ii region (Kitayama et al. 2004; Yoshida et al. 2007a).

At the end of the life of the Pop III star after ∼2.7 Myr, when it leaves a massive BH remnant behind, the surrounding gas has been evacuated from the minihalo, thus preventing the relic BH from accreting any cold gas for an extended period of time. After the radiation from the central star shuts down, the relic H ii region starts cooling partly due to adiabatic expansion cooling, and partly due to atomic hydrogen line cooling, facilitated by collisional excitations from the enhanced abundance of free electrons. The elevated electron fraction catalyzes the formation of H₂ and HD molecules in the relic H ii region, thus providing additional molecular cooling (Ricotti et al. 2001; Johnson et al. 2007). Roughly 10 Myr after the radiation from the Pop III star turns off, the temperature decreases from 10⁴ K to 10³ K and the H₂ abundance increases from 10⁻⁸ to 10⁻⁴. While H₂ and HD molecules are susceptible to photodissociation from LW radiation produced in neighboring star-forming halos, the relatively short time required for the re-formation of H₂ and HD, compared to the average time between the formation of Pop III stars, allows their abundances to remain high.

Figure 4 shows the distance between newly formed Pop III stars and the central BH in the range of redshifts z = 10–20 from three simulations. A total of ∼50 Pop III stars have been formed in the BHN and BHS cases, accompanied by individual H ii regions, according to our criterion that stellar radiative feedback is taken into account only for those stars formed within a 10 kpc radius from the center of the emerging protogalaxy. We find that for over 250 Myr after the seed BH formed, there has been no further star formation within the host halo in simulation BHS. This is because it takes time for the gas expelled by the BH progenitor star to be reincorporated into the halo, and the modest feedback from the BH prevents the gas from cooling. In simulation BHN, on the other hand, stars inside the host halo continue to form already much earlier, demonstrating that X-ray heating from the accreting BH in simulation BHS implies a strong local negative feedback.

Zoom In Zoom Out Reset image size

Note the intense burst of star formation taking place within the virial radius of the protogalactic halo at z = 11.5–10.5 in the BHS. This starburst is fueled by a massive infall of gas. This gas is able to cool not only via molecular hydrogen cooling but also via atomic hydrogen line cooling. A fraction of the gas mass moving toward the central BH is thus consumed by star formation rather than being accreted onto the BH.

In Figure 5, we show the evolution of the stellar-mass density and the star formation rate density (SFRD). The comoving stellar-mass density is calculated using the stellar masses that have been formed within the 10 kpc radius from the center of the emerging protogalaxy at a given redshift. For the SFRD, we take the time derivative of the stellar-mass density at the redshift when the Pop III stars formed. For comparison, we overplot SFRDs using analytic fitting formulae (Hernquist & Springel 2003), for higher-mass halos where atomic hydrogen cooling is dominant, and for minihalos, assuming that one 100 M_☉ Pop III star forms per system via H₂ cooling (Yoshida et al. 2003).

Zoom In Zoom Out Reset image size

The estimated SFRDs from our simulations lie between the two analytical fits. Evidently, at z = 20 the SFRD in the BHB run is higher than in the BHN and BHS cases, by a factor of about three. This is a consequence of the positive feedback, where gas collapse into distant minihalos is facilitated via H₂ cooling promoted by the strong X-ray emission from the HMXB. As density fluctuations grow and star formation is enhanced, the local negative radiative feedback from Pop III stars begins to dominate, mitigating the positive HMXB feedback effect. We have not considered cooling due to metals and dust, both produced in SN explosions. At low densities, n < 500 cm⁻³, and at low metallicities, Z < 10⁻³ Z_☉; however, H₂ cooling is expected to dominate over metal-line cooling (Jappsen et al. 2007, 2009a, 2009b).

The characteristic property to distinguish the first galaxies from lower mass minihalos is their ability to cool through atomic hydrogen line emission, depending on the virial temperature of the hosting halo with a virial mass M_vir:

$$\begin{equation} T_{\rm vir} \sim 10^4 \, \mathrm{K} \left(\frac{M_{\rm vir}}{5\times 10^7{\,M_\odot }}\right)^{2/3} \left(\frac{1+z}{10}\right) \mbox{.} \end{equation} \tag{ 18 }$$

Above T_vir ≳ 10⁴ K, the gas within the halos is able to cool mainly via atomic hydrogen. In Figure 6, we show the evolution of the virial mass for the most massive halo, which will host the first galaxy at z ∼ 10, as well as of the total and cold gas masses for the three simulations presented in this work. The virial mass of the DM halo is estimated as the mass within a sphere with average DM density ρ ∼ 200 ρ₀(z), where ρ₀ is the mean cosmic density at a given redshift. Cold gas is identified by the condition that the gas temperature is less than half of the halo virial temperature, T < 0.5 T_vir(z). We also indicate the critical mass required for the onset of atomic hydrogen cooling at a given redshift.

Zoom In Zoom Out Reset image size

We find that at z ≳ 18, the halo is dominated by hot gas, exceeding the amount of cold gas by an order of magnitude. As time passes on, the cold gas mass increases, eventually accounting for ≳ 80% of the total gas mass in simulations BHN and BHS. This trend can be understood by the vulnerability of the halo gas to stellar radiative feedback. The corresponding evacuation of gas from the halo is very strong at high redshifts, z ≳ 18, because the halo potential wells were not yet deep enough to retain photoheated gas. We note that the sharp dips in the amount of cold gas are due to star formation inside the halo itself (at z ≲ 12), as well as in neighboring halos that are sufficiently close, within ≲ 5 kpc. The ionizing feedback from the accreting BH starts operating at z < 18, above which the effect is too weak as a result of the initially very low BH accretion rates.

While the total amount of gas is not sensitive to the BH feedback, the reduction in cold gas mass by a factor of ∼5 indicates that the additional heating from this feedback, on top of the stellar feedback, has a significant impact on the gas in the center of the forming galaxy. As the halo grows further via smooth accretion and mergers with minihalos, however, at z ∼ 13, both the total gas mass and mass of cold gas are no longer sensitive to the BH radiative feedback. At z ∼ 11.5, the condition for an atomic cooling halo is satisfied.

For simulation BHB, the heating from the HMXB is so strong that all gas particles have temperatures T > 0.5 T_vir(z), over the entire range of simulated redshifts z ≳ 15, and the total gas mass is reduced by an order of magnitude by photoevaporation, as is clearly evident in the dot-dashed lines shown in Figure 6. This result would imply that if an HMXB existed within a minihalo at high redshifts, it would take significantly longer for the halo to reassemble the lost gas and to eventually evolve into a primordial galaxy.

In Figure 7, we show the evolving conditions of the gas in the vicinity of the accreting BH in simulation BHS, while the first galaxy is assembled. Initially, the BH is located at the center of the H ii region produced by the progenitor star. As filamentary structure develops, gas infall and mergers with minihalos frequently occur along the filaments (see the middle panel of Figure 7). Cold accretion along the filaments, which are dense enough to allow molecule re-formation, in turn leading to enhanced cooling, efficiently delivers cold gas into the center of the halo. Roughly 300 Myr after the seed BH formed, the conditions for the onset of atomic hydrogen cooling are met. The virial radius of the halo at this time is r_vir ∼1 kpc.

Zoom In Zoom Out Reset image size

Figure 8 shows the evolution of the accretion rate onto the BH, the temperature and density of the neighboring gas, as well as the BH mass for simulations BHN and BHS. For the temperature and density, we average over the 50 closest SPH particles. The BH remnant, originating in the death of a Pop III star, initially resides within a medium of high temperature, ∼10⁴ K, and very low density, n_H ∼ 0.1 cm⁻³, preventing the gas from cooling and hence accreting onto the BH. At this stage, the accretion rate is only $\dot{M}_{\rm BH}\sim 10^{-11}{\,M_\odot }$ yr⁻¹, and the corresponding accretion luminosity, assuming a 10% radiative efficiency, is $L_{\rm acc}=\epsilon \dot{M}_{\rm BH} c^2 \sim 5.7 \times 10^{33}$ erg s⁻¹. This value is seven orders of magnitude below the Eddington luminosity, and BH growth and feedback are negligible in the beginning. As the potential well of the host halo becomes deeper, however, the amount of infalling gas increases with time, boosting the accretion rate, especially at z ∼ 17–19, as can be seen in Figure 8. Consequently, the accretion luminosity also becomes large enough to influence the surrounding gas, keeping it at a temperature that is an order of magnitude higher than in the BHN comparison simulation.

Zoom In Zoom Out Reset image size

The accretion rate in the BHN simulation is already comparable to the Eddington value, depicted as dashed lines in Figure 8 (top panel), at z ∼ 15.5, while it is still an order of magnitude lower in the BHS case. Occasionally, star formation takes place very close to the BH, e.g., ∼1 kpc away at z ∼ 14 (see Figure 4). The radiative feedback from this event acts to compound the heating effect from the BH accretion, thus rendering the removal of gas out of the shallow potential well more effective. Such gas evacuation at z ∼ 14 due to both feedback effects is clearly seen in Figure 6.

The combined stellar and BH radiative feedback results in an accretion rate that is on average four orders of magnitude below the Eddington value at z = 14–13. Even 300 Myr after the formation of the BH, its mass has increased by only 1.5% in the BHS simulation, whereas there is a two orders of magnitude growth in the BHN case. This indicates that the feedback from a stellar-mass BH is sufficiently strong to prevent significant growth, suggesting a very important constraint on SMBH formation scenarios. We infer that the radiative feedback from an accreting BH might be partly responsible for the low density of quasars at redshifts z ∼ 6, by suppressing early BH growth.

In Figure 9, we compare the temperature, electron fraction, H₂ fraction, and HD fraction of the primordial gas within the virial radius of the first galaxy halo as a function of number density for simulations BHN and BHS. We see that at low number densities, n_H = 10⁻³ − 1 cm⁻³, corresponding to the outskirts of the host halo, for both the BHN and BHS simulations the gas temperature is high, T ∼ 10³–10⁴ K, and the electron fraction is enhanced, $f_{e^{-}}=10^{-2}$, due to the accumulated relic H ii regions produced by previous star formation inside and outside of the halo and due to heating in accretion shocks.

Zoom In Zoom Out Reset image size

In the BHS run, additional heating from the BH allows the gas within the halo to maintain a relatively high temperature over the entire range of densities. There is a substantial difference seen in the electron fraction: while varying over the large range of $f_{e^{-}}=10^{-6}-1$ in the BHN simulation, in the BHS case most gas within the virial radius has electron fractions above ≳ 10⁻², and there is no gas present with fractions of <10⁻⁴. In response to the substantially increased electron fraction, the formation of H₂ and HD is catalyzed, leading to molecule abundances that are higher by two orders of magnitude compared to the BHN run.

In the absence of feedback from a BH in simulation BHN, the molecular hydrogen fraction converges to the well-known asymptotic value of $f_{\rm H_2}\sim 10^{-3}$. This behavior reflects a freeze-out process (Oh & Haiman 2002), where the molecular hydrogen abundance no longer evolves once the recombination and cooling timescales become less than the H₂ formation and dissociation timescales, t_rec, t_cool ≪ t_form, t_diss. Here, the recombination time is $t_{\rm rec}\sim 1/(\alpha _B n f_{e^-})$, and the H₂ formation and dissociation timescales are $t_{\rm form}\sim n_{\rm H_2}/\dot{n}_{\rm H_2,form}$ and $t_{\rm diss}\sim n_{\rm H_2}/\dot{n}_{\rm H_2,diss}$, respectively.

To illustrate the basic difference between the BHN and BHS simulations, we calculate a simple one-zone model, using the same chemistry module as in the simulations, but focusing on a single, representative fluid element as it falls into the center of the emerging galaxy. In the model, the radiation from the BH is imposed by a 1/r² radiation field, adopting a typical accretion rate of 10⁻⁸ M_☉ yr⁻¹, and we assume that the fluid element has a fixed density of n = 10 cm⁻³, placed r = 100 pc away from the BH at z = 11.6. Our fluid element is initially hot, T = 10⁴ K, and highly ionized, typical for a location inside a relic H ii region. The results of the one-zone model are shown in the top panel of Figure 10. To guide the interpretation of these results, we also show, in the bottom panel of Figure 10, the cooling time and the timescales for formation and dissociation of molecular hydrogen, and for the evolution of the electron fraction, $t_{\rm e}\equiv n_{e}/\dot{n}_{\rm e}$.

Zoom In Zoom Out Reset image size

In the case without BH feedback (solid lines), the H₂ formation and dissociation timescales are significantly longer than both the timescales for the evolution of the electron fraction and the cooling rate once the temperature has dropped below a few 10³ K, establishing the freeze-out value of $f_{\rm H_2}\sim 10^{-3}$. In the presence of BH feedback (dashed lines), on the other hand, the additional ionization from the BH X-ray emission significantly increases the fraction of free electrons at low temperatures. The enhanced electron fraction promotes the formation of H⁻, and this, at low temperatures, in turn implies an increase in the H₂ fraction. We note that similarly high H₂ fractions have been found in the related but different context of the non-equilibrium chemistry and cooling of primordial gas behind structure formation shocks under external irradiation (Shapiro & Kang 1987).

The overall distribution of the HD abundance is similar to that of H₂. In the BHS run, the HD fraction increases to f_HD ∼ 10⁻⁶, which exceeds the critical level needed for efficient cooling to the cosmic microwave background (CMB) temperature in local thermodynamic equilibrium, f_{HD, crit} ∼ 10⁻⁸, by two orders of magnitude (Johnson & Bromm 2006), as can be seen in Figure 9. Nevertheless, HD cooling does not succeed in tying the temperature to the CMB because of the continuous heating from the BH.

We plot the properties of the gas in simulation BHS in Figure 11, but this time as a function of distance from the central BH. The most noticeable feature here is the presence of the two distinct branches toward the center r <100 pc which is explained by the BH feedback. To better understand the BH feedback effect, in Figure 12 we follow the evolution of one representative SPH particle for ∼60 Myr. Initially, the tracked particle is located at 100 pc away from the central BH, where BH feedback is negligible. The sudden jumps in both temperature and electron fraction indicate a number of star formation events inside and outside of the halo. Soon after such events, the gas begins to cool back to a few ∼10³ K, as is typical for relic H ii regions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As the gas particle moves within 100 pc from the BH, however, the continuous heating from the accreting BH starts to operate, resulting in temperatures of ∼10³ K and an electron fraction of $f_{e^-}\sim 0.1$. As a consequence, H₂ molecule formation is catalyzed, leading to fractions of $f_{\rm H_2}\gtrsim 10^{-2}$. Finally, when the gas particle is located too close to the BH, r ≲ 10 pc, the strong X-ray flux heats up the gas again to ∼10⁴ K, leading to the dissociation of molecules and producing the two distinct branches within the central 100 pc, clearly seen in Figure 11. We point out that the dense clumps around 700–800 pc from the BH in Figure 11 correspond to the middle branch in the H₂ and HD molecule fractions, evident in the right panel of Figure 9.

In simulation BHB, where the continuous feedback from an HMXB is included, a dramatic change is seen in the gas properties surrounding the central HMXB. Figure 13 shows the comparison between simulations BHN and BHB at z = 15, plotted in the same fashion as in Figure 11, but now extending to 30 kpc from the HMXB to illustrate the effect of feedback not only on the gas within the host halo but also on the distant IGM. We carefully select the snapshots to guarantee that we are seeing the feedback effect from the HMXB without confusion from Pop III stellar radiation. In the BHB run, the gas temperature within the central 1 kpc reaches ∼10⁴ K, or even higher. The H₂ and HD molecule abundances are consequently suppressed via effective collisional dissociation, leading to values of $f_{\rm H_2}\sim 10^{-6}$ and f_HD ∼ 10⁻¹¹, respectively. The mean density of the gas in the central 100 pc is n_H ∼ 0.01 cm⁻³ which is four orders of magnitude smaller than that in the BHN run, showing the strong photoevaporation due to HMXB radiative feedback.

Zoom In Zoom Out Reset image size

In both simulations, part of the IGM outside a 1 kpc radius is heated and ionized by an I-front which had been generated 6.5 Myr ago by Pop III stars formed ∼5 kpc away from the center of the halo hosting the HMXB. The HMXB feedback effect on the IGM is most clearly evident in gas which has never experienced any stellar radiation feedback, as the latter could masquerade as the former. The distant IGM in the BHN simulation is in a cold and neutral state, while in the BHB run, the IGM unaffected by any stellar radiative feedback is in a warm and partially ionized phase because of the pervasive effect of the strong HMXB emission. These are favorable conditions for H₂ and HD molecule formation, promoting their abundance by an order of magnitude over the BHN case.

To more clearly distinguish the HMXB feedback on the first galaxy and on the distant IGM, we compare in Figure 14 the volume-averaged temperature, electron fraction, and H₂ fraction for simulations BHN and BHB inside and outside of a 15 kpc radius, beyond which our algorithm does not trigger the buildup of H ii regions around newly formed Pop III stars. Hence, any variations in gas properties in the far zone are caused by the HMXB feedback alone. Here, we adopt a 15 kpc cutoff radius rather than 10 kpc, the region inside which stars are endowed with an H ii region (see Section 2.4), because an ionized region generated by a Pop III star formed near r ∼ 10 kpc from the center of the emerging galaxy will extend out to ∼15 kpc. The corresponding X-ray heating increases the volume-averaged temperature within 15 kpc by only a factor of three, while the temperature in the distant IGM is enhanced by two orders of magnitude. The difference in the near zone is relatively small because the gas within the 15 kpc radius is affected by both the feedback from Pop III stars and from the HMXB, indicating that the overall gas properties are dominated by the local stellar feedback. In the absence of any local stellar feedback, as is the case in the region beyond the 15 kpc radius, on the other hand, the IGM is substantially pre-heated and partially ionized when HMXB feedback is present. These conditions facilitate molecule formation, enhancing abundances by an order of magnitude. Our results are in good agreement with the study of Kuhlen & Madau (2005), given that the X-ray flux in our simulation BHB is comparable to that for their power-law case.

Zoom In Zoom Out Reset image size

We assume that the HMXB source has been turned on for ∼150 Myr, which is much longer than any realistic lifetimes. For example, the main-sequence lifetime of a donor star of 30 M_☉, close to the predicted typical Pop III mass (Clark et al. 2011b; Hosokawa et al. 2011; Stacy et al. 2012), would be ∼20 Myr. Thus, our work can be considered to represent the case where at least one HMXB exists continuously over 150 Myr, somewhere in the Lagrangian volume of the emerging protogalaxy. A preheated IGM would provide us with a promising observational window into the early universe. Under these conditions, the IGM would be sufficiently hot to be decoupled from the CMB temperature, T_CMB = 2.73 K (1 + z), but not hot enough to be ionized, leaving the IGM substantially neutral. Such a preheated, neutral IGM is likely detectable in the emission of redshifted 21 cm radiation, which is produced by the spin-flip transition between the singlet and triplet hyperfine-structure levels of neutral hydrogen. In principle, future 21 cm observations could provide constraints on the character of early feedback processes, and in particular on the importance of a contribution from Pop III miniquasars (see, e.g., Furlanetto et al. 2006).

In the previous section, we discussed the situation established under radiative feedback from an accreting isolated BH and an HMXB. Effectively, we have made the assumption that all Pop III stars end their lives as a BH, without any SN explosion preceding their deaths. If a Pop III star, however, has a mass in the range of 140–260 M_☉, then it would die as a PISN, exerting a strong mechanical feedback on the surrounding primordial gas, enriching it with the metals produced before and during the explosion. We are able to directly compare our results to those of Greif et al. (2010), where PISN enrichment of the IGM, stellar radiative feedback, chemical mixing, and metal-line cooling were included, because all simulations start from the same initial conditions. The underlying DM structure was thus identical, providing us with an ideal laboratory to study the impact of variations in the baryonic physics on the assembly process of the first galaxies.

We find that the gas properties within the virial radius of the host halo in our BHN simulation at z ∼ 10 are very similar to those in the Greif et al. (2010) PISN simulation, except that the central gas in the latter simulation is enriched to Z ∼ 10⁻³ Z_☉. It is, however, still unclear whether such enrichment level is sufficient to enable the transition from Pop III to Pop II star formation, because the critical metallicity threshold for this transition is still a matter of debate. The additional heating from the central BH in our BHS run, on the other hand, significantly enhances molecule formation, possibly by up to two orders of magnitude, and the gas might be able to further fragment. Which feedback is dominant in shaping second-generation star formation in the first galaxies? The answer will strongly depend on the IMF of the first generation of stars and on the detailed physics of the transition from the predicted top-heavy IMF for Pop III to the more normal IMF in already metal-enriched gas.