THE FIRST GALAXIES: ASSEMBLY UNDER RADIATIVE FEEDBACK FROM THE FIRST STARS

We use a modified version of the N-body/TreePM SPH code gadget (Springel 2005; Springel et al. 2001b; Schaye et al. 2010) to perform a suite of zoomed cosmological hydrodynamical simulations of the assembly of a halo that reaches a virial mass M_vir ≈ 10⁹ M_☉ at redshift z = 10. The simulations are initialized at redshift z = 127 in a box of size 3.125 h⁻¹ Mpc comoving. Initial particle positions and velocities are obtained by applying the Zeldovich approximation (Zeldovich 1970) to particles arranged on a Cartesian grid. We adopt a transfer function for matter perturbations generated with cmbfast (version 4.1; Seljak & Zaldarriaga 1996).

We first perform a low-resolution simulation without star formation down to redshift z = 10, and locate the most massive halo, with virial mass M_vir ≈ 10⁹ M_☉ and virial radius r_vir = 3.1 kpc. We then trace the particles found within 3r_vir from the most bound particle of this halo back to their locations at the start of the simulation. We hierarchically refine the initial particle setup using a nested sequence of cubical patches ("zooms") centered on the traced particles, inside which we increase the mass resolution by successive factors of eight with each increasing level of zoom. In the zoom with the highest mass resolution, which entirely contains the traced particles and which we refer to as the refinement region, gas (dark matter) particles have masses m_g = 484 M_☉ (m_DM = 2350 M_☉). The simulations are performed with the gravitational forces softened over a sphere of Plummer-equivalent radius = 0.1 h⁻¹ kpc comoving applied to all particles.

The particle dynamical variables, such as position, velocity, and density, are evolved in time using individual particle gravito-hydrodynamical time steps determined by the smaller of the dynamical time step and the Courant time step (e.g., Equation (16) in Springel 2005), which is the standard gadget time stepping scheme. We do not explicitly limit the particle gravito-hydrodynamical time steps by either the chemical time or the radiative cooling or radiative heating time. However, chemistry, radiative cooling, and radiative heating described below are solved by subcycling the smallest gravito-hydrodynamical time step among all particles in the simulation on the relevant timescales (see Appendix A.1.4).

We assume that the gas is of primordial composition with hydrogen mass fraction X = 0.75 and a helium mass fraction Y = 1 − X. We use a modified version of the implicit solver dvode (Brown et al. 1989) to follow the non-equilibrium chemistry and cooling of H₂, D, HD, D⁺, H⁺, H, D, and He, and we include H⁻ and $\rm H_2^+$ assuming their collisional equilibrium abundances (Johnson & Bromm 2006; Greif et al. 2010). We consider all relevant radiative cooling processes: cooling by collisional ionization, collisional excitation of atomic and molecular lines, the emission of free–free and recombination radiation, and Compton cooling by the cosmic microwave background. Once stars form and emit radiation, the chemical and thermal evolution of the gas is also affected by the photodissociation of molecular hydrogen and deuterium, and photoionization of hydrogen and helium, as we describe in Sections 2.5 and 2.6. It should be kept in mind that at high gas densities, a Jeans floor employed to avoid artificial fragmentation artificially increases the gas temperature and affects the distribution and dynamics of the gas (see Equation (1) in Pawlik et al. 2011).

Star formation is a complex astrophysical phenomenon, many details of which remain to be understood (for a review see, e.g., McKee & Ostriker 2007). In the MW and in nearby galaxies, star formation is observed to occur inside a hierarchy of clouds with masses ∼10⁴–10⁶ M_☉. The clouds are transformed into stars at a rate $\dot{\rho }_\star = \rho _{\rm cl}/\tau _\star$, where τ_⋆ = τ_ff/_ff is the star formation timescale, τ_ff = [3π/(32Gρ_cl)]^1/2 is the free-fall time at the characteristic density ρ_cl of the star-forming clouds, and _ff is the star formation efficiency per free-fall time. Star formation is observed to be a slow process with a typical efficiency per free-fall time of only _ff ∼ 0.01, independent of the characteristic densities of the star-forming clouds in the range n_{H, cl} ≡ ρ_clX/m_H ∼ 10–10⁵ cm⁻³ (Krumholz & Tan 2007).

While gas masses ∼10⁴–10⁶ M_☉ that characterize star-forming clouds in the nearby universe are resolved, our simulations lack the resolution and physical detail to follow the formation of individual stars. Hence, we cannot exploit our simulations to estimate the SFRs inside individual star-forming clouds from first principles. Instead, we adopt a phenomenological model that specifies how quickly gas turns into stars. The model is motivated by and consistent with investigations of star formation in the nearby universe. Prescriptions for treating star formation in simulations of high-redshift galaxies are commonly calibrated with relations obtained from the local universe. This practice is necessitated by the current lack of direct observations of star formation at high redshifts. The phenomenological approach is sufficient to allow us to investigate the effect of stellar dissociating and photoionizing radiation on the interstellar gas and the IGM.

We restrict star formation to occur only in regions with gas densities exceeding a threshold density, n_H ⩾ n_SF, where we set n_SF ≡ 500 cm⁻³. This is somewhat lower than the densities ∼10⁴ cm⁻³ of metal-free clouds on the verge of collapse to form stars (e.g., Bromm et al. 2002; Abel et al. 2002). It is also lower than the similar densities of metal-free clouds able to shield their H₂ from external dissociating radiation (e.g., Safranek-Shrader et al. 2012). Because our simulations focus on the assembly of galaxies inside more massive halos, we must adopt, for reasons of computational viability, a lower resolution than employed in these previous works, preventing us from adopting still larger star formation threshold densities. Our choice for the star formation threshold density is close to the lowest threshold density n_H = 1000 cm⁻³ explored by Muratov et al. (2012) above which the time of the formation of the first star formed inside high-redshift minihalos was found to be insensitive to an increase in the threshold density (see also, e.g., Maio et al. 2009; Wise et al. 2012).

We adopt a star formation timescale that depends solely on the threshold density for star formation n_SF,

$$\begin{equation} \tau _\star \equiv \frac{\tau _{\rm ff} (n_{\rm SF})}{\epsilon _{\rm ff}} = 201.2 \,\mbox{Myr}\left(\frac{n_{\rm SF}}{500 \,\mbox{cm}^{-3}}\right)^{-1/2}, \end{equation} \tag{ 1 }$$

where we have set _ff = 0.01. Hence, the star formation timescale is independent of the gas density n_H ⩾ n_SF. This amounts to assuming that all star formation occurs inside clouds with characteristic density n_{H, cl} = n_SF, and is consistent with both our limited resolution and the observations in the nearby universe. A similar star formation recipe using a density-independent star formation time has been employed in, e.g., Kravtsov (2003). However, we caution that the SFRs of simulated galaxies depend on the specific choice for the value of the star formation efficiency, at least unless star formation is self-regulated by feedback (e.g., Ricotti et al. 2002; Haas et al. 2012).

Our numerical implementation of the star formation law is identical to that of Schaye & Dalla Vecchia (2008). The star formation law is interpreted stochastically, and the probability that a star-forming gas particle is turned into a star particle in a time interval Δt is given by min (Δt/τ_⋆, 1). Gas particles are converted to star particles assuming a conversion efficiency of 100%, i.e., the masses of the star particles are identical to those of the gas particles from which they are formed. The implied relatively large mass of star particles is consistent with our resolution, but sets a lower limit on the stellar mass fractions and hence ionizing and LW luminosities in the simulated galaxies to which the effects of radiative feedback may be sensitive. We impose an upper limit T_⋆ ≡ max (1.5 × 10⁴ K, 10^0.5 × T_floor) on the temperature at which gas is allowed to form stars, where T_floor is the minimum temperature set by the Jeans floor.

We interpret the star particles in our simulations as simple stellar populations, i.e., instantaneous stellar bursts that are characterized by an initial mass function (IMF), metallicity, and age. We compute the time-dependent hydrogen and helium ionizing luminosities $Q_{{\rm H\,\scriptsize{I}}}$, $Q_{{\rm He\,\scriptsize{I}}}$, and $Q_{{\rm He\,\scriptsize{II}}}$, as well as the luminosities Q_LW in the LW band with energies 11.2–13.6 eV, of these star formation bursts using the population synthesis models from Schaerer (2003). The models assume a power-law IMF $p(m_\star) \propto m_\star ^{-\alpha }$ with the Salpeter (1955) exponent α = 2.35 but allow for a variation of the range of the stellar masses sampled from the power-law distribution. We describe the stellar bursts using the Schaerer (2003) zero-metallicity models with initial masses in the range 50–500 M_☉. The age of a burst is the time difference between the simulation time at which the star particle was created and the current simulation time. For reference, the luminosities at zero-age main sequence are $Q_{{\rm H\,\scriptsize{I}}} = 10^{47.98}$, $Q_{{\rm He\,\scriptsize{I}}} = 10^{47.80}$, $Q_{{\rm He\,\scriptsize{II}}} = 10^{47.05}$, and $Q_{\rm LW} = 10^{46.96} \,\hbox{photons} \,M_{\odot }^{-1}\,\mbox{s}^{-1}$. We only consider star particles inside the refinement region, and we do not follow the propagation of radiation outside this region.

Molecular hydrogen H₂ and deuterated hydrogen HD are photodissociated upon absorption of radiation in the LW band belonging to the photon energy range 11.2–13.6 eV. The rate of photodissociation of H₂ is (e.g., Abel et al. 1997)

$$\begin{eqnarray} k_{\rm LW} &=& 1.1 \times 10^8 \,\mbox{s}^{-1}\frac{F_\nu (\nu _{\rm LW})}{\,\mbox{erg}\,\mbox{Hz}^{-1}\,\mbox{s}^{-1}\,\mbox{cm}^{-2}}\nonumber \\ &=& 1.38\times 10^{-12} \,\mbox{s}^{-1}J_{21}, \end{eqnarray} \tag{ 2 }$$

where ν_LW ≡ 12.4 eV is the characteristic LW frequency, F_ν(ν_LW) is the LW flux, J₂₁ = J_ν/(10⁻²¹ erg Hz⁻¹ s⁻¹ cm⁻² sr⁻¹) is the normalized LW mean intensity, and J_ν = F_ν(ν_LW)/4π. We set the rate for photodissociation of HD identical to that of H₂ (e.g., Glover & Jappsen 2007; Wolcott-Green & Haiman 2011).

We write the LW intensity J₂₁ as the sum of the cosmological LW background intensity $\bar{J}_{21}$ and the local LW intensity $J_{21}^{\rm loc}$ produced by the stellar populations represented by the star particles in the refinement region of our simulations, $J_{21} = \bar{J}_{21} + J_{21}^{\rm loc}$. The refinement region is too small to follow the buildup of the cosmological LW background in a self-consistent manner. We therefore treat the intensity of the LW background as a free parameter, and approximate its evolution using

$$\begin{equation} \bar{J}_{21}(z) = \bar{J}_{21,0} \times 10^{-(z - z_0)/5}, \end{equation} \tag{ 3 }$$

where $\bar{J}_{21,0}$ is the intensity of the LW background at z₀. Setting $\bar{J}_{21,0}= 1$ and z₀ = 10, Equation (3) provides a good fit to the LW background evolution presented in Figure 1 of Greif & Bromm (2006), and is also consistent with computations of the LW background in other works (e.g., Wise & Abel 2005; Ahn et al. 2009). The use of an LW background that is a fixed function of redshift does not allow us to capture the effects of self-regulation of star formation inside minihalos by LW feedback (e.g., Ahn et al. 2012), at least until the intensity of LW radiation emitted by the stars inside the simulation box has become higher than the intensity of the background.

For reasons of computational efficiency, we determine the contributions from the N_⋆ star particles in the refinement region to the intensity $J_{i}^{\rm loc} = \sum _{j = 1}^{N_\star } J^\star _{j}$ evaluated at the location of gas particle i in the optically thin approximation (e.g., Wise & Abel 2008a), i.e.,

$$\begin{eqnarray} J_{j}^{\star } &=& \frac{h_{\rm P} \nu _{\rm LW}}{4 \pi \delta \nu _{\rm LW}} \frac{Q_{{\rm LW}, j} m^\star _j}{4 \pi r_{ij}^2} , \end{eqnarray} \tag{ 4 }$$

or, in units of 10⁻²¹ erg Hz⁻¹ s⁻¹ cm⁻² sr⁻¹,

$$\begin{eqnarray} J^\star _{21, j}&\approx & 1 \left(\frac{Q_{{\rm LW}, j}}{10^{47}\,M_{\odot }^{-1}\,\mbox{s}^{-1}}\right) \left(\frac{m^\star _j}{500 \,M_{\odot }}\right) \left(\frac{r_{ij}}{1 \,\mbox{kpc}}\right)^{-2}. \end{eqnarray} \tag{ 5 }$$

Here, Q_{LW, j} is the photon luminosity of star particle j per unit mass in the LW frequency band, δν_LW = (13.6 eV–11.2 eV)/h_P is the width of the LW band, h_P is Planck's constant, $m^\star _j$ is the mass of star particle j, and r_ij is the distance between the gas and the star particle.

For simplicity, we approximate the time-dependent LW luminosities Q_LW computed in Section 2.4 by their zero-age main-sequence values, and we assume that the stars emit LW radiation for 3 Myr (Schaerer 2003). We approximately account for radiative transfer (RT) effects by attenuating the total LW intensity J₂₁ by a local self-shielding factor (Equation (10) in Wolcott-Green et al. 2011 with α = 1.1). The self-shielding factor depends on the column density of molecular hydrogen, which we compute in the local Jeans approximation (e.g., Shang et al. 2010).

When computing the contribution to the LW intensity from local sources using Equation (4), we ignore the absorption of LW radiation along the way to the absorbing gas particle. Near emitters, this is justified because molecular hydrogen is efficiently collisionally dissociated inside the H ii regions surrounding them (Johnson et al. 2007; but see Ricotti et al. 2002). Outside the H ii regions, this approximation is good as long as the intergalactic molecular hydrogen fraction does not substantially exceed its primordial value, $\eta _{\rm H_2} \sim 10^{-6}$. In this case, the optical depth for absorption in the LW band remains insignificant out to distances ≳ 1 Mpc comoving (e.g., Ricotti et al. 2001, their Figure 12; see also Glover & Brand 2003, their Figure 7; Ahn et al. 2009), much larger than the size of the refinement region. Finally, we also ignore the absorption of LW photons by atomic hydrogen series lines, which also is not significant on the small scales simulated here (e.g., Haiman et al. 2000; Ahn et al. 2009). Our implementation of photodissociation is similar to that in Johnson et al. (2013), but we do not account for the photo-detachment of H⁻.

In this section, we describe our implementation of photoionization by stellar radiation. This implementation involves two main steps. First, we transport ionizing photons radially from the star particles through the simulation box, and compute the fraction of the photons absorbed by atomic hydrogen and helium. Second, we infer the associated photoionization and photoheating rates. Here we will only give a brief overview. We present a detailed description of the implementation of the RT and the computation of the photoionization and photoheating rates in the Appendix. There we will also discuss tests of this implementation. The coupling of the RT with the hydrodynamical evolution is achieved by passing the photoionization and photoheating rates, along with the photodissociation rates described in the previous section, to the non-equilibrium solver for the chemical and thermal evolution of the gas described in Section 2.2, and this coupling is described in Appendix A.1.4.

We transport the ionizing radiation emitted by the star particles using the multi-frequency RT code traphic (Pawlik & Schaye 2008, 2011). traphic solves the time-dependent RT equation by tracing photon packets emitted by source particles at the speed of light and in a photon-conserving manner through the simulation box. The photon packets are transported directly on the spatially adaptive, unstructured grid traced out by the SPH particles, which allows one to exploit the full dynamic range of the SPH simulations. A directed radial transport of the photon packets from the sources is accomplished despite the irregular distribution of SPH particles by guiding the photon packets inside cones. A photon packet merging technique renders the computational cost of the RT independent of the number of ionizing sources. The transport of photons is discretized in RT time steps Δt_r, after each of which the chemical and thermal evolution of the gas is advanced based on the number of absorbed photons.

We make the following approximations specific to the current work. Each star particle emits ionizing photon packets to its neighboring SPH particles once per RT time step in a set of tessellating emission cones centered around eight different directions. The effective angular sampling of the surrounding volume is larger than implied by this number of directions thanks to the splitting of photon packets among neighbors inside the same emission cone, and because of the randomization of the emission directions at each RT time step. The photons are transported radially away from the star particles by tracing them downstream inside transmission cones with solid angle 4π/128. This angular resolution is sufficiently high to track the delay of the ionization fronts around individual halos by dense filaments, giving rise to the typical "butterfly" shape of ionized regions, as shown in Figure 1.

Zoom In Zoom Out Reset image size

To reduce the computational cost, we limit the propagation of photons to at most a single inter-particle distance per RT time step, which approximates the full time-dependent RT in the limit of small RT time steps. Ionizing photons are transported using a single frequency bin, and the absorption of the photons by neutral hydrogen and neutral and singly ionized helium is computed in the gray approximation. The gray approximation does not allow us to capture effects of spectral hardening, and hence we may underestimate the effects of photoionization and photoheating near and ahead of ionization fronts, such as the enhanced formation of molecular hydrogen (Ricotti et al. 2001, 2002). These approximations are discussed in further detail in the Appendix, in particular in Appendix A.2.3.

Figure 2 shows the formation history of the dwarf galaxy in simulation LW+RT, including both photodissociating and photoionizing radiation (blue curves). For comparison, the figure also shows the formation histories of the dwarf galaxy in the simulation in which star particles were sources of LW but not of photoionizing radiation (LW; red curves), and in the simulation in which star particles remained dark and the intensity of the LW background was set to zero (NOFB; black curves). The formation history is obtained by using subfind to locate the dwarf halo progenitor that contains most of the 50 most bound particles of the dwarf halo at the final simulation redshift z = 11, and then repeating this procedure to find the progenitor of this progenitor and so on, tracing the halo assembly back to z = 25. The quantities displayed in Figure 2 are obtained from the properties of the particles inside the virial radius. The figure shows that the dark matter halo hosting the emerging dwarf galaxy grows by about three orders of magnitude in mass in about 300 million years, consistent with expectations (Pawlik et al. 2011, their Figure 1).

Zoom In Zoom Out Reset image size

The evolution of the galaxy in simulation LW+RT proceeds in several main phases which will be discussed below: (1) the assembly of a dark matter minihalo with mass ∼10⁶ M_☉ and the accretion and condensation of gas inside it, leading to the formation of stars just below z = 23, (2) the subsequent accretion of gas under feedback from star formation inside the dark matter minihalo, (3) the evolution of this minihalo into an atomically cooling halo at z ≈ 16, during which the properties of the galaxy become insensitive to the inclusion of LW radiation, (4) the ensuing growth into a dwarf halo, during which the properties of the galaxy become robust against feedback from photoionization, and, finally, (5) the formation of two nested rotationally supported gaseous disks below z ≲ 13.5. Our discussion of these phases will be complemented by comparisons with the evolution of the dwarf galaxy in simulations LW and NOFB.

During the first phase, and as the dark matter halo grows in mass and accretes gas, both the baryon fraction, which is initially slightly smaller than the cosmic baryon fraction Ω_b/Ω_m ≈ 0.17, and the central gas densities increase. By z = 25, the molecular hydrogen fraction has significantly departed from its initial value, enabling the minihalo gas to cool efficiently. The average gas temperature inside the virial radius is initially slightly higher than the virial temperature (dotted curve). Note that in simulation NOFB, as the molecular hydrogen fraction builds up, this relation then reverses and the virial temperature becomes higher than the average gas temperature (e.g., O'shea & Norman 2007). At z ≈ 23, the central gas densities become larger than the threshold density for star formation. A single gas particle is turned into a star particle, triggering the emission of LW and ionizing radiation from the associated stellar burst. At that time, the halo has reached a mass of ∼2 × 10⁶ M_☉.

Photoionization from the radiation emitted by the first stellar burst almost instantly increases the average gas temperatures to ≳ 10⁴ K. The associated increase in thermal pressure pushes the gas away from the minihalo center, and removes a fraction of it from inside the virial radius. As a consequence, the baryon fraction is reduced to about f_bar = 0.15. The reduction in the baryon fraction is consistent with but slightly smaller than that found in previous simulations of the assembly of minihalos under feedback from star formation (e.g., Wise & Abel 2008a; Wise et al. 2012). This may be because the minihalo simulated here is fed by dense filaments, and the inflow of gas along these filaments provides a strong obstacle for photoheating to drive the gas beyond the virial radius, and it replenishes the photoevaporated regions with fresh gas (e.g., Abel et al. 2007). Nevertheless, the reduction in the central gas mass is sufficient to drive the central gas densities below the threshold density for star formation, and the stellar burst shuts itself off. Figure 1 shows images of the gas density, temperature, and ionized and molecular hydrogen fraction around the minihalo at the end of the stellar burst.

After the first stellar burst is shut off, the average mass-weighted molecular hydrogen fraction approaches ∼2 × 10⁻³, as expected inside the relic H ii region (Oh & Haiman 2002). The increased fraction of molecular hydrogen enables the gas to cool quickly, and the central gas densities increase. A few tens of Myr after the end of the first burst, the gas has become sufficiently cold and dense for another stellar burst to be ignited (e.g., O'shea et al. 2005; Alvarez et al. 2006). Photoionization heating from this second burst again lowers the gas densities, but this time there is no strong decrease in the baryon fraction. The accreting minihalo is thus massive enough to retain most of the gas inside the virial region. Because the negative feedback from photoheating is now less strong, this second starburst is more extended in time than the first one, and it involves the conversion of several gas particles to star particles. However, the combined feedback from the stellar clusters represented by the star particles eventually shuts off star formation. But already a few tens of Myr later the central gas densities have again increased above the SF threshold density, and the galaxy continues to form stars.

By redshift z ≈ 16, the halo has reached virial temperatures T_vir ≳ 10⁴ K. Consequently, a significant fraction of the atomic hydrogen is collisionally excited, and its radiative de-excitation endows the gas with an additional channel to lose its thermal energy. Feedback from photoionization heating continues to keep the gas inside the halo at relatively low densities. However, the galaxy now forms stars continuously, albeit at a rate significantly smaller than in the absence of radiative feedback. The comparison with the quickly rising SFRs in simulation LW that included emission of LW radiation but not that of ionizing photons shows that the feedback from the photodissociation of molecular hydrogen by LW radiation alone becomes inefficient in preventing the gas from forming stars as the halo mass approaches the atomic cooling limit, in good agreement with previous works (e.g., Ahn & Shapiro 2007; Wise & Abel 2007a; O'shea & Norman 2008). The transformation of the minihalo into an atomic cooling halo is accompanied by a significant decrease of the baryon fraction by ≈20% at z ≈ 15–16 and in the specific angular momentum of the gas j_gas ≡ J_gas/M_gas, where J_gas is the total gas angular momentum with respect to the motion of the most bound particle and M_gas the total gas mass.

Below z ≲ 15, photoionization heating becomes increasingly inefficient at reducing the density of the halo gas, and the gas then condenses up to ≳ 10⁵ cm⁻³. As a result, the expansion of H ii regions is impeded by the increased recombination rates. The properties of the galaxy then become rapidly insensitive to the inclusion of ionizing radiation. The average mass-weighted ionized fraction, which has been increasing until then, remains roughly constant at around $\langle \eta _{{\rm H\,\scriptsize{II}}} \rangle \lesssim 0.1$. As a result of the large densities, the SFRs approach those seen in the simulation without radiation. The average mass-weighted molecular hydrogen fraction reaches values in excess of $\langle \eta _{\rm H_2} \rangle \gtrsim 2 \times 10^{-3}$ as H₂ forms efficiently not only in fossil H ii regions (e.g., Oh & Haiman 2002; Johnson et al. 2007) but also in the ionization fronts of the H ii regions under irradiation from the ionizing stars (Ricotti et al. 2002; see also Shapiro & Kang 1987; Jeon et al. 2012).

A major merger at redshift z ≈ 15 is accompanied by a strong increase in the specific angular momentum of the halo gas. The gas then settles in a disk, which is essentially in place at z ≈ 13.5. The merger and the formation of this disk are shown in Figure 4, which display a sequence of snapshots of the gas density inside the forming dwarf galaxy. The disk is surrounded by a second larger-scale disk at z ≈ 12. The evolution and the properties of these two disks will be discussed in more detail in Section 4.3. After the formation of the disks, the difference between virial and average gas temperatures that has amplified since the halo has become an atomic cooling halo continues to increase. In the atomic cooling halo this difference arises primarily because a fraction of the gas does not shock-heat to the virial temperature upon accretion onto the halo. Instead, the gas can efficiently cool inside the dense filaments that reach into the virial region and supply the halo center with gas (Pawlik et al. 2011). A qualitatively similar difference between average and virial temperatures has been seen in other simulations of high-redshift low-mass galaxies (e.g., O'shea & Norman 2007; Wise & Abel 2007b; Greif et al. 2008). At z = 11, the final simulation redshift, the SFRs reach ∼0.2 M_☉ yr⁻¹, consistent with the SFRs found in previous simulations of high-redshift low-mass galaxies (e.g., Wise & Cen 2009; Razoumov & Sommer-Larsen 2010; Yajima et al. 2011; Wise et al. 2012).

Radiative feedback also affects the properties of the dark matter halo hosting the emerging dwarf galaxy. Figure 3 compares the evolution of the dark matter density profile of the galaxy in simulation NOFB with that in simulation LW+RT (left panel), and with that in simulation LW (right panel). The bottom panels show the dark matter density profiles for each pair of simulations at four representative redshifts. The profiles are scaled by dividing by a singular isothermal profile $\rho _{\rm iso} (r) = \rho _{\rm vir} r^2_{\rm vir} / r^2$. In the last expression, r_vir is the virial radius of the halo, ρ_vir = 200ρ_crit, and ρ_crit is the critical density. The top panels show the ratios of the density profiles of the simulations including radiation and the simulation without radiation shown in the bottom panels, which helps to illustrate the effects of stellar feedback. In simulation NOFB, i.e., in the absence of radiative feedback (dashed curves), the dark matter density profile is approximately singular isothermal at all redshifts.

Zoom In Zoom Out Reset image size

The central dark matter densities in the simulations that included ionizing and/or LW radiation (solid curves in each of the bottom panels) are initially significantly lower, by up to factors ∼5, than those in the simulation without radiation. The dark matter density profiles in these simulations thus do not follow a singular isothermal shape but show a spatially resolved dark matter "core," extending to radii significantly larger than the gravitational softening scale (leftmost vertical lines). The reduction in the central dark matter densities is more distinct and exists down to lower redshifts in simulation LW+RT than in simulation LW. The difference in the central dark matter densities originates in the difference in the distribution of the gas inside the assembling dwarf galaxy. In simulation LW, gas cannot cool and condense as efficiently as in simulation NOFB because molecular hydrogen, the main coolant in low-mass primordial galaxies, is photodissociated by LW radiation. In simulation LW+RT, the central gas densities are, on average, further reduced as photoionization heating drives the gas away from the halo center. The radiative feedback on the distribution of baryons implies a significant change in the gravitational potential and, in turn, in the gravitational pull on the dark matter, which hence remains less centrally concentrated.

The comparison of the dark matter density profiles in the simulations with and without feedback demonstrates that the ability of gas to cool and condense to high densities is crucial for establishing the singular isothermal density profile seen in the simulation without feedback (e.g., Wise & Abel 2008b; Zemp et al. 2012). That gas condensation can lead to a more concentrated dark matter density distribution is well known and is usually described using the framework of halo contraction models (Blumenthal et al. 1986; Gnedin et al. 2004, 2011). One may then expect the reverse of this process, i.e., the removal of gas from the halo center, to lead to a reduction in the concentration of the dark matter. Indeed, previous works have demonstrated the ability of SN explosions to lower the central dark matter densities in dwarf galaxies with masses ≳ 10⁹ M_☉ at z ≲ 10 (e.g., Navarro et al. 1996; Mashchenko et al. 2006, 2008; Governato et al. 2012).

Here we have shown that radiative feedback can have a qualitatively similar effect on the dark matter distribution in high-redshift minihalos. The reduction in the central dark matter densities due to radiative feedback is similar to that seen in the adaptive mesh refinement simulations of Wise & Abel (2008b). In our simulation LW+RT, the centrally suppressed dark matter profiles shown in the left panel of Figure 3 can be approximated by a Navarro et al. (1997) profile and concentration parameter ∼2 at z = 18.5 and 15.5. The effect of radiative feedback on the dark matter profiles of high-redshift low-mass halos was also investigated by Ricotti (2003). The shape of dark matter profiles in a cosmological simulation including gas dynamics and radiative feedback was found to be similar, on average, to the shape in a simulation that was identical except that it only treated the dynamics of the dark matter. Our finding that radiative feedback can counteract the effects of gas condensation on the dark matter density profile is consistent with the results in Ricotti (2003) and Wise & Abel (2008b). However, Ricotti (2003) also points out that there are significant statistical variations in the shape of the dark matter profile of halos at fixed mass and redshift in the simulation that included gas physics. We therefore caution against overinterpreting our conclusions drawn from the investigation of a single dwarf halo.

Figure 4 shows a sequence of snapshots of the gas density distribution centered on the dwarf galaxy in simulation LW+RT. At the final simulation redshift, the gas at the halo center is organized in a massive central, spherical clump, an inner compact disk, and an outer extended disk. The panels can be compared with Figure 8 of Pawlik et al. (2011), which shows the gas densities in our earlier simulation identical to simulation LW+RT discussed here except that it did not include star formation or radiation. The comparison shows that the radiative feedback is weak and the inclusion of star formation and LW and ionizing radiation does not prevent the formation of the disks seen in that earlier simulation.

Zoom In Zoom Out Reset image size

The sequence of events leading to the formation of the two disks is very similar with and without the inclusion of radiation, and the reader may therefore refer to Pawlik et al. (2011) for additional details. A major merger at redshift z ≲ 15 channels gas in the halo center. This leads to the formation of the first gaseous disk by z = 13.5. The disk subsequently develops spiral arms. The spiral arms exert torques that imply an outward transport of angular momentum in the disk gas. As a consequence, the disk shrinks in size. At z ≈ 12, a sequence of minor mergers replenishes the halo center with gas, producing the second gaseous disk. This disk remains spatially extended until the end of the simulation. The second disk surrounds the first disk, and the two have orientations tilted with respect to each other. This tilt is an interesting consequence of the hierarchical assembly of the emerging dwarf galaxy (e.g., Roškar et al. 2010; Pawlik et al. 2011). The formation redshifts of the disks are similar in all simulations and hence insensitive to the inclusion of radiation.

Photoheating creates low-density regions in the disk gas, leaving behind a disk morphology more complex than in the absence of photoheating (compare, e.g., the bottom right panel of Figure 4 with the corresponding panel in Figure 8 of Pawlik et al. 2011). These regions remain, however, locally confined, and the disks remain, as a whole, intact and relatively unaffected by photoheating. This insensitivity of the gas to photoheating is consistent with the fact that at the time of formation of the disks, i.e., at z ≲ 15, the galaxy halo has already entered the atomic cooling regime. Hence, its gravitational potential is sufficiently deep to confine photoheated gas in its interior. A positive feedback loop which consists of high gas densities confining the photoheated H ii regions allows the gas to collapse to increasingly higher densities, further confining the photoheated gas.

The left panel of Figure 5 shows the Toomre parameter, averaged in annular bins, of the disk gas at the final simulation redshift, Q = c_sκ/(πGΣ), where c_s is the adiabatic sound speed, κ = (4Ω² + rdΩ²/dr)^1/2 is the epicyclic frequency, and Ω is the angular velocity of the disk gas (Toomre 1964). In all three simulations Q ≳ 1 at radii larger than the gravitational softening length, implying that the disk configuration is stable against fragmentation. Note that the inclusion of LW and ionizing radiation increases disk stability by increasing the sound speed. The inclusion of ionizing radiation further helps to preserve the disks because photoheating evaporates the gas from the low-mass halos merging with the dwarf galaxy. In the absence of feedback, on the other hand, baryon-rich low-mass halos that pass through the disks can disturb the disks significantly (e.g., Figure 3 in Pawlik et al. 2011).

Zoom In Zoom Out Reset image size

The middle and right panels of Figure 5 show spherically averaged profiles of the SFR surface densities (middle), and the Kennicutt–Schmidt relation (right), i.e., the relation between SFR surface density and gas surface density, at the final simulation redshift. The Kennicutt–Schmidt relation is characterized by a power-law behavior and a suppression of star formation below gas surface densities ≲ 10² M_☉ pc⁻². The suppression is a result of the star formation recipe, which limits star formation to densities above 500 cm⁻³. At the final simulation redshift, such densities are realized both in the central region and in the disks, and the galaxy shows a spatially extended morphology of star formation. In contrast, star formation in the galaxy before disk formation is limited to the central region (see Figure 4). The Kennicutt–Schmidt power-law behavior can be understood by writing $\Sigma _{\rm SFR} = \tau _{\rm gas}^{-1}\Sigma _{\rm gas}$, where $\tau _{\rm gas} = \rho / \dot{\rho }_{\star }$ is the gas consumption time. Using τ_gas = τ_⋆ (Equation (1)), we find $\Sigma _{\rm SFR}= \tau _\star ^{-1} \Sigma _{\rm gas}$. This relation is shown with the dotted line in Figure 5, and the results from the simulation are in close agreement with it.

Taken at face value, the large Toomre Q values that indicate stability against disk fragmentation seem inconsistent with the non-zero SFRs of the disks. However, our simulations do not have sufficient resolution or the physical detail required to study fragmentation of the disks into individual stars. We have therefore employed a phenomenological model for star formation according to which stars form from gas with densities larger than the adopted threshold density for star formation, an approach employed in most galaxy formation simulations. It remains open if, at higher resolution or greater physical detail, the disks in our simulation would fragment to form stars, or if they would remain stable, possibly feeding a central massive black hole (e.g., Eisenstein & Loeb 1995; Koushiappas et al. 2004; Lodato & Natarajan 2006). Addressing these issues in cosmological simulations such as the simulations here is computationally challenging. Simulations of isolated disk galaxies have shown that star formation is slower in disks that are more stable as quantified by the smallest Toomre Q value in the disk (Li et al. 2005). It may therefore be that star formation in high-redshift low-mass disk galaxies is less efficient than implied by our simulations.

Figure 7 quantifies the ionization and thermal history (left) as well as the evolution of LW intensities and of the molecular hydrogen fraction (right) around the dwarf galaxy in simulation LW+RT (blue curves). The neighborhood considered here is defined, for simplicity of geometry and to reduce boundary artifacts, as the sphere of comoving radius 5r_{vir, com}(z = 10), where r_{vir, com}(z = 10) ≈ 34 kpc comoving is the virial radius of the dwarf galaxy at z = 10. We center the sphere on the comoving position of the most bound particle of this galaxy at z = 10, and this is the same position on which the slices in Figure 6 are centered. The spherical neighborhood is contained in the refinement region at all redshifts. Our qualitative discussion below is not sensitive to the precise definition of the neighborhood adopted here. Volume-weighted (mass-weighted) averages are obtained by weighting with the SPH kernel volume (particle mass).

Zoom In Zoom Out Reset image size

The fraction of the volume ionized (solid) increases with decreasing redshift until z ≲ 15, after which it remains approximately constant at $\langle \eta _{{\rm H\,\scriptsize{II}}}\rangle _{V} \approx 0.5$. The fraction of the mass ionized (dashed) shows a qualitatively similar behavior, but it slightly decreases after z ≲ 15 and reaches $\langle \eta _{{\rm H\,\scriptsize{II}}} \rangle _{m} \approx 0.2$ at the final simulation redshift. The redshift below which the ionized fractions no longer increase is similar to the redshift at which the galaxy evolves into an atomically cooling object. The fraction of the mass ionized is initially close to but slightly larger than the fraction of the volume ionized, but this changes at z ≲ 19 after which the fraction of ionized mass becomes increasingly lower than the fraction of the ionized volume. Reionization thus proceeds from the inside-out, from the dense halo gas into the diffuse IGM, as expected (e.g., Wise & Abel 2008a). The mass- and volume-weighted average gas temperature of the ionized gas, which we have defined here as gas with ionized hydrogen fraction $\eta _{{\rm H\,\scriptsize{II}}} > 10^{-3}$, fluctuates between 3000 and 10, 000 K, in good agreement with the results in Figures 12 and 13 in Wise & Abel (2008a).

Figure 7 shows that the ionization of the region around the main galaxy in simulation LW+RT is accompanied by an increase in the average LW intensities. Before the formation of the first star, the average LW intensity is close to but slightly smaller than the intensity of the imposed LW background (dotted line), a consequence of self-shielding. Emission of stellar radiation increases the LW intensity above the background, and the latter quickly becomes unimportant. At z ≲ 16 the LW intensities increase more rapidly and approach the intensities found in simulation LW that only included LW radiation (red curves), and the LW radiation efficiently destroys the molecular hydrogen. However, the mass-weighted fraction remains significantly larger than the volume-weighted fraction, mostly because of self-shielding. The comparison with simulation LW shows that the inclusion of ionizing radiation promotes the formation of molecular hydrogen (e.g., Haiman et al. 1996a; Ricotti et al. 2001; Oh & Haiman 2002).

The dwarf galaxy highly but not fully ionizes its neighborhood by z = 11. This is likely due to a combination of reasons. First, as the galaxy evolves into an atomically cooling object, photoheating is no longer able to substantially reduce the gas densities inside it (see Section 2). Ionizing photons emitted by the stellar sources are then efficiently consumed by recombinations in the dense gas, thus limiting the ionizing impact on the surrounding IGM. Second, as we will discuss in Section 5.2, star formation in the dwarf galaxy and its neighboring galaxies exerts a strong negative feedback on neighboring low-mass galaxies, suppressing star formation inside them. These galaxies thus cannot significantly contribute to reionizing the gas. Finally, the considered volume lacks more massive galaxies, which would be robust against the feedback from the dwarf galaxy and could potentially help reionizing the gas. This is an artifact of the refinement region being too small to contain these rarer halos. Note also that the galaxy resides in a highly overdense region, and our computation of the mean ionized fraction includes contributions from the self-shielded neutral gas in halos inside this region.

Figure 8 shows the evolution of the minimum virial mass of star-forming halos in the refinement region, as found in simulation LW+RT in which stars emitted both LW and ionizing radiation (LW+RT; blue). The corresponding evolutions in the simulation in which stars emitted LW but not ionizing radiation (LW; red), and in the simulation in which no radiation was present (NOFB; black), are also shown. We computed the minimum mass both by considering halos that form stars for the first time (crosses), which is hereafter referred to as the minimum collapse mass, and by considering all halos, including those that have previously formed stars (boxes with matching but lighter colors). Our definition of the minimum collapse mass is insensitive to our choice of the threshold density for star formation n_{H, SF} = 500 cm⁻³, at least in the absence of feedback from star formation, because the transition from densities n_H ≳ 10 cm⁻³ to n_H ≲ 10³ cm⁻³ which bracket the adopted star formation threshold density then occurs within a narrow range of halo masses (see the bottom left panel in Figure 9).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the absence of LW or ionizing radiation, the minimum collapse mass is consistent with the mass of halos at virial temperature T_vir ≈ 2200 K (Equation (6) with μ = 1.22) independent of redshift. This mass is consistent with but slightly larger than the masses ∼5 × 10⁵–10⁶ M_☉ at the time of gas collapse in previous simulations of minihalos (e.g., Yoshida et al. 2003; O'shea & Norman 2007; Wise & Abel 2008a; Greif et al. 2008). The small difference is probably a numerical artifact of our limited resolution, which is lower than those realized in the previous works. Differences are also expected because the mass of the halo that forms the first star is subject to statistical uncertainties related to the low number of investigated halos (e.g., Greif et al. 2008; O'shea & Norman 2008), and because of differences in the strength of dynamical heating from tidal interactions and mass accretion (Yoshida et al. 2003; O'shea & Norman 2007). The inclusion of LW radiation increases the minimum collapse mass. This increase in the minimum collapse mass due to photodissociation by LW radiation has been investigated in detail in a number of previous works (e.g., Machacek et al. 2001; Mesinger et al. 2006; O'shea & Norman 2008; Safranek-Shrader et al. 2012).

The inclusion of both LW and ionizing radiation increases the minimum collapse mass in a similar manner but often by a smaller factor, demonstrating a positive feedback from the enhanced formation of molecular hydrogen in ionization fronts and fossil H ii regions. Eventually, however, the negative feedback from photoheating outweighs the positive feedback, and there is no new star-forming halo inside the refinement region below z ≲ 13. Star formation can still proceed in halos that have previously formed stars, but the masses of these halos are significantly larger than those in simulation NOFB. This is the result of the Jeans filtering of the IGM, which impedes the accretion of gas on low-mass halos. The increased scatter in the minimum collapse mass in the presence of photoionization is an expression of the local nature of feedback from photoheating, which is limited to inside the H ii regions.

The top panels of Figure 9 show the specific SFRs, i.e., the SFRs divided by virial mass, both as a function of virial mass (left) and of stellar mass (middle), and the stellar mass fractions (right) for the halos inside the high-resolution region at the final simulation redshift, z = 11. The bottom panels of Figure 9 show the maximum gas densities inside the virial radius (left) and the baryon mass fractions as function of virial (middle) and stellar mass (right) for these halos.

In the simulation without radiation (NOFB), the specific SFR ∼ 2 × 10⁻¹⁰ yr⁻¹ is nearly independent of halo mass in the range ∼10⁷–10⁹ M_☉. The specific SFR is reduced in a fraction of the halos with masses ≲ 10⁷ M_☉ just above the minimum collapse mass, which we attribute primarily to the effects of dynamical heating during the gravitational collapse of these halos (e.g., Yoshida et al. 2003). Indeed, the fact that the transition to high central gas densities occurs within a finite range of halo masses shows that dynamical heating plays a non-negligible role, and in the lowest mass halos this may be amplified by the limited mass resolution we afford. Note that some of the halos show an increased baryon fraction f_bar ≳ 0.2 as a result of dynamical interaction and ongoing mergers with other halos.

The inclusion of LW radiation implies a complete suppression of star formation only in halos with masses ≲ 2 × 10⁷ M_☉, corresponding to virial temperatures significantly below 10⁴ K (vertical dashed line). This is in qualitative agreement with the results from previous high-resolution simulations of the collapse of minihalos in the presence of an LW radiation background. These works demonstrated that as the minihalo mass approaches the atomic cooling limit, the presence of LW radiation cannot prevent the buildup of molecular hydrogen, because it is catalyzed by the elevated electron fraction inside central structure formation shocks (e.g., Wise & Abel 2007a; O'shea & Norman 2008), and also because of self-shielding (e.g., Ahn & Shapiro 2007; Susa 2007). In contrast, the additional inclusion of ionizing radiation can potentially suppress star formation in all halos below the atomic cooling limit by evaporating the gas from the halo centers. However, as we have already discussed above in Figure 8, inspection of the stellar mass fractions reveals that in simulation LW+RT, star formation is possible down to lower halo masses than in simulation LW, a consequence of the positive feedback from ionization ahead of ionization fronts and recombinations of H ii regions on the formation of molecular hydrogen (Ricotti et al. 2002, 2008).

The simultaneous inclusion of both LW and ionizing radiation reduces the average baryon fractions in nearly the full range of simulated halo masses. The reduction of the baryon fraction is strongest, on average, for the lowest-mass halos with masses ≲ 10⁶ M_☉, and for halos in the intermediate-mass regime, with masses in the range 10⁷–10⁸ M_☉. Because halo masses ≲ 10⁶ M_☉ are below the minimum collapse mass, halos in this mass range do not form stars, and hence their baryon fraction is reduced with respect to that in simulation NOFB due to the effects of radiation from external sources, and due to Jeans filtering in the photoheated IGM. More massive halos may accrete gas despite external feedback and Jeans filtering, which explains why the baryon fraction increases, on average, in the halo mass range ∼10⁶–10⁷ M_☉. Halos in the range 10⁷–10⁸ M_☉ are efficiently forming stars, and therefore their baryon fractions are strongly reduced by radiative feedback from internal sources. Finally, halos with masses ≳ 10⁸ M_☉, corresponding to virial temperatures ≳ 10⁴ K, are robust not only against Jeans filtering and feedback from external sources, but also against the feedback from internal sources. The baryon fraction of these halos is reduced primarily because there was not yet enough time for accretion to compensate for the mass loss in past episodes of gaseous outflows.

Ricotti et al. (2008) investigated radiative feedback from the first galaxies using cosmological simulations of comoving size 1–2 h⁻¹ Mpc including processes such as chemical enrichment not treated here. Utilizing resolution similar to that realized here, they demonstrated that the positive feedback due to enhanced molecular hydrogen formation near ionization fronts and inside recombining H ii regions can be very strong and compensate for the negative feedback due to photodissociation of molecular hydrogen (see also Ricotti et al. 2002). Our simulations may underestimate this positive feedback because we do not capture the spectral hardening of the ionizing radiation, implying less ionization to stimulate the formation of hydrogen ahead of ionization fronts. Ricotti et al. (2008) further demonstrated that the stellar mass fractions of the lowest mass halos are characterized by a large scatter at fixed halo mass. Our simulations do not find as large a scatter, which may be a consequence primarily of the larger star particle mass and statistical limitations caused by the small size of the refinement region in our simulations, although other differences between the simulations may contribute. We do however find a similar large scatter in the baryon fraction at low halo masses, and we also agree with Ricotti et al. (2008) that this scatter is driven by radiative feedback from both external and internal sources.

A.1.1. traphic

A.1.2. Absorption of Ionizing Photons by Hydrogen/Helium

A fraction 1 − exp [ − τ(ν)] of each photon packet in frequency bin ν emitted or transmitted from particle i at position r_i to SPH neighbor j at position r_j, separated by the propagation distance d_ij, is absorbed. The optical depth τ(ν) is the sum τ(ν) = ∑τ_α(ν) of the optical depths τ_α(ν) of each absorbing species $\alpha \in \lbrace {\rm H\,\scriptsize{I}, He\,\scriptsize{I}, He\,\scriptsize{I}}\rbrace$, and $\tau _\alpha (\nu) = \int _{\mathbf {r}_i}^{\mathbf {r}_j} dr\ \sigma _\alpha (\nu) n_\alpha (\mathbf {r}) \approx \sigma _\alpha (\nu) n_\alpha (\mathbf {r}_j) d_{ij}$. The number of photons absorbed by species α is δN_{abs, α}(ν) = [w_α(ν)/∑_αw_α(ν)]δN_abs(ν), where the weights w_α(ν) = τ_α(ν) (Pawlik & Schaye 2011) and δN_abs(ν) = ∑_αδN_{abs, α}(ν) is the total number of absorbed photons.

In this work we choose, for reasons of computational efficiency, to transport hydrogen-ionizing radiation using a single frequency bin, i.e., we set N_ν = 1. In this bin, we adopt frequency-averaged photoionization cross-sections 〈σ_γα〉 for the absorption by the individual species α using the gray approximation (e.g., Mihalas & Weibel Mihalas 1984),

$$\begin{equation} \langle \sigma _{\gamma \alpha } \rangle \equiv \int _{\nu _\alpha }^\infty d\nu \frac{4 \pi J_{\nu }(\nu)}{h_{\rm P} \nu } \sigma _{\gamma \alpha } (\nu) \times \left[\int _{\nu _{{\rm H\,\scriptsize{I}}}}^\infty d\nu \frac{4 \pi J_{\nu }(\nu)}{h_{\rm P}\nu } \right]^{-1}, \end{equation} \tag{ A1 }$$

where J_ν(ν) is the mean intensity of the radiation field and h_Pν_α is the ionization potential. We assume that the mean intensity J_ν(ν) is characterized by a blackbody spectrum of temperature T_BB = 10⁵ K, appropriate for the emission of radiation by the first stars (e.g., Schaerer 2003), which gives $\langle \sigma _{\gamma {\rm H\,\scriptsize{I}}} \rangle = 1.63 \times 10^{-18} \,\mbox{cm}^{2}$, $\langle \sigma _{\gamma {\rm He\,\scriptsize{I}}} \rangle = 2.28 \times 10^{-18} \,\mbox{cm}^{2}$, and $\langle \sigma _{\gamma {\rm He\,\scriptsize{II}}} \rangle = 6.19 \times 10^{-20} \,\mbox{cm}^{2}$. We have employed the fits to the photoionization cross-sections from Verner et al. (1996).

At fixed mean intensity J_ν(ν), the gray approximation implies photoionization rates identical to those inferred from a full multi-frequency treatment. However, because of the use of only a single frequency bin, our simulations ignore the hardening of the radiation spectrum with distance from the source caused by the preferential absorption of lower energy photons with larger absorption cross sections (see, e.g., the discussion in Pawlik & Schaye 2011).

A.1.3. Photoionization and Photoheating Rates

At the end of each RT time step the photoionization and photoheating rates are computed. We determine the photoionization rates directly from the total number of photons N_{abs, α}(ν) = ∑δN_{abs, α}(ν) absorbed by a given SPH particle during the RT time step Δt_r in frequency bin ν,

$$\begin{equation} \Gamma _{\gamma \alpha } = \frac{\sum _\nu N_{\rm abs, \alpha }(\nu)}{\eta _\alpha N_{\rm H} \Delta t_{\rm r}}, \end{equation} \tag{ A2 }$$

where N_H ≡ m_gasX_H/m_H is the number of hydrogen atoms associated with the SPH particle of mass m_gas. This ensures that the same number of photons that have been removed from the simulation is used to determine the ionization balance and temperature of the ionized gas, i.e., photon conservation (Abel et al. 1999).

The heating rate per atom due to photoionization of species α, assuming a single frequency bin in the gray approximation, is $\mathcal {E}_{\gamma \alpha } = \Gamma _{\gamma \alpha } \langle \epsilon _\alpha \rangle$, where

$$\begin{eqnarray} \langle \epsilon _\alpha \rangle &=& \int _{\nu _\alpha }^\infty d\nu \frac{4 \pi J_{\nu }(\nu)}{h_{\rm P} \nu } \sigma _{\gamma \alpha } (\nu) (h_{\rm P}\nu - h_{\rm P}\nu _\alpha) \nonumber\\ && \times \left[ \int _{\nu _\alpha }^\infty d\nu \frac{4 \pi J_{\nu }(\nu)}{h_{\rm P} \nu } \sigma _{\gamma \alpha } (\nu) \right]^{-1} \end{eqnarray} \tag{ A3 }$$

is the average energy of the absorbed ionizing photons in excess of the photoionization threshold. Because we assume that the mean intensity J_ν(ν) is characterized by a blackbody spectrum of temperature T_BB = 10⁵ K, the average excess energies for photoionization of hydrogen and helium are $\langle \epsilon _{{\rm H\,\scriptsize{I}}}\rangle = 6.32 \,\mbox{eV}$, $\langle \epsilon _{{\rm He\,\scriptsize{I}}}\rangle = 8.70 \,\mbox{eV}$, and $\langle \epsilon _{{\rm He\,\scriptsize{II}}}\rangle = 7.88 \,\mbox{eV}$.

A.1.4. Radiation–Hydrodynamical Coupling

We have previously described a number of tests of our implementation of traphic in gadget for the RT on static density fields. In Pawlik & Schaye (2008), we carried out monochromatic RT simulations of increasing complexity, assuming that the gas is composed only of hydrogen and has a fixed temperature. By comparing with reference solutions such as those published in Iliev et al. (2006), we demonstrated the ability of traphic to accurately capture the evolution of ionization fronts, and to reproduce the ionized fractions. We also showed that traphic is able to produce sharp shadows behind opaque absorbers, which gives rise to the typical "butterfly" shape of the ionized regions (e.g., Abel et al. 1999). In Pawlik & Schaye (2011), we extended our implementation of traphic to enable the transport of multi-frequency radiation, to account also for the ionization of helium, and to compute the evolution of the gas temperature due to radiative cooling and photoionization heating. Results from test simulations with this new implementation were in excellent agreement with reference solutions, such as those in Iliev et al. (2006).

The implementation of traphic in gadget that we use in this work differs from that described in Pawlik & Schaye (2008, 2011) in two main respects. First, we have substituted the explicit solver for the evolution of the chemistry and temperature of primordial atomic gas in the presence of photoionizations used and tested in Pawlik & Schaye (2011) with the implicit solver used in Pawlik et al. (2011) and described in Section 2.2. This solver accounts for additional species, such as molecular hydrogen, and employs different rates for the atomic and molecular physics. It has been described and extensively tested in a number of publications (e.g., Johnson & Bromm 2006; Greif et al. 2010). However, it has not yet been tested in combination with traphic. In Appendix A.2.1, we therefore repeat the RT test simulation on a static cosmological density field from Pawlik & Schaye (2011). The results of this test are in excellent agreement with our previous results, demonstrating the validity of our implementation.

Second, the current work is the first to employ traphic in radiation–hydrodynamic simulations that account for the feedback of photoionization heating on the gas dynamics. We have performed the radiation–hydrodynamical tests from Iliev et al. (2009), and we have achieved excellent agreement with the published reference results in all of them. In Appendix A.2.2 we provide a discussion of one of the tests relevant to the present work, the simulation of the evaporation of a minihalo by an internal ionizing source.

Finally, in Appendix A.2.3, we investigate how our adoption of a limited photon propagation speed and a finite angular sampling affects the early evolution of the minihalo in simulation LW+RT analyzed in the main text.

A.2.1. Radiative Transfer and Chemical and Thermal Evolution

In this appendix, we repeat Test 4 of the cosmological RT code comparison project (Iliev et al. 2006) that we have previously discussed in Pawlik & Schaye (2011), using a chemistry solver different from the one employed here. The test involves the simulation of H ii regions around multiple ionizing sources in a static cosmological density field. It was designed to capture important aspects of state-of-the-art simulations of hydrogen reionization, such as the delayed propagation in or trapping of ionization fronts by dense gas.

The setup of this test is identical to that of Test 4 in Pawlik & Schaye (2011), to which we refer the reader for a detailed description. Briefly, the initial conditions are provided by a snapshot (at redshift z ≈ 8.85) from a cosmological N-body and gas-dynamical uniform-mesh simulation. The simulation box is L_box = 0.5 h⁻¹ Mpc comoving on a side, where h = 0.7, and is uniformly divided into N_cell = 128³ cells. We Monte Carlo sample this density field to replace the mesh cells with N_SPH = N_cell = 128³ SPH particles. The gas consists purely of atomic hydrogen and is assumed to be initially neutral at temperature T = 100 K. The ionizing sources are chosen to correspond to the 16 most massive halos in the box, and are assumed to have blackbody spectra with temperature T_BB = 10⁵ K.

For comparison with Pawlik & Schaye (2011) we solve the time-independent RT equation with an angular resolution of N_c = 32, set the number of neighbors to which sources emit radiation to $\tilde{N}_{\rm ngb}=32$, employ a time step Δt_r = 10⁻⁴ Myr, and transport photons only over a single inter-particle distance per time step. We transport radiation using a single frequency bin, employing the gray photoionization cross-section $\langle \sigma _{\gamma {\rm H\,\scriptsize{I}}}\rangle = 1.63 \times 10^{-18} \,\mbox{cm}\,\mbox{s}^{-1}$. We assume that each photoionization adds $\langle \epsilon _{{\rm H\,\scriptsize{I}}} \rangle = 6.32 \,\mbox{eV}$ to the thermal energy of the gas, as appropriate for the adopted blackbody spectrum.

In Figure 11, we show results of this simulation and compare them with those presented in Pawlik & Schaye (2011) and also with the reference results reported in Iliev et al. (2006). The agreement is very good if one accounts for the differences in the atomic physics and the approximations employed (see Pawlik & Schaye 2011 for a detailed discussion).

Zoom In Zoom Out Reset image size

A.2.2. Radiation–Hydrodynamical Coupling

In this appendix, we perform Test 6 of the cosmological RT code comparison project (Iliev et al. 2009). This test consists of the simulation of the expansion of an H ii region driven by an ionizing point source at the center of a spherically symmetric region with a steeply declining gas density profile. In contrast to the test described in the previous appendix, the gas distribution is not forced to remain static but may evolve in response to photoionization heating from the central source. The test situation shares some of the main features of the evolution of gas densities inside high-redshift cosmological minihalos under photoionization from a central massive metal-free star.

The physical setup of the test is as follows. The initial gas densities are described by a cored isothermal density profile,

$$\begin{equation} n_{\rm H}(r) = \left\lbrace \begin{array}{@{}ll}n_0 & \quad {{\rm if}\ r \le r_0,}\\ n_0 (r/r_0)^{-2} & \quad {\rm otherwise.}\\ \end{array} \right. \end{equation} \tag{ A4 }$$

We follow Iliev et al. (2009) and set n₀ = 3.2 cm⁻³ and r₀ = 91.5 pc, and we assume that the gas, consisting solely of atomic hydrogen, is initially neutral at temperature 100 K. The central source has a constant ionizing luminosity of 10⁵⁰ photons s⁻¹ throughout the simulation, and is characterized by a blackbody spectrum with temperature 10⁵ K. Following Iliev et al. (2009), all relevant cooling processes are included, except for Compton cooling. In this test, gravitational forces are ignored.

We implement the power-law density profile by applying a local radial stretch along the direction toward the central source to the inter-particle distances of SPH particles that are initially distributed uniformly with density $\tilde{n}_{\rm H}$. Denoting the initial particle positions with $(\tilde{r}, \theta , \phi)$, such a stretching is expressed by the transformation $(\tilde{r}, \theta , \phi) \rightarrow (r, \theta , \phi)$, where θ and ϕ are the usual two angles needed to specify a position in spherical coordinates. Mass conservation then requires the new coordinates (r, θ, ϕ) to satisfy

$$\begin{equation} n_{\rm H}(r) r^2 \sin \theta dr d\theta d\phi =\tilde{n}_{\rm H} \tilde{r}^2 \sin \theta d\tilde{r} d\theta d\phi. \end{equation} \tag{ A5 }$$

Substituting a power-law density profile, n_H(r)∝rⁿ, the last equation can be integrated to yield $r \propto \tilde{r}^{3/(3-n)}$. In the present case, n = 2, and hence $r \propto \tilde{r}^3$. Therefore, a singular isothermal profile can be obtained by stretching the initial inter-particle distances between uniformly distributed particles along the radial directions toward the center by a factor $\tilde{r}^2$. The initial uniform density field is generated by placing particles randomly in the simulation box. To reduce the associated shot noise, the resulting particle distribution is regularized by evolving the particles under the influence of a reversed-sign (i.e., repulsive) gravitational force until the particles settle down into a glass-like quasi-equilibrium (White 1996). After generating the singular isothermal density profile, we replace the central region within r₀ with a uniform glass-like particle distribution to introduce the central core.

The numerical parameters used in this test are as follows. We place the ionizing source in the center of a box with linear size 1.6 kpc. We sample the gas inside the core with radius r₀ using 30,000 particles, which corresponds to an equivalent uniform grid resolution of 160³ cells. This resolution is slightly higher than the uniform grid resolution employed in Iliev et al. (2006), which was 128³ cells. We use an angular resolution N_c = 8, a number of neighbors for the emitting source $\tilde{N}_{\rm ngb}=32$, employ an RT time step Δt_r = 10⁻⁵ Myr, and transport photons only over a single inter-particle distance per time step. We limit the sizes of the individual particle hydrodynamical time steps to be less than 0.1 Myr. The radiation is transported using a single frequency bin and assuming a gray photoionization cross-section $\langle \sigma _{\gamma {\rm H\,\scriptsize{I}}}\rangle = 1.63 \times 10^{-18} \,\mbox{cm}\,\mbox{s}^{-1}$, as well as that each photoionization adds $\langle \epsilon _{{\rm H\,\scriptsize{I}}} \rangle = 6.32 \,\mbox{eV}$ to the thermal energy of the gas. The parameter values are chosen to make contact with the conditions in the simulations presented in the main part of this paper, and to ensure numerical convergence. The test results are not critically sensitive to the adopted parameter values. However, we note that because of our choice of emitting photon packets only once per RT time step, the angular sampling depends on the RT time step, and a larger RT time step would imply an increased scatter in the ionized fraction, temperature, and density profiles (see Appendix A.2.3).

The results of the test are shown and discussed in Figure 12. The results obtained with traphic are in very good agreement with the reference results published in Iliev et al. (2009).

Zoom In Zoom Out Reset image size

A.2.3. Effect of Limited Photon Packet Propagation Speed and Finite Angular Sampling

In this appendix, we discuss the effects of our approximations regarding the limited photon propagation speed and the emission of photon packets by ionizing sources along a finite set of random directions per RT time step. To this end, we have repeated the RT around the first stellar burst in the minihalo progenitor of simulation LW+RT. The first test simulation that we have performed, hereafter test-base, employed RT parameters identical to those employed in LW+RT. In particular, in this simulation ionizing sources emitted photon packets along N_EC = 8 random directions per RT time step, and the photon propagation was done at the speed of light but limited to a single inter-particle distance per RT time step. To investigate the effect of the limitation in the propagation speed, simulation test-base is compared with a simulation in which photon packets propagate at the speed of light with no additional limitation regarding the number of inter-particle distances per RT time step (hereafter test-speed; and with ionizing sources emitting photons into N_EC = 8 random directions per RT time step as in test-base). To investigate the effect of the angular sampling, simulation test-speed is compared with a simulation in which sources emit photon packets along 10³ times more random directions per RT time step (hereafter test-sampling; and transporting photons at the speed of light as in test-speed). The latter is implemented by repeating the emission of radiation along the N_EC = 8 random directions 10³ times per RT time step, with each emission using a different set of N_EC = 8 random directions and emitting photon packets that contain by a factor 10³ fewer photons.

Figure 13 shows the gas temperature in a thin slice through the center of the minihalo progenitor at z = 22, at the end of the stellar burst, in the three test simulations test-base (left), test-speed (middle), and test-sampling (right). The final extent of the photoheated region is similar in simulations test-base and test-speed despite the difference in propagation speeds. This is because the RT with traphic conserves the number of ionizing photons that are transmitted and absorbed, and because the final extent of photoheated H ii regions is set primarily by this number. However, limitations of the speed at which photons are propagated can lead to differences in extent and shape of the developing H ii regions early during their evolution (e.g., Pawlik & Schaye 2008). The increased angular sampling in simulation test-sampling reduces the noise that characterizes the shape of the photoheated region in simulation test-base, as expected. However, this noise is already small in simulation test-base, which hence yields a photoheated region similar in shape to that in simulation test-sampling. We note (but do not show) that the evolution of the properties of the minihalo, such as, e.g., the maximum gas density, is nearly identical in all three test simulations.

Zoom In Zoom Out Reset image size