CONFINED POPULATION III ENRICHMENT AND THE PROSPECTS FOR PROMPT SECOND-GENERATION STAR FORMATION

Jeremy S. Ritter; Chalence Safranek-Shrader; Orly Gnat; Miloš Milosavljević; Volker Bromm

doi:10.1088/0004-637X/761/1/56

1. INTRODUCTION

The first stars forming from metal-free gas in the early universe, the Population III (Pop III), profoundly transformed their cosmic environment. Radiation from these stars first ionized the otherwise neutral, chemically pristine gas. In their final demise, likely involving supernovae or the collapse into a black hole, the Pop III stars injected momentum, energy, and nucleosynthetic products into the environment (Karlsson et al. 2012). The amplitude and character of these effects depended sensitively on the Pop III stellar masses, as well as on theoretically uncertain aspects of stellar evolution. Early numerical investigations of the formation of Pop III stars in metal-free cosmic halos concluded that the stars were isolated and rather massive, with masses ≳ 100 M_☉ (Bromm et al. 1999, 2002; Abel et al. 2000, 2002; Gao et al. 2007; O'Shea & Norman 2007). Recent numerical investigations, however, have shown that angular-momentum-aided fragmentation (Stacy et al. 2010; Clark et al. 2011; Greif et al. 2011, 2012) and radiative termination of accretion (Hosokawa et al. 2011; Stacy et al. 2012a; see also McKee & Tan 2008 for an analytical approach) during protostellar accretion could limit the masses to as low as a few tens of solar masses. The final stellar mass range remains uncertain but could plausibly be ∼1–50 M_☉ with a characteristic mass ∼10 M_☉.

For a wide range of stellar masses and ambient gas densities, Pop III stars have sufficient luminosities and effective temperatures to form extended H ii regions, unless they form in especially high mass halos (e.g., Whalen et al. 2008) or with unusually high accretion rates during the protostellar phase (e.g., Hosokawa et al. 2012). The raised pressure of the interior ionized gas produces a hydrodynamic response that reduces the circumstellar gas density, and it is the density profile modified by H ii region dynamics that the supernova blast wave expands into. The density structure imprinted by the H ii region can interact with the blast wave to drive turbulence and compositional mixing of the supernova ejecta.

Because the stellar ionizing photon production rate is a strong function of stellar mass (e.g., Bromm et al. 2001; Schaerer 2002), the properties of H ii regions of the first stars will depend on the mass. A metal-free star with mass >100 M_☉ could emit >10⁵⁰ ionizing photons per second, which would ionize not only the parent cosmic minihalo but a larger volume of the primordial intergalactic medium (IGM). The pressure of the photoionized gas drives a supersonic, radially propagating pressure wave, which ultimately shocks (Whalen et al. 2004, 2008; Kitayama et al. 2004; Alvarez et al. 2006; Abel et al. 2007; Yoshida et al. 2007). This expels a fraction of the baryons from the halo even before the supernova explodes. Only the densest clumps residing in the filaments streaming into the halo from the cosmic web are able to avoid ionization (Abel et al. 2007). Furthermore, the H ii regions may also be prone to violent dynamical instabilities that could create a clumpy medium (Whalen & Norman 2008), with interesting consequences for supernova blast-wave dynamics. Most of these features should carry over to lower-mass Pop III stars with lower ionizing photon production rates, but the magnitude of the effects is weaker: the temperature of the photoionized gas (at a fixed distance from the star) is lower, the pressure wave is less pronounced, and the final ionized mass fraction inside the halo and the ionizing photon escape fraction from the halo are lower.

Similarly, the outcome of the stellar demise is expected to depend sensitively on stellar mass. Restricting to nonrotating, metal-free stars with minimum mass loss (e.g., Heger & Woosley 2002; Heger et al. 2003), it is expected that the stars with masses 15–50 M_☉ undergo core collapse, leaving behind a neutron star, or with fallback, a black hole. The stars with masses 50–140 M_☉ can produce a black hole directly, though the stars above 100 M_☉ may be affected by pair-pulsational instability. The stars with masses 140–260 M_☉ explode as pair-instability supernovae (PISNe), and even more massive stars can also directly collapse into a black hole. Stellar rotation introduces new effects and significant theoretical uncertainties, e.g., the stars might first form a central black hole, and the subsequent accretion of the rotating stellar envelope onto the black hole might power an outflow that could produce an explosion (see, e.g., MacFadyen & Woosley 1999; MacFadyen et al. 2001; Kohri et al. 2005; Milosavljević et al. 2012; Lindner et al. 2012), with a peculiar nucleosynthetic imprint that can be sought in especially metal-poor stars (e.g., Iwamoto et al. 2005).

Much of the theoretical work on the formation of the first stars and galaxies (e.g., Wise & Abel 2008a; Wise et al. 2012; Greif et al. 2010) and on strategies for detecting metal-free stellar populations in the observations of galaxies and supernovae (e.g., Schaerer 2003; Johnson et al. 2009; Raiter et al. 2010; Zackrisson et al. 2011) has postulated the massive character of metal-free stars. The prospect of PISNe has been particularly interesting in view of the large expected explosion energy ∼10⁵³ erg and large nucleosynthetic yield, and many investigations have focused on this type of explosion (see, e.g., Pan et al. 2012; Hummel et al. 2012; Whalen et al. 2012). Such an explosion further reduces the residual baryonic content of the halo (Bromm et al. 2003), and the metals are dispersed far outside, thus enriching a relatively large region of the surrounding medium (e.g., Wise & Abel 2008b; Greif et al. 2010; Maio et al. 2010; Wise et al. 2012). While this enrichment has only a minor effect on the ultimate metal budget of the IGM (Scannapieco et al. 2002; Yoshida et al. 2004), it has the potential to raise the metallicity of larger collapsing halos above a "floor" (e.g., Bromm & Loeb 2003; Smith et al. 2008), itself a function of the dust content (e.g., Schneider et al. 2012), thought to be necessary for fragmentation into low-mass stars.

There has been much less work done on the cosmic imprint of Pop III stars with masses in the range tentatively favored by recent investigations (Stacy et al. 2010, 2012a; Clark et al. 2011; Greif et al. 2011; Hosokawa et al. 2011). The ionizing luminosities of these moderate-mass Pop III stars could be up to an order of magnitude lower than those of their massive counterparts, the supernova energies could be ∼10⁵¹ erg, the nucleosynthetic yields could be at most a few solar masses per supernova, and the characteristic abundance patterns could be more consistent with currently available observational data (e.g., Tumlinson 2006; Vangioni et al. 2011). Thus, the evolution of minihalos forming only moderate-mass Pop III stars should differ in critical ways from that of objects forming very massive ones.

Here, we report the first results from a program to study the detailed transport of the nucleosynthetic output of the first stars, from the point of injection in supernovae, until at least a fraction of the new metals have recondensed in descendent cosmic halos, where they may contribute to the formation of second-generation, metal-enriched stars. The specific aim of the paper is to investigate the implications of the hypothesis that the first stars had moderate masses–masses that are insufficient to produce extraordinarily extended H ii regions and ultra-energetic supernovae. We present a cosmological simulation tracking the supernova ejecta of one such intermediate-mass Pop III star. It is expected on theoretical grounds that the character of the new star formation in the descendent halo will depend on the rate of gas inflow $\dot{M}_{\rm inf}$ into the nucleus of the halo, where the gas can rapidly cool and fragment into stars. Also, the cooling rate of the gas and the compressibility of the turbulence in the flow depend on the average metallicity Z_inf and the degree of metallicity homogeneity in the gas. The simulation allows us to measure these parameters and thus shed light on the conditions for the formation of metal-enriched stars in the aftermath of a Pop III supernova.

Our focus is the gross dynamics of the supernova remnant and to trace the dispersal of the metals produced. We here do not attempt to make detailed predictions for chemical abundance patterns in second-generation stars that may in part survive to the present day and that can be compared with the extensive data on metal-poor stars in the Milky Way (e.g., Cayrel et al. 2004; Beers & Christlieb 2005; Frebel et al. 2005; Lai et al. 2008). It has been argued that the observed abundances in these ancient stars can be explained with Pop III supernovae from moderate-mass (∼15–40 M_☉) progenitors (e.g., Joggerst et al. 2010; Tumlinson 2010). This interpretation resonates with the basic assumption on the Pop III mass scale made here. In view of possibly strong selection effects in the existent metal-poor abundance data, it is important, however, to keep an open mind in this regard. One key selection bias might result from the high level of enrichment in the immediate neighborhood of a PISN (Wise & Abel 2008b; Greif et al. 2010), such that second-generation stars could already exhibit overall metallicities of 10⁻³ to 10⁻² Z_☉; most surveys targeting the extremely metal-poor tail would then miss such PISN-enriched fossils (Karlsson et al. 2008). Future surveys, such as Gaia, should provide us with essentially unbiased samples and an opportunity to constrain theoretical models.

The paper is organized as follows. In Section 2, we describe our numerical algorithm. In Section 3, we present our results. In Section 4, we discuss implications for the formation of metal-enriched star clusters in the early universe, and in Section 5, we present our main conclusions. We adopt the cosmological parameters consistent with the Wilkinson Microwave Anisotropy Probe seven-year data (Komatsu et al. 2011). We explicitly indicate whether comoving or physical length units are utilized; all other quantities, such as density, are expressed in physical units.

2. SIMULATING A SUPERNOVA IN A MINIHALO

2.1. Cosmological Initial Conditions and Gravity

We initialize our simulation in a periodic cosmological box of volume 1 Mpc³ (comoving) at redshift z = 146, large enough not only to permit the formation of such a star but also to allow the subsequent growth of the star's host halo, by accretion and merging, into a larger object that is the potential site of second-generation star formation. The initial conditions are generated with the multiscale cosmological initial conditions package GRAPHIC2 (Bertschinger 2001) with two levels of nested refinement on top of the 128³ base grid. The highest refined region had an effective resolution of 512³, corresponding to a dark matter particle mass of 230 M_☉. While we will focus on the dynamics of an initially ∼10⁶ M_☉ minihalo at redshift z ∼ 20, the region of high refinement was centered on the overdensity that first collapses to form a more massive ∼10⁸ M_☉ halo at a lower redshift. The time integration was carried out with the multipurpose astrophysical adaptive mesh refinement (AMR) code FLASH (Fryxell et al. 2000), version 3.3, employing the direct multigrid Poisson gravity solver of Ricker (2008). The AMR refinement criteria, which we describe in Section 2.3 below, are set to ensure that relevant baryonic processes are sufficiently well resolved in the simulation. In collapsing regions where baryons dominate the gravitational potential, the AMR resolution length quickly becomes smaller than the dark matter particle separation. There, we employ a mass- and momentum-conserving quadratic-kernel-based softening procedure developed in Safranek-Shrader et al. (2012) which is applied to the dark matter density variable during mapping onto the computational mesh, to render the dark matter particle contribution to the gravitational potential smooth on the scale of the local AMR grid.

2.2. Chemistry and Cooling

Prior to the initialization of the H ii region, we integrate the full nonequilibrium chemical network for hydrogen, helium, deuterium, and their chemical derivatives. The chemical and radiative cooling updates, which include, among other processes, the cooling or the heating by the cosmic microwave background (CMB), are operator split from the hydrodynamic update and are subcycled within each cell. The chemical reaction rates are the same as we have summarized in Safranek-Shrader et al. (2010), except that the chemical and radiative cooling update is isochoric in the present simulation. The network is integrated with the Bulirsch–Stoer-type semi-implicit extrapolation method of Bader & Deuflhard (1983). We test the thermodynamic evolution of the gas in Safranek-Shrader et al. (2012).

After turning on the ionizing source, which we describe in Section 2.4 below, we continue to integrate the nonequilibrium chemical network in the neutral gas, but in the ionized gas, we adopt a scheme in which the gas temperature and chemical abundances relax exponentially toward a target photoionization equilibrium. The equilibrium gas temperature and chemical abundances are functions of the photoionization parameter ξ ≡ 4πF/n_H only, where F is the radiative energy flux and n_H is the number density of hydrogen nuclei, and are precomputed with the same chemical network that we integrate in the neutral gas. The relaxation timescale is chosen to be the true, nonequilibrium thermal timescale t_therm ≡ (dln T/dt)⁻¹, defined by the instantaneous local radiative cooling and photoionization heating rates.

After inserting the supernova remnant, the gas can contain an arbitrary admixture of metals in addition to the primordial elements. The radiative point source is removed, thus discontinuing the photoionization equilibrium calculation. We revert to an integration of the nonequilibrium chemical network for the primordial species and calculate the ionization state of the metals and their contribution to the cooling by interpolation from precomputed tables. At temperatures above 8000 K, the collisional ionization equilibrium state of each of the metal species, including C, N, O, Ne, Mg, Si, S, and Fe, is calculated separately using methods described in Gnat & Sternberg (2007) and Gnat & Ferland (2012). Below 8000 K but above 200 K, we ignore the metallic contribution to the free electron abundance and tabulate a single density-independent cooling function (such that the volumetric radiative loss rate scales with the square of the density), appropriate for a gas at relatively low densities ≲ 10³ cm⁻³, using the code CLOUDY (Ferland et al. 1998). The tabulated data extend to a minimum metallicity of 10⁻³ Z_☉; below this metallicity, we set the metal cooling rate in this range of temperatures to zero. Below 200 K, we disable cooling by metal lines altogether; however, because low-temperature thermodynamics is not the focus of this work, this has no effect on the forthcoming results. The treatment of primordial species takes into account the electrons provided by the metals and vice versa, but we do not include collisional charge exchange reactions between the primordial species and the metal species. We assume that the metals have solar abundance ratios. This is an arbitrary, unnatural choice given that our metals are being produced by a core-collapse supernova, but this choice should have only a very minor effect on the thermodynamic evolution of the supernova remnant.

2.3. Adaptive Mesh Refinement Strategies

The simulation is set to adjust the FLASH refinement level ℓ based on the local baryon density ρ_b to ensure that denser regions are more highly refined. The grid separation (cell size) is related to the refinement level via Δx = 2^{−ℓ + 1 − 3}L, where L is the size of the computational box.³ As gravitational collapse drives the local fluid density ρ_b to increasingly higher values, the local refinement level is adjusted to satisfy $\rho _{\rm b}< 3\bar{\rho }_{\rm b} \, 2^{3(1+\phi)(\ell -\ell _{\rm base})}$ at all times, where ℓ_base = 5 is the refinement level corresponding to the base grid of size $128=2^{\ell _{\rm base}+2}$ , and we choose the scaling parameter to equal ϕ = −0.3 (value consistent with Safranek-Shrader et al. 2012). We insert a star particle when the hydrogen number density first exceeds 1000 cm⁻³; this happens at redshift z ≈ 19.7 and refinement level ℓ = 14, corresponding to cell size Δx ≈ 0.74 pc (physical).

At the end of the star's life and just prior to inserting the young supernova remnant, we force the simulation to refine further to level ℓ = 21 at the location of the star particle, corresponding to cell size Δx ≈ 0.006 pc (physical). This allows us to insert the supernova remnant well within the free-expansion phase, before the ejecta have started to interact with the circumstellar environment. Then, we disable automatic derefinement and allow the simulation to refine additional computational cells, up to the maximum level ℓ_max initially equal to 21, based on the standard second derivative test in FLASH, but in this case, the second derivative test is computed using the cell metallicity. This ensures that the entire region enriched with metals will be resolved at the refinement level ℓ_max. As the supernova remnant expands, we gradually and cautiously lower ℓ_max, and thus coarsen the resolution of the metal-enriched region, to allow us to use ever-longer time steps and integrate the simulation for 40 Myr past the insertion of the supernova.

2.4. The First Star and Ionizing Radiation Transfer

When the hydrogen number density first exceeds 1000 cm⁻³ near the center of the collapsing minihalo, we insert a collisionless star particle.⁴ The star is free to move in the combined gravitational potential of the dark matter and baryons but does not contribute to the potential, i.e., it is treated as massless. We calculate the transfer of the star's ionizing radiation and the response of the circumstellar medium to the associated photoheating in the H ii region as follows. We assume that it is an isotropic source of hydrogen-ionizing radiation and treat the radiation as monochromatic with photon energy hν = 16 eV and photon emission rate Q(H) = 6 × 10⁴⁹ photons s⁻¹. This ionizing photon rate is about three times as high as has been estimated for a 40 M_☉ metal-free star without mass loss, Q(H, 40 M_☉) ≈ 2 × 10⁴⁹ photons s⁻¹ (Schaerer 2002). We deliberately overestimate the ionizing luminosity to allow for the presence of additional, lower-mass, longer-lived stars also producing additional ionizing photons in a compact cluster around the primary 40 M_☉ star. Also, boosting the ionizing luminosity allows the H ii region to break out more easily at the moderate spatial resolution at which we insert the star particle and perform ionizing radiation transfer.

Our intention was to calibrate the monochromatic ionizing photon energy to a 40 M_☉ star, but the prescription we used to arrive at the adopted value of 16 eV suffers from lack of physical consistency. In retrospect, the inconsistency could have been avoided by carrying out photon-number-flux-weighted averaging of photoionization cross sections and heating rates as in, e.g., Appendices A.2 and A.3 of Pawlik et al. (2012), assuming a spectral energy distribution containing a single 40 M_☉ Pop III star with an effective temperature ∼8 × 10⁴ K (Schaerer 2002), as well as perhaps several lower-mass companion Pop III stars with lower effective temperatures. This would have resulted in higher average photoelectron energies than in our monochromatic treatment. The unintended consequence of our crude prescription is that in the simulation, photoionization will heat the gas only to a maximum temperature ∼14, 000 K, whereas for an effective photospheric temperature of, e.g., ∼10⁵ K (similar in moderate-mass and very massive metal-free stars), the gas at small distances ≲ 1–10 pc from the star should be heated to higher temperatures ∼30, 000 K, especially as the ionized gas density drops (see, e.g., Yoshida et al. 2007). An artifact of this is a potential underestimate of the amplitude of the outward-propagating pressure wave resulting from the intense photoionization heating in the H ii. We address this issue further in Section 3.1 below.

To carry out ionizing radiation transfer, we center a spherical coordinate system on the current position of the star particle and pixelize the angular coordinate using the HEALPix scheme of Górski et al. (2005) with N_pix = 3072 pixels. Within each angular pixel, we split the radial coordinate in the range 4 × 10²⁰ cm ⩽ r < 10²⁴ cm (comoving), corresponding to 0.7 pc ≲ r < 17 kpc (physical), into N_rad = 103 logarithmic radial bins. For each volume element E_{p, b} defined by an angular HEALPix pixel p and a radial bin b, and for each mesh cell E_c defined by the cell index c, we compute the volume of the intersection of the two elements V_{p, b; c} ≡ vol(E_{p, b}∩E_c) ⩾ 0. In general, only a fraction of the cell lies within E_{p, b}. In this case, we recursively split the cell into sub-cells via an octree subdivision procedure and compute the contribution of the subcells lying within E_{p, b} to the intersection volume; the recursion is continued until a desired accuracy is achieved. Armed with the intersection volumes, we compute the average case-B hydrogen recombination rate $\dot{n}_{{\rm rec},p,b}$ evaluated under the assumption of full ionization (as in the standard Strömgren calculation) within each E_{p, b} via

$\begin{equation} \dot{n}_{{\rm rec},p,b} = \frac{1}{V_{p,b}}\sum _{c} \alpha _{{\rm rec},B}(T_c) n_c^2 V_{p,b;c}, \end{equation} \tag{ 1 }$

where V_{p, b} = vol(E_{p, b}) ≈ ∑_cV_{p, b; c}. In Equation (1), T_c and n_c are, respectively, the temperature and density of hydrogen nuclei within the cell, and α_{rec, B}(T) is the hydrogen recombination coefficient. Within each bin, we compute the approximate integral of $4\pi r^2 \dot{n}_{{\rm rec},p}(r)$ along the radial direction (now replacing the bin index b with the continuous radial coordinate r) and determine the pixel-specific Strömgren radius R_{S, p} that solves the equation

$\begin{equation} Q ({\rm H}) = \int _0^{R_{{\rm S},p}} 4\pi r^2 \dot{n}_{{\rm rec},p}(r) dr. \end{equation} \tag{ 2 }$

We assume that no ionizing radiation penetrates to radii r ⩾ R_{S, p}. At smaller radii, we compute the local ionizing photon number flux (in units of photons cm⁻² s⁻¹) using

$\begin{equation} f_r(r) = \frac{1}{4\pi r^2} \left[Q ({\rm H}) - \int _0^{r} 4\pi r^2 \dot{n}_{{\rm rec},p}(r^{\prime }) dr^{\prime } \right]. \end{equation} \tag{ 3 }$

The local ionizing flux is then used to compute the local photoionization parameter for cell c lying within the pixel p with a center at a radius r_c via ξ_p(r) = 4πhνf_p(r_c)/n_c. Given a value of the photoionization parameter, the photoionization equilibrium temperature of the gas and the chemical abundances within the cell are evaluated via table lookup (see Section 2.2). The computation is fully parallel and is repeated at every hydrodynamic time step. We do not compute the radiation pressure from Thomson scattering and photoionizations. This radiative transfer scheme is not explicitly photon conserving and is not designed to correctly reproduce the ionization front kinematic, especially when an R-type ionization front is expected, but it is sufficient to reproduce the basic hydrodynamic response of the gas to photoionization heating.

2.5. Supernova and Metal Transport

We insert the supernova 3 Myr after inserting the star. The supernova is initialized in the free-expansion phase, when the remnant is about t − t_SN = 35 yr old, well before the freely expanding ejecta have started to interact with the primordial circumstellar medium. This early insertion is possible thanks to our ability to drastically increase the local mesh resolution at a specific time in the simulation—just prior to supernova explosion—and then gradually degrade the resolution as the remnant expands. The ejecta are initially cold and have a total mass M_SN = 40 M_☉ and a total kinetic energy E_SN = 10⁵¹ erg. The metal yield is assumed to be Z_SN = 0.1 and uniformly premixed in the ejecta, implying 4 M_☉ of metals. The velocity of the ejecta is linear in the distance from the center of the explosion, and the density is uniform. The initial radius of the remnant is chosen such that its radius is 10% of the radius of the remnant at which the swept-up mass equals the ejecta mass, or in this case, R_SN = 0.075 pc. The simulation carries out advection of the absolute metallicity Z by treating it as a passive mass scalar quantity. This advection is subject to undesirable numerical diffusion, which is unavoidable in Eulerian solvers of this type (see, e.g., Plewa & Müller 1999). Therefore, in parallel with passive scalar advection, we also insert a number N_part = 5 × 10⁵ of passive Lagrangian tracer particles, originally uniformly distributed in the ejecta. The simulation advances particle positions by integrating the first-order ordinary differential equation dX_i/dt = v(X_i), where X_i is the location of particle i and v(X_i) is the baryon fluid velocity, quadratically interpolated from the computational grid, at the location of the particle. We do not take into account the possibility that the supernova may be a source of molecules and dust. Both are undoubtedly produced in the free-expansion phase, but their survival of the reverse shock hinges on small-scale clumpiness of the ejecta and on other complex, incompletely understood physics (see, e.g., Cherchneff & Dwek 2009, 2010).

3. RESULTS

In Figure 1, we show the redshift evolution of the radius R₂₀₀ within which the average halo density equals 200 times the critical density, as well as the corresponding enclosed halo mass M₂₀₀. The star particle is inserted at z ≈ 19.7 and the supernova at z ≈ 19.5. The halo radius doubles from R₂₀₀ ≈ 160 pc to 300 pc in the course of the simulation, and the halo mass quadruples from M₂₀₀ ≈ 10⁶ M_☉ at star insertion to 4 × 10⁶ M_☉ at the end of the simulation. Figure 2 shows the density–temperature phase diagram for gas within 1 kpc (physical) from the center of the halo, immediately prior to the formation of the star particle, after the supernova explosion, and finally at the end of the simulation when the halo mass has quadrupled. Before stellar insertion, the thermodynamic structure of the halo is that of an archetypical minihalo in which gas cools by rovibrational transitions of molecular hydrogen. The thermodynamic structure at the end of the simulation is similar, which illustrates that the halo relatively quickly "recovers" from the radiative and kinetic impact of the first star and resumes central gravitational collapse. Photoionization and the supernova blast wave, which we discuss in detail in Sections 3.1 and 3.2, respectively, only temporarily disrupt the cooling and condensation of gas in the halo center. Then, in Section 3.3, we study the transport of the metals injected by the explosion.

**Figure 2.** Density–temperature phase plot for the gas within 1 kpc (physical) from the center of the minihalo at redshifts, from left to right, z = 19.7 (right before the insertion of the star), z = 19.5 (in the Sedov–Taylor phase of the supernova, 8 kyr after the explosion), and z = 17 (end of the simulation, 40 Myr after the explosion). By the end of the simulation, the gas heated by photoionization and the supernova has cooled and the halo has resumed central gravitational collapse, albeit at four times higher halo mass and now a nonzero metallicity. The dashed line shows the CMB temperature.
Download figure:
Standard image High-resolution image

3.1. The H ii Region

The photoionization raises the temperature of the interior of the H ii region to ≈10⁴ K with a maximum temperature ≈1.4 × 10⁴ K; the latter temperature ceiling is an artifact of our photoionization prescription (see Section 2.4). The increase in gas pressure drives an outflow in the ionized gas, thus lowering the central density. In Figure 3, we show the star-centered, spherically averaged gas density and spherically averaged, mass-weighted radial velocity profile just before inserting the star and 3 Myr later, just prior to inserting the supernova. During this interval, the spatial resolution in the vicinity of the star is Δx ≈ 0.75 pc. We separately examine the density and velocity profiles of all gas and of the ionized gas only. The ionized gas has a maximum outward velocity of ∼9 km s⁻¹, which can be compared to the maximum gas inflow velocity ∼ − 7 km s⁻¹ before photoionization. The central density drops from ∼2000 cm⁻³ to ∼4 cm⁻³. Interestingly, at radii as small as 2–3 pc, the densest gas clumps (n ∼ 10–1000 cm⁻³) resist photoionization and do not acquire the outward velocity of the ionized gas.

**Figure 3.** Spherically averaged density (upper panel) and radial velocity (lower panel) profiles for all baryons (blue line) and the ionized gas (red line), before the insertion of the star at z = 19.7 (dashed lines) and before the insertion of the supernova at z = 19.5 (solid lines). For reference, we also show the density profile of gravitational mass density (including dark matter) in the units of the proton mass (black line). The neutral gas density peak located inside the H ii region is a clump of dense neutral gas that has withstood photoionization. The gas with positive (negative) radial velocity is flowing away from (toward) the star.
Download figure:
Standard image High-resolution image

An examination of the three-dimensional geometry of the ionized and neutral phases in Figure 4 shows that the sustained neutral (i.e., self-shielding) clumps are associated with the photoablated terminus of dense gas filaments feeding into the halo from the cosmic web. This outcome was already observed by Abel et al. (2007) in a similar simulation of a primordial H ii region with more accurate radiation transfer. It is worth stressing that our simulation may not properly resolve the process of the ionization of dense cloud cores because the method is not explicitly photon conserving, and because it may not resolve the shock that is expected to precede the D-type ionization front in the cold gas (we find that the ionized gas near the D-type front is slightly overpressured relative to the cold gas). For these reasons, we are not certain that the simulation has converged in terms of the rate of photoablation of dense neutral clumps, but qualitative consistency with previous work is reassuring. Continuous photoablation presents a source of fresh, dense ionized gas in the vicinity of the star. This, in addition to the temperature ceiling mentioned above, explains why the central ionized gas density in our simulation remains higher than the value ∼0.1 cm⁻³ characteristically observed in one-dimensional models.

We do not observe a clear indication of the expected outward-propagating pressure wave and the associated shock wave, which would have manifested as a dense, supersonically propagating shell. Such a shell could become important at later times as it would modify the dynamics of the supernova blast wave; e.g., it could drive a reverse shock into the supernova ejecta, which could trigger hydrodynamical instability and mixing of the ejecta with the swept-up circumstellar medium (Whalen et al. 2008). The absence of the pressure wave can be attributed to our approximate treatment of ionizing radiation transfer where the ionization front speed is not properly constrained by the finite ionizing photon consumption rate and is resolution dependent, and also because of the artificial temperature ceiling. However, the ionized gas does acquire a positive radial velocity out to radii ∼120 pc, consistent with expectations for our adopted source luminosity.

3.2. Blast-wave Evolution and Fallback

Because of the low energy and compactness of the remnant, the cooling by inverse Compton scattering off the CMB is relatively unimportant, different from the case of ultra-energetic Pop III supernovae, where this cooling channel dominates the early evolution (Bromm et al. 2003). The blast wave also largely avoids the dense neutral clumps resisting photoionization, and they remain a repository of low-entropy gas in the halo. This is easily seen in Figure 5, which shows mass-weighted projected density, temperature, and metallicity 8.5 Myr after the explosion. Outside of the filaments, after 40 Myr, the effects of photoionization and the blast wave have largely been erased; the thermal phase structure within the halo is similar to that of a higher-mass unperturbed minihalo. The appearance of warmer (≲ 10⁴ K) gas at low densities reflects the higher virial mass of the halo at the end of the simulation and the residual entropy left by the fossil H ii region. The presence of metals at local metallicities Z ≳ 10⁻⁵ to 10⁻³ Z_☉ and an enhanced abundance of hydrogen–deuteride now facilitate more rapid cooling toward low temperatures. The enhanced cooling is now especially apparent at high densities n_H ≳ 100 cm⁻³ that are found in the central region. The dense gas is insufficiently well resolved at the end of the simulation; the resolution limit hinders compression of the fluid to densities n_H ≳ 1000 cm⁻³. Nevertheless, it is clear that the dense gas resides at temperatures a factor of a few lower than the gas at the same densities prior to the insertion of the ionizing source. We do not observe the anticipated gas cooling to the CMB temperature, through either metal line cooling or the cooling by hydrogen–deuteride (e.g., Johnson & Bromm 2006), because we enforced an artificial temperature floor at 200 K.⁵

**Figure 5.** Baryon mass-weighted projection of hydrogen density, temperature, and fluid-based metallicity in units of the solar metallicity 8.5 Myr after the explosion. The circle shows the virial radius of the halo, which is here 180 pc (physical).
Download figure:
Standard image High-resolution image

An insight into the rate of recovery of the halo can be obtained from Figure 6, which shows the Lagrangian radii that track the radial location of a grid of mass coordinates of the spherically averaged baryon and supernova ejecta mass as a function of time from the explosion. In this period, we gradually coarsen the spatial resolution from Δx ≈ 0.006 pc, a high value achieved by telescopic AMR at the location of the explosion, to Δx ≈ 14 pc at the end of the simulation, 40 Myr later. It is worth keeping in mind that the underlying flow is complex and contains simultaneous inflows and outflows that are averaged over; each of the radii reflects only the net displacement of the mass coordinate. The baryonic Lagrangian radii show expansion, driven by the supernova blast wave, starting at ≈20 kyr at the innermost radius shown, and appearing at progressively larger mass coordinates until outward-directed motion becomes evident near the virial radius at ≈2 Myr after the explosion. Each of the radii corresponding to mass coordinates inside the halo, however, shows a reversal from a net outflow to a net inflow. The reversal commences at ≈400 kyr at the innermost radius shown and reaches progressively larger radii until ≈20 Myr, when the flow at the virial radius turns into an inflow. Outside of the virial radius, cosmic infall persists largely unabated throughout the remnant's history.

An examination of the evolution of the actual three-dimensional structure of the baryonic flow (see the projection in Figure 5) reveals that the prompt reversal can only in part be attributed to a fallback of the mass that the blast wave has swept up and to which it has transferred momentum. Instead, the filamentary inflow of dense, cold, molecular gas from the nearby cosmic web directly into the center of the minihalo contributes a large fraction of the mass inflow rate. The baryonic filaments are seen diagonally from top left to bottom right in Figure 7, which shows a slice of the density field. Because the filaments survive photoionization and the explosion largely intact, they can sustain a net inflow into the center of the halo even when a majority of the ejecta and the swept-up mass are moving outward.

The Lagrangian radii for the supernova ejecta in Figure 6, right panel, show first signs of interaction of the ejecta with the circumstellar medium around 1–10 kyr after the explosion. What follows is a similar picture of expansion eventually reversing into fallback. Fractions of (3.125%, 6.25%, 12.5%, 25%, 50%) of the ejecta mass turn around and start falling back at around (0.2, 1.0, 4.0, 8.0, 30.0) Myr after the explosion. The corresponding turnaround radii are (3, 8, 30, 60, 200) pc. The outer ∼50% of the ejecta passes the virial radius and continues to travel outward until the end of the simulation. The net outflow of ejecta at r ≳ R₂₀₀ can be contrasted with the net inflow of baryons at the same radii, underscoring the three-dimensional nature of the flow and an inadequacy of one-dimensional integrations in treating the long-term dynamics of supernova remnants in cosmic minihalos. Spherical symmetry implicit in one-dimensional integrations forces the ejecta and the swept-up circumstellar medium into a spherical thin shell at late times. It does not allow for the presence of the dense, low-entropy clouds that enter the halo by infall from the filaments of the cosmic web. The expanding shell may pass and leave behind such clouds. The spherical symmetry also does not allow radial segregation of the ejecta as a result of the instability of the decelerating thin shell.

Rayleigh–Taylor fingering is evident in Figure 5 at 8.5 Myr after the explosion. The long-term evolution of the fragmented shell is seen in Figure 7, left panel, which shows the ejecta tracer particles in projection, superimposed on a slice of the density field through the center of the halo. It is also clear in these figures that the blast wave has expanded biconically perpendicular to the sheet-like baryonic overdensity deriving from the cosmic web; this pristine gas, unpolluted by the ejecta, is still able to stream into the halo center. In the right panel of Figure 7, showing the central 360 pc, pristine baryonic streams arrive diagonally into the central few tens of parsecs, where they form a two-armed spiral and join a self-gravitating central core. At this point in the simulation, the spatial resolution is insufficient to resolve the internal structure of the core.

Strikingly, we observe filamentary accretion of ejecta-enriched gas perpendicular to the pristine sheet. This is the return of the ejecta from the supernova remnant back into the center of the halo. A closer examination of the time evolution of the ejecta kinematics shows that in the snowplow phase, the ejecta are compressed into a thin shell, as expected. The thickness of the shell is itself not well resolved; thus, we are not in the position to explore its susceptibility to direct gravitational fragmentation (e.g., Mackey et al. 2003; Salvaterra et al. 2004; Machida et al. 2005; Chiaki et al. 2012). The Rayleigh–Taylor instability buckles and redistributes momentum in the shell. The highest momentum fragments continue semi-ballistically outward, consistent with the observation by Whalen et al. (2008) that the blast wave of even a moderate 10⁵¹ erg explosion has enough momentum to escape a standard minihalo. The low-momentum fragments, however, do not have sufficient momentum to escape; they turn around, fall back, and eventually join the central unresolved core. The falling ejecta make thin multiply folded streams and filaments, indicating what seems to be a certain lack of mixing with the halo gas. The pristine gas accreting through the filaments of the cosmic web and the ejecta-enriched gas accreting from the supernova remnant meet in the central core. If the core is turbulent—poor resolution at this stage of the simulation would not allow us to detect such turbulence—turbulent mixing of the ejecta, now containing some swept-up pristine gas, with a much larger mass of pristine gas would ensue.

3.3. Metal Enrichment

In Figure 8, we show the net rates of baryon and metal outflow from or inflow into the central 20 pc of the halo. Some baryonic inflow is already evident early on, after 50 kyr from the explosion, and after ≈1.5 Myr, baryons transition to permanent inflow. After ≈5 Myr, the inflow rate settles at $\dot{M}_{\rm b}\approx 0.002\,M_\odot \,\hbox{yr}^{-1}$ . The metal flow transitions from outflow to inflow at ≈4 Myr after the explosion and fluctuates in the range $\dot{M}_Z\sim (0.5\hbox{--}5)\times 10^{-7}\,M_\odot \,\hbox{yr}^{-1}$ . The average asymptotic absolute metallicity of the inflowing material is

$\begin{equation} Z_{\rm inflow}= \frac{\dot{M}_Z}{\dot{M}_{\rm b}} \sim 10^{-4}\left(\frac{Z_{\rm SN}}{0.1}\right) \sim 0.005 \,Z_\odot \left(\frac{Z_{\rm SN}}{0.1}\right), \end{equation} \tag{ 4 }$

where we are referring to the total metal mass ignoring fractional abundances of individual metal species.

**Figure 8.** Metal (lower curve; $\dot{M}_Z$ ) and total baryon (upper curve; $\dot{M}_{\rm b}$ ) net mass flow rate through a sphere of radius 20 pc centered on the gravitational potential minimum. Dashed lines indicate outflows and solid lines inflows. The total mass of ejecta in the simulation is M_SN = 40 M_☉, and the metal mass is *M_Z* = Z_SNM_SN = 4 M_☉. Net outflow reverses into an inflow earlier in the baryons because of the presence of cold filaments delivering metal-free gas from the cosmic web into the halo center.
Download figure:
Standard image High-resolution image

**Figure 8.** Metal (lower curve; $\dot{M}_Z$ ) and total baryon (upper curve; $\dot{M}_{\rm b}$ ) net mass flow rate through a sphere of radius 20 pc centered on the gravitational potential minimum. Dashed lines indicate outflows and solid lines inflows. The total mass of ejecta in the simulation is M_SN = 40 M_☉, and the metal mass is *M_Z* = Z_SNM_SN = 4 M_☉. Net outflow reverses into an inflow earlier in the baryons because of the presence of cold filaments delivering metal-free gas from the cosmic web into the halo center.
Download figure:
Standard image High-resolution image

We proceed to examine the evolution of the metal-mass-weighted metallicity probability density function (PDF). The metallicity PDF is computed in three different ways. In the first, "fluid–fluid" method, the metallicity is obtained by summing the amplitudes of the passive scalars defining metal abundances on the computational grid. The metal mass at metallicity Z is then calculated by computing a volume integral of Zρ over the cells with metallicities in a narrow logarithmic bin containing Z. The shortcoming of this method is that it yields, as a consequence of a well-known shortcoming of passive scalar transport schemes, nonvanishing, small metallicities even in many cells that ejecta could not have reached on hydrodynamic grounds.⁶ With this in mind, the first method is more accurate where the metallicity is near its peak value and is highly inaccurate elsewhere. Our remaining two methods for computing the metallicity PDF are designed to alleviate the impact of spurious diffusion.

In the second, "fluid–particle" method, the metallicity Z of a computational cell is calculated in the same way as in the first method, but the metal mass is calculated by summing the mass of the ejecta tracer particles in the cell and multiplying by the ejecta absolute metallicity Z_SN = 0.1. The tracer particle mass inside a cell is calculated with the cloud-in-cell method, with the weighting function taken to correspond to spreading the particle mass over a cubical kernel identical to the host cell, but centered on the particle. In the second method, only the cells having nonvanishing overlap with at least one tracer particle cloud-in-cell kernel contribute to the metallicity PDF. It should be kept in mind, however, that a metal particle displaced from the center of its host cell contributes mass to the neighboring cell that may be separated from the host cell by a shock transition or a contact discontinuity.

In the third, "particle–particle" method, the metallicity inside the cell is itself computed by summing up the metal mass inside the cell (again, as in the second method, computed by adding up tracer particle mass inside the cell and multiplying by Z_SN) and dividing this by the sum of the metal mass and the fluid hydrogen and helium masses in the cell. Since the combined abundance of hydrogen and helium is near unity, the third method completely ignores the fluid metallicity as defined by the passive scalar advected on the computational grid. Ejecta particles as tracers of metallicity obviously do not suffer from the same spurious diffusion as the fluid metallicity, but they are affected by systematic issues of their own. In convergent flows such as near the insufficiently well resolved snowplow shell, metal particles can get trapped in discontinuous flow structures, which can produce a spurious increase of local metallicity as computed by the third method. Thus, while the fluid metallicity is bounded from above by its maximum value, Z_SN = 0.1, at the point of injection, the particle-based metallicity can exceed this value and even approach unity.

In Figure 9, we show the metal-mass-weighted metallicity PDF at four different times: first, in the Sedov–Taylor phase, and then at approximately 1, 13, and 40 Myr after the explosion. It is immediately clear that already at early times, the metallicity PDF exhibits a tail dM_Z/dln Z∝Z extending toward very low metallicities, where differences between the PDFs computed with the three methods became more severe. Spurious diffusion across compositional discontinuities in the fluid-based methods and cloud-in-cell smearing in particle-based methods both contribute to the tail. We also notice that the particle–particle method has developed a tail extending into the forbidden region, Z > Z_SN, a clear signature of spurious particle trapping.

As the supernova remnant ages, metal mixing in the complex interior flow, which is numerical in the simulation but can be facilitated by turbulence in nature (e.g., Pan et al. 2011), broadens the metallicity peak from its initial location at Z ≈ Z_SN to an approximate range 10⁻⁵ ≲ Z ≲ Z_max(t), where Z_max(t) is a time- and PDF-extraction-method-dependent maximum absolute metallicity. For fluid-based metallicity PDFs, we find that Z_max ≈ (0.07, 0.02, 0.002) at t − t_SN ≈ (1, 13, 40) Myr. We also notice a prominent narrow metallicity peak, built from cells belonging to the central, unresolved core that is accreting both from the pristine filaments of the cosmic web and from the ejecta-enriched supernova remnant fallback. The peak is migrating slowly toward lower metallicities, covering Z_core ≈ (2, 0.7, 0.3) × 10⁻⁴ in the same three epochs. The metallicity of the unresolved core is consistent with the average metallicity of the accreting fluid given in Equation (4).

Finally, we turn to the energetics of the metal-enriched material. In Figure 10, we show the evolution of the total kinetic and internal (thermal) energy of fluid elements having absolute metallicity Z ⩾ 10⁻⁶. At 1 kyr after the explosion, the kinetic energy decreases below its initial value when the ejecta start interacting with the ambient medium and the well-resolved reverse shock facilitates kinetic-to-internal conversion. Adiabatic and radiative cooling both contribute to the rapid decrease of internal energy at 0.01–1.0 Myr after the explosion. The renewed rise of both energies in the old remnant reflects the fallback of ejecta-enriched fluid into the halo center.

4. DISCUSSION

Existing treatments of metal injection by Pop III supernovae have followed the excursion of the metal-endowed ejecta into the cosmic environment and their possible fallback into dark matter halos, but usually with one-dimensional simulations (Kitayama & Yoshida 2005; Whalen et al. 2008) or relatively low resolution, three-dimensional, grid-based or smoothed particle hydrodynamic simulations (Nakasato & Shigeyama 2000; Bromm et al. 2003; Greif et al. 2007; Wise & Abel 2008b). These treatments point to a prompt and widespread dispersal of supernova ejecta into the primordial IGM if one assumes that the supernovae were extraordinarily energetic, e.g., consistent with the popular hypothesis that the Pop III stars exploded as PISNe. However, in particularly high mass halos or those hosting moderate-mass stars, the detailed dynamics of the ejecta is sensitive to the pre-supernova evolution of the H ii region. With one-dimensional hydrodynamic simulations, Whalen et al. (2008) showed that the H ii region of a moderate-mass star at the center of the minihalo may remain confined within the minihalo. They calculated that the blast wave of a 10⁵¹ erg supernova in ≲ 10⁶ M_☉ minihalos does, after ∼10 Myr, reach the virial radius of the halo and concluded that this results in the "disruption" of the baryonic component of the halo. They also considered high-mass ∼10⁷ M_☉ minihalos in which a strong photodissociating background had prevented star formation (e.g., O'Shea & Norman 2008) and found that there the H ii region may not break out and the supernova blast wave may be confined at very small radii. They speculated that the subsequent prompt fallback of metal-enriched supernova ejecta into the center of the halo could lead to the formation of a globular-cluster-like stellar system, similar to the suggestion of Wada & Venkatesan (2003) that a relatively small number of supernovae exploding in a more massive ∼10⁸ M_☉ cosmic halo could also lead to a fallback and prompt star cluster formation. This is a particularly interesting possibility, in view of the unexplained origin of metal-poor globular clusters (Brodie & Strader 2006): do they form in the centers of their own low-mass dark matter halos (e.g., Peebles 1984; Rosenblatt et al. 1988; Bromm & Clarke 2002; Boley et al. 2009; Griffen et al. 2010), or do they form off-center in the merging and fragmentation of more massive protogalactic disks (e.g., Kravtsov & Gnedin 2005, and references therein)?

It is important to realize that the fate of the supernova ejecta may be influenced by dynamical processes operating on relatively small spatial scales and strongly violating spherical symmetry. In addition to the likely intrinsic asphericity of the explosion itself, one such process is the Rayleigh–Taylor instability of the decelerating contact discontinuity after the ejecta have started to interact with the ambient medium. The ambient medium into which the supernova expands is bound to be inhomogeneous, especially if the luminosity of the progenitor stars and other hot stars that have formed with it is insufficient to completely photoionize the densest clumps remaining from the progenitor gas cloud. On the somewhat larger ∼100 pc scales of the host minihalo, a pronounced departure from spherical symmetry will arise from the filamentary geometry of the cosmic web. Such asymmetry can interact with the expanding blast wave. Further, as an aging supernova remnant forms a thin, momentum-conserving snowplow shell, the decelerating shell is itself unstable to the Rayleigh–Taylor instability and forms fingers. Finally, while the blast-wave instabilities lead to partial entrainment of the ejecta and the shocked primordial gas, any ejecta-enriched gas falling back toward the center of the halo, where it may engage in a new round of star formation, may not necessarily be chemically homogeneous (full chemical mixing requires at least a few turbulent eddy turnover times, which may not be available in the blast wave). Most of these phenomena remain to be explored in the context of Pop III stars; the results we have presented here are one of the first steps in that direction.

4.1. Metal-enriched Star Formation: Continuous, Bursty, or Self-limiting?

The prompt resumption of baryonic infall into the center of the cosmic minihalo suggests that star formation could resume on similarly short, ∼1–5 Myr, timescales. In the simulation presented here, the metals that the supernova has synthesized begin returning to the halo center ∼4 Myr after the explosion. The mean metallicity of the inflow, averaged over the primordial streams and the metal-enriched collapsing supernova remnant streams converging in the center of the halo, is ∼0.005 Z_☉ assuming a net supernova metal yield of 4 M_☉. If the metal-rich gas successfully mixes with the primordial gas, then, according to a prevailing belief, the metallicity of the combined streams is well sufficient to ensure that new star formation in this gas should be producing normal, low-metallicity Pop II stars. The Pop III to II transition cannot be ascertained in the present simulation because we do not adequately resolve the cold gas at densities n_H ≳ 100 cm⁻³. This gas belongs to the central core and is colder than primordial gas at the same densities preceding the insertion of the first star. It is at unresolved densities that the metallicity of the gas strongly affects the thermodynamic evolution and thus governs the outcome of fragmentation (e.g., Omukai et al. 2005; Santoro & Shull 2006; Schneider et al. 2006; Smith et al. 2009; Safranek-Shrader et al. 2010, assuming low metallicities, Z ≲ 0.01 Z_☉).

The gaseous core in the center of the minihalo that is fed by metal-enriched streams can be compared to the star-forming core of a giant molecular cloud (GMC). We attempt such comparison with an eye on the overall hydrodynamics of the gas flow, putting aside any differences in chemical composition, dust content, CMB temperature, etc., in the two systems. Dense GMC cores grow by accreting from larger-scale, lower-density cloud environments that are supersonically turbulent. Smaller cores and their embedded protostar clusters merge with each other to form larger associations. Star formation in the cores is thought to be self-regulated, with the "feedback" from ongoing star formation controlling the amount of star-forming gas available and eventually leading to GMC dissolution (e.g., Murray et al. 2010; Murray 2011; Goldbaum et al. 2011; Lopez et al. 2011). The feedback entails, e.g., the pressure of the H ii regions and the pressure of direct stellar radiation. While the baryonic mass of a cosmic minihalo is perhaps smaller than that of a GMC, the centrally concentrated underlying dark matter halo renders the gravitational potential depth in a minihalo's central gaseous core similar to that in a GMC's star-forming core. The underlying dark matter gravitational potential may make a cosmic minihalo less susceptible to star-formation-induced dissolution than a GMC, assuming that the two systems are forming stars at the same rate.

GMCs are thought to sustain star formation over tens of millions of years, in spite of the H ii regions, supernovae, and other feedback activity that commences quickly following the formation of the first massive stars. The feedback's impact preceding the eventual dissolution of the cloud complex may be limited to reducing the average star formation rate. We do not know whether star formation in metal-enriched minihalos and their descendent ≳ 10⁷ M_☉ halos can be sustained as in GMCs in spite of the feedback, or whether the feedback from second-generation stars is so effective as to evacuate baryons from the halo and prevent further star formation on cosmological timescales (see, e.g., Mori et al. 2002; Frebel & Bromm 2012). It seems that these two scenarios can be evaluated only with attention to the detailed, small-scale mechanics of star formation and its feedback. A corollary of the first scenario would be the presence of multiple, chemically distinct stellar populations in the stellar system that forms this way.

The statistics of low-luminosity stellar systems in the Local Group in conjunction with theoretical modeling of structure formation in a ΛCDM universe suggests that star formation must have, to some extent, been suppressed in low-mass cosmic halos (e.g., Susa & Umemura 2004a; Ricotti & Gnedin 2005; Gnedin & Kravtsov 2006; Moore et al. 2006; Madau et al. 2008; Bovill & Ricotti 2009, 2011; Koposov et al. 2009; Muñoz et al. 2009; Salvadori & Ferrara 2009; Busha et al. 2010; Gao et al. 2010; Griffen et al. 2010; Macciò et al. 2010; Okamoto et al. 2010; Font et al. 2011; Lunnan et al. 2012; Rashkov et al. 2012). This suppression, however, could have arisen both internally, as in GMCs, and externally, where the suppression is driven by the heating of the IGM by the background ionizing radiation coming from other, more massive but less common halos (e.g., Thoul & Weinberg 1996; Kepner et al. 1997; Barkana & Loeb 1999; Shapiro et al. 2004; Susa & Umemura 2004b; Mesinger et al. 2006; Mesinger & Dijkstra 2008; Okamoto et al. 2008). The heating could have shut off the inflow of fresh gas into the halo, eventually stunting star formation. The prevalence of either suppression mechanism as a function of halo properties cannot be understood without first understanding the detailed mechanics of star formation and feedback following the initial fallback of Pop III supernova ejecta.

4.2. Uncertainties Related to the Nature of First Stars

A source of uncertainty regarding the formation and nature of the first metal-enriched star clusters is the incomplete understanding of the precise character of metal-free Pop III stars and of the energetics and nucleosynthetic output of the associated supernovae. Theoretical investigations of Pop III star formation have not yet converged, and the precise functional form of the resulting initial mass function is not yet known. The evolution of metal-free stars is very sensitive to their poorly understood internal rotational structure starting with the protostellar phase (see, e.g., Stacy et al. 2011, 2012b) and the potential presence of a close companion. Similarly, the explosion mechanism in metal-poor stars is not known. Each of the suggested mechanisms, including the standard delayed-neutrino mechanism hypothesized to drive Type II explosions, the pair-instability and pulsational pair-instability mechanisms expected in relatively massive metal-free stars (see, e.g., Heger & Woosley 2002; Chatzopoulos & Wheeler 2012; Yoon et al. 2012, and references therein), and the more exotic "hypernova"-type mechanisms involving a (possibly magnetized) outflow from a central compact object, implies its own characteristic explosion energy and nucleosynthetic imprint.

Another aspect of the final stellar demise that should be very sensitive to both the stellar mass and rotation is the nature of the compact remnant, if any, that it leaves behind. While the pair instability disperses the star completely, leaving no remnant, core collapse in stars with initial masses ≳ 15–25 M_☉ (the precise threshold mass depending on unknown aspects of the explosion; see, e.g., Heger et al. 2003; Zhang et al. 2008) leaves behind a black hole, either via direct collapse into one or by fallback of the ejecta onto the neutron star. The black hole may be produced with an initial kick that would eject it from the dark matter halo (Whalen & Fryer 2012), but if the kick is small and the black hole is retained near the gravitational center, it could accrete the gas collecting in this region. The energy liberated in this accretion could be another, potentially significant source of feedback. The accretion rate and the luminosity of this source would be subject to complex radiation–hydrodynamic couplings that are only understood under the most idealized assumptions, such as in an absence of shadowing by or radiation trapping within the innermost accretion flow (see, e.g., Milosavljević et al. 2009a, 2009b; Park & Ricotti 2011, 2012; Li 2011). To date, only very few cosmological simulations investigating the impact of the radiation emitted by accreting first-star remnants on hydrodynamics in the host halo have been carried out (e.g., Kuhlen & Madau 2005; Alvarez et al. 2009; Johnson et al. 2011; Jeon et al. 2012), and only a partial picture of the feedback on second-generation star formation is available. Jeon et al. (2012) show that X-rays from the remnant impact the cosmic neighborhood in complex ways but do not entirely prevent (and can actually promote) gas cooling and star formation. Much remains to be understood about the role of compact remnants in the evolution of the first galaxies.

Many investigations of the formation of first stars and galaxies and their role in reionization have explored implications of the hypothesis that Pop III stars had high masses, explosion energies, and metal yields. We find that the opposite limit of moderate masses, explosion energies, and metal yields might imply sharply divergent outcomes, not only for the parent cosmic minihalos but also for the overall pace of transformation of the early universe. A related issue is the nature of the first galaxies (for a review, see Bromm & Yoshida 2011). Their assembly process is influenced by the feedback from the Pop III stars that formed in the minihalo progenitor systems, and this feedback in turn sensitively depends on the mass scale of the first stars. In the case of extremely massive stars, negative feedback is very disruptive, completely sterilizing the minihalos after only one episode of star formation. Only after ∼10⁸ yr did the hot gas sufficiently cool to enable the recollapse into a much more massive dark matter halo, where a second-generation starburst would be triggered (Wise & Abel 2008b; Greif et al. 2010). This separation into two distinct stages of star formation would be smoothed out in the case studied here. Indeed, the first galaxies may then already have commenced their assembly in massive minihalos, as opposed to the atomic cooling halos (e.g., Oh & Haiman 2002) implicated in the previous scenario.

5. CONCLUSIONS

We have carried out a cosmological hydrodynamical simulation designed to investigate the evolution of a cosmic minihalo in the aftermath of the formation of the first, Pop III star, assuming that the star has a moderate mass and explodes as a moderate-energy Type II-like supernova, consistent with the recent downward revision of Pop III stellar mass estimates (Stacy et al. 2010, 2012a; Clark et al. 2011; Greif et al. 2011, 2012; Hosokawa et al. 2011). We analyzed the dynamics of supernova ejecta and gas flow inside the minihalo. These can be compared to similar studies involving massive Pop III stars exploding as ultra-energetic PISNe. Our main conclusions are as follows.

The moderate-mass star and the moderate-energy supernova it produces inflict significantly smaller damage to the host minihalo than that inflicted by a massive star exploding as an ultra-energetic supernova. The star only partially photoionizes the host halo. The densest gas flowing from the filaments of the cosmic web into the star's parent cloud remains neutral and survives the supernova blast wave. This dense gas resumes accretion into the center of the minihalo only ∼1 Myr after the explosion. After ∼20 Myr, all spherically averaged baryonic mass coordinates are moving inward and the accretion rate at the center of the halo is ∼0.002 M_☉ yr⁻¹.

Following instability of the blast wave during the snowplow phase, a fraction of the supernova ejecta starts falling back to the center of the halo. The ejecta fallback accretion rate reaches a steady state after ∼5 Myr. The accreting ejecta are incompletely mixed with the primordial gas and are confined in thin sheets and filaments. Less than half of the ejecta escape the virial radius; the escaping ejecta can be traced back to the outermost Rayleigh–Taylor fingers.

The average metallicity of the accreting matter is in the range 0.001–0.01 Z_☉. This result depends on the assumed supernova energy and metal yield, but in this study, we did not pursue parameter space exploration. The metallicity of the accreted gas is sufficient to ensure that new stars forming in the central core, unresolved in the simulation, will have even lower characteristic masses.

These results bring into focus metal-enriched star formation in cosmic minihalos and their immediate descendent halos. The character of post-supernova metal transport and gravitational fragmentation of the metal-enriched gas and of the impact of subsequent star formation and supernovae are both worth addressing in follow-up investigations. The early universe likely exhibited a range of different explosion settings, from modest to ultra-energetic. It will be intriguing to see how these high-redshift systems are matched with possible local fossils, such as the ultra-faint dwarf galaxies (e.g., Bovill & Ricotti 2009; Salvadori & Ferrara 2009; Tumlinson 2010) and metal-poor globular clusters (e.g., Brodie & Strader 2006). From what we have already learned, Pop III supernova feedback appears to play a key role in shaping early cosmic history and the nature of relic stellar systems.

We gratefully acknowledge stimulating discussions with J. Scalo and helpful technical consultation with A. Dubey, C. Daley, T. Plewa, P. Ricker, and M. Ruszkowski. We thank R. Banerjee for providing part of the software used to visualize this work. J.S.R. was supported in part by the John W. Cox Endowment for the Advanced Studies in Astronomy. This research was supported by NSF grants AST-0708795 and AST-1009928 and NASA ATFP grant NNX09AJ33G. The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing computing and visualization resources that have contributed to the research results reported within this paper. The software used in this work was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago.

CONFINED POPULATION III ENRICHMENT AND THE PROSPECTS FOR PROMPT SECOND-GENERATION STAR FORMATION

Article metrics

Permissions

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION

2. SIMULATING A SUPERNOVA IN A MINIHALO