THE FIRST GALAXIES: CHEMICAL ENRICHMENT, MIXING, AND STAR FORMATION

We perform our simulations in a cosmological box of size 1 Mpc³ (comoving), initialized at z = 99 with a fluctuation power spectrum corresponding to a concordance Λ cold dark matter (ΛCDM) cosmology with matter density Ω_m = 1 − Ω_Λ = 0.3, baryon density Ω_b = 0.04, Hubble parameter h = H₀/100 km s⁻¹Mpc⁻¹ = 0.7, where H₀ is the present Hubble expansion rate, spectral index n_s = 1.0, and a normalization σ₈ = 0.9 (Spergel et al. 2003). We take the chemical evolution of the gas into account by following the abundances of H, H⁺, H⁻, H₂, H⁺₂, He, He⁺, He⁺⁺, and e⁻, as well as the three deuterium species D, D⁺, and HD. We include all relevant cooling mechanisms, i.e., H and He collisional ionisation, excitation and recombination cooling, bremsstrahlung, inverse Compton cooling, and collisional excitation cooling via H₂ and HD (Glover & Jappsen 2007). We explicitly include H₂ cooling via collisions with protons and electrons, which is important for the chemical and thermal evolution of gas with a significant fractional ionization (Glover & Abel 2008).

In a preliminary run with 64³ gas and dark matter (DM) particles, we locate the formation site of the first ≃10⁸ M_☉ halo. This object is massive enough to cool through Lyα emission and fulfil our prescription for a "first galaxy" (Greif et al. 2008). We subsequently carry out a standard hierarchical zoom-in procedure to achieve high mass resolution in the comoving volume of the galaxy (e.g., Navarro & White 1994; Tormen et al. 1997; Gao et al. 2005). We apply three consecutive levels of refinement, such that a single parent particle is replaced by a maximum of 512 child particles. The particle masses for DM and gas in the highest resolution region containing the comoving volume of the galaxy with an extent of ≃300 kpc (comoving) are ≃33 M_☉ and ≃5 M_☉, respectively. The mass resolution is ≃400 M_☉, so that the Jeans mass at a density of n_H = 100 cm⁻³ and temperatures near the cooling threshold set by the cosmic microwave background (CMB) is well resolved.

The first minihalo collapses at a redshift z ≃ 30 near the location of the nascent galaxy. Once its central density exceeds n_H = 100 cm⁻³, we assume that a single Pop III star with 200 M_☉ forms, near the middle of plausible values determined in numerical studies (O'Shea & Norman 2007). However, for the purpose of determining the properties of the ensuing H ii region, we use the ionization, dissociation, and heating rates for a 50 M_☉ Pop III star (Schaerer 2002), which limits the severity of the radiative feedback. We note that the SN remnant will likely propagate somewhat further in this case, since pressure equilibrium with the ambient medium will be achieved at later times (see also Greif et al. 2007). On the other hand, the photoheating from nearby stars may mitigate this effect and lead to only minor deviations from the true metal distribution. The corresponding surface temperature, luminosity, and stellar lifetime may be found in Schaerer (2002), while details of the ray-tracing algorithm are described in Greif et al. (2009b). Two key assumptions of this model are that the gas is optically thin to molecule-dissociating radiation, and that the effects of stellar evolution are negligible, which would otherwise lead to a decrease in the number of ionizing photons emitted at the end of the main sequence (Marigo et al. 2001; Schaerer 2002).

At the end of its brief lifetime, the progenitor is predicted to explode as a PISN, ejecting of order 100 M_☉ in metals and 10⁵² erg of kinetic energy (Heger & Woosley 2002). In the simulation, we assign the above energy and metal content to the central 200 M_☉ of gas as described in Greif et al. (2007). The explosion is initialized at the end of the free expansion phase, where the momentum of the swept-up mass leads to a transition to the Sedov–Taylor phase. This simplified treatment is roughly valid as long as the central density is low, which is established by the previous photoheating from the star (e.g., Johnson et al. 2007). All particles within the initial ejecta are assigned a metal yield of 50%, consisting of C, O, and Si and relative abundances corresponding to a solar composition. Although this is somewhat inconsistent with the results of Heger & Woosley (2002), who show that the actual C and Si abundances differ roughly by a factor of four from the solar values, the oxygen abundance is reproduced correctly. Since this is the most common element, we expect this caveat to not qualitatively affect our conclusions. We note that other elements such as Mg, S, Ca, and Fe are produced as well, but are negligible in terms of providing additional cooling (Santoro & Shull 2006).

A key difficulty in modeling chemical mixing in SPH simulations is that no inherent mass flux between particles exists. Most simulations simply work around this caveat (e.g., Scannapieco et al. 2005; Kobayashi et al. 2007), while others implement ad hoc methods to smooth the metals within the kernel (e.g., Wiersma et al. 2009). A more realistic treatment clearly requires a physical model for the distribution of metals between particles. We have recently developed a diffusion-based algorithm to accomplish this task (Greif et al. 2009a). On all scales of interest, chemical mixing proceeds via turbulent motions that cascade down to smaller and smaller scales. By quantifying the degree of turbulence on the smallest resolved scale, one can estimate the local degree of mixing. The exchange of metals between particles is then modeled as a diffusion process, with the diffusion coefficient set by the smoothing length and the velocity dispersion within the kernel. The SPH equations for diffusion are well known (see Monaghan 2005), and only minor modifications are required to apply this machinery to chemical mixing. We note that a similar method has recently been used in the context of IGM enrichment by massive galaxies (Shen et al. 2009).

As a consequence of a more accurate metallicity estimate, we can model the additional cooling provided by heavy elements. For this purpose we assume that C, O, and Si are produced with solar relative abundances. This is a relatively good approximation for the cooling provided by the heavy elements ejected by a 200 M_☉ PISN. For these species, our model distinguishes between two temperature regimes. At temperatures below 2 × 10⁴ K, we use a slightly modified version of the chemical network presented in Glover & Jappsen (2007). This treatment follows the chemistry of C, C⁺, O, O⁺, Si, Si⁺, and Si⁺⁺, in addition to the primordial species mentioned above, and accurately models the effects of fine structure cooling from C, C⁺, O, Si, and Si⁺. It does not include the effects of cooling from molecules such as CO or H₂O, which is unimportant at the densities probed in this study (Glover & Jappsen 2007). A modification to the Glover & Jappsen (2007) treatment is the use of more accurate rate coefficients for the collisional excitation of the fine-structure lines of C and O by atomic hydrogen, taken from the recent calculations of Abrahamsson et al. (2007). The second, more important modification is that we keep C and Si in their first ionized states, since they are easily ionized by a soft UV background, which is likely to be present at z ≲ 20 (e.g., Bromm & Loeb 2003).

At temperatures above 2 × 10⁴ K, a full non-equilibrium treatment of metal chemistry, such as we use in cold gas, becomes significantly more costly, owing to the increasing number of ionization states that become available. We therefore adopt a cooling rate derived from the values given in Sutherland & Dopita (1993) for gas in collisional ionization equilibrium. Sutherland & Dopita (1993) give total cooling rates (including hydrogen and helium line cooling, and bremsstrahlung) for metallicities Z = 0, 10⁻³, 10^−2.0, 10^−1.5, 10^−1.0, 10^−0.5, 1.0, and 10^0.5 Z_☉. By subtracting off the values for the zero metallicity case from the other cases, we can isolate the contribution of the metals. For metallicities in the range 10⁻³ Z_☉ < Z < 10^0.5 Z_☉, we can then determine the metal cooling rate by interpolation of these values. For metallicities Z < 10⁻³ Z_☉ and Z > 10^0.5 Z_☉, we must extrapolate from these values, which is inevitably more inaccurate, although we note that in the former case metal cooling is ineffective compared to H and He cooling.

In view of the significant uncertainties that exist regarding the nature and the quantity of the dust produced by Pop III SNe (Salvaterra et al. 2004; Bianchi & Schneider 2007; Nozawa et al. 2007; Cherchneff & Dwek 2009), we do not include dust in our present simulation. We do not expect this to lead to a significant error: although dust cooling may be extremely important during the late stages of protostellar collapse in low-metallicity gas (see, e.g., Omukai et al. 2005; Clark et al. 2008), it is not an effective coolant at the gas densities resolved in this study.

The thermal and chemical evolution of the IGM is complicated by the fact that minihalos in the vicinity of the SN remnant collapse and form stars, which photoheat the gas that is eventually incorporated into the galaxy. Additional feedback by accretion from stellar remnant black holes or subsequent SN explosions may become important as well. However, in the present work, we choose to not include more than one SN, since we believe it is essential to understand the effects of a single SN before attempting to follow the evolution of multiple, interacting SNe. We therefore assume that all stars formed in minihalos collapse directly to black holes after they die. Recent investigations have shown that accretion onto these stellar remnant black holes is inefficient as a result of the low density of the gas after photoheating (Johnson & Bromm 2007; Alvarez et al. 2009), at least for regions of the universe that are not part of an unusually large overdensity (e.g., Gao et al. 2005). This leads us to neglect radiative feedback from early miniquasars (e.g., Kuhlen & Madau 2005).

In light of the complexity induced by multiple star formation sites, we have chosen to run two distinct simulations: in simulation A (Sim A), we insert a 50 M_☉ Pop III star whenever the density in a minihalo exceeds n_H = 100 cm⁻³, and subsequently solve for the H ii region according to the prescription in Section 2.2. After to the order of 50 stars have formed, we stop the ray-tracing routine and instead form sink particles. This ensures that the simulation does not become exceedingly complex. Once a sink particle is created, all gas particles within the smoothing length corresponding to the above density threshold are immediately accreted. Further accretion is governed by the criterion that a particle must fall within this smoothing length (details on the implementation may be found in Jappsen et al. 2005). In simulation B (Sim B), we always form sink particles and do not include any radiative feedback beyond that of the very first star. This distinction allows us to investigate the effects of photoheating on the distribution of metals and the assembly of the galaxy.

The expansion of the SN remnant into the IGM is depicted in Figures 1 and 2, where we show the central 100 kpc (comoving) of Sim A and B. After ≃15 Myr, the SN remnant has propagated to ≃1 kpc in radius, roughly the size of the original H ii region. Both simulations are identical at this point, since no other minihalos have collapsed yet. In a few cases, the shock ahead of the SN remnant has disrupted neighboring minihalos that are just in the process of collapsing. However, most minihalos remain unaffected since they are either too far away or have collapsed to high enough densities (Cen & Riquelme 2008). The distinction between disrupted and largely unaffected minihalos leads to an interesting result that will be discussed in Section 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

After roughly 100 Myr, or a Hubble time at z ≃ 30, the SN remnant comes into pressure equilibrium with the surrounding medium and stalls. During this period, multiple neighboring minihalos undergo runaway collapse and form stars. In Sim A, the ensuing H ii regions then partially overlap with the SN remnant and the comoving volume of the nascent galaxy. At this point, the complexity of the simulation dramatically increases and a multi-phase medium forms. In Figure 3, we show the distribution of the gas in density and temperature space within the high-resolution region at z ≃ 10, excluding gas accreted onto sink particles. It occupies more than 6 orders of magnitude in density, ranging from n_H = 4 × 10⁻⁵ to 100 cm⁻³, and more than 3 orders of magnitude in temperature, ranging from T = 10 to 4 × 10⁴ K. Most of the gas is shock-heated to the virial temperature and cools once the density exceeds n_H ∼ 1 cm⁻³, but some parcels of gas in Sim A are photoheated to 10⁴ K and have cooled down to T ≃ 1000 K by z ≃ 10. This warm, diffuse medium resides in the low-density IGM at n_H ∼ 10⁻³ to 10⁻⁴ cm⁻³ and will not be available for further star formation for ≳10⁸ yr, the approximate IGM free-fall time. A third phase, consisting of gas shock-heated by the SN remnant, has cooled down from an initial temperature of ∼10⁸ K to ∼10⁴ K and has become indistinguishable from the relic H ii region gas.

Zoom In Zoom Out Reset image size

Concomitant to the onset of star formation in other minihalos, the potential well of the nascent galaxy assembles at the center of a massive filament (see bottom panels of Figures 1 and 2). Accreted halos are subject to tidal stripping and rapidly lose their mass, similar to the disruption of satellite galaxies entering the halo of the Milky Way. In some cases, they complete a few orbits before being disrupted, or undergo a major merger. Interestingly, virial shocks around halos in Sim B are hotter and more distinct than in Sim A. The additional photoheating in Sim A increases the pressure at larger radii, reducing the infall velocity and, in consequence, the post-shock temperature of the gas. For the same reason, structures in Sim B are generally more pronounced, which is particularly true at the center of the galaxy (see inlays).

During the first few million years after the SN explosion, the metals are distributed into the IGM by the bulk motion of the blast wave, which leads to a tight correlation between temperature and metallicity. This becomes evident from the top panels shown in Figures 1 and 2. Once the SN remnant stalls, mixing is facilitated by agents that act on various scales. On the largest scales, filamentary accretion dominates (e.g., Wise & Abel 2007; Greif et al. 2008), while on smaller scales photoheating by neighboring stars as well as dynamical interactions associated with the virialization of the galaxy become important. As can be seen in the middle and bottom panels of Figures 1 and 2, these all interact to create a complex hierarchy of turbulent motions that efficiently mix the gas. In Figure 4, we show the resulting density–metallicity relation for the three output times in Figures 1 and 2. At early times, a distinct correlation arises: hot, underdense regions of the IGM tend to be enriched to Z ∼ 10⁻³ Z_☉ or higher, while the densest regions, such as nearby minihalos, remain almost devoid of metals. Once the potential well of the galaxy assembles, metal-rich gas accumulates at high densities and leads to a flattening of the relation.

Zoom In Zoom Out Reset image size

A more detailed view of the metal distribution is shown in Figures 5 and 6, where we plot the enriched mass as a function of gas overdensity and metallicity. Early on, the initial ejecta have not propagated very far and the metals occupy a small region in density and metallicity space, located at the relic H ii region density of ∼0.1 cm⁻³ and at supersolar metallicities. Over time, the SN remnant then distributes its metals to the surrounding medium and the average metallicity decreases, while the total enriched mass increases. Once the potential well of the galaxy has become deep enough, a large fraction of the hot, metal-rich gas residing in the IGM has recollapsed to high densities. By z ≃ 10, the enriched mass contained within the densest parcels of gas in the galaxy has increased to roughly 10⁵ M_☉, with a clear peak in the distribution at Z ∼ 10⁻³ Z_☉. Interestingly, this exceeds any critical metallicity quoted in the literature (Bromm & Loeb 2003; Schneider et al. 2006), and we expect that a stellar cluster with a normal IMF will form (e.g., Omukai et al. 2005; Schneider et al. 2006; Clark et al. 2008, 2009). It therefore appears that a single PISN is sufficient to trigger a transition from Pop III to Pop II star formation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, we note that the additional photoheating in Sim A has an interesting effect on the distribution of metals. As is evident from Figures 1 and 2, the enriched region is generally larger and better mixed in Sim A. In some cases, the photoheating even ejects significant amounts of enriched gas out of the potential well of the galaxy, which can be seen in the case of the metal bubble extending to the lower left of the simulation box. This behavior is intriguing in light of claims that the IGM is substantially enriched at high redshifts (Songaila & Cowie 1996; Schaye et al. 2003), although this is generally attributed to SN feedback instead of photoheating, which acts on objects with small circular velocities (Thoul & Weinberg 1996; Dijkstra et al. 2004). Comparing Figures 5 and 6, we also find more quantitative evidence of this effect. From top left to bottom right, an increasing amount of enriched gas accumulates at low densities, which then accretes onto the nascent galaxy.

In both simulations, we have included fine-structure cooling for C, O, and Si at low temperatures and collisional excitation and recombination cooling at high temperatures. Depending on metallicity, the latter can dominate over H and He cooling at very early times. Such a phase exists during the first few million years after the SN explosion, when the temperature of the SN remnant exceeds 10⁵ K and the metallicity within the remnant is well above solar. By performing test runs with and without metal line cooling in a smaller simulation box, we have found that metals temporarily enhance the net cooling rate by a factor of a few. However, even for the case of a PISN, the influence on the dynamical evolution of the SN remnant remains small.

At temperatures below 10⁴ K, fine-structure cooling provided by heavy elements takes over. However, we find that once the gas has cooled to the regime where fine-structure cooling becomes important, the metallicity has dropped such that cooling is instead dominated by H₂. This was already found by Jappsen et al. (2007, 2009a), who argued that dust is responsible for changing the cooling and fragmentation properties of the gas. Since cooling via dust grains kicks in at densities much higher than we can resolve here, a detailed treatment is beyond the scope of this work, although we expect that it will play a crucial role for the further evolution of the gas (e.g., Omukai et al. 2005). We note that the relative importance of the different cooling channels could be altered in the presence of a strong Lyman–Werner (LW) radiation background, which may lead to the dissociation of H₂ (e.g., Haiman et al. 1997; Johnson et al. 2007; Ahn et al. 2009).

Concomitant to the assembly of the galaxy, nearby minihalos undergo runaway collapse and form stars. Depending on distance and central density, their collapse could be unaffected, delayed, or completely inhibited in the presence of an ionization front or the SN remnant. Recent detailed radiation-hydrodynamics simulations performed by Whalen et al. (2008a, 2008b, 2010) have revealed that for a large region in parameter space most minihalos remain unaffected. However, in distinct cases the ionization front or the SN shock can disrupt the core of the halo, heating and ionizing the gas and possibly leading to metal enrichment. This in turn may trigger the formation of less massive so-called Pop III.2 stars (O'Shea et al. 2005; Johnson & Bromm 2006; Yoshida et al. 2007b), or normal Pop II stars if enough metals are brought to the center of the halo. The amount of fragmentation sensitively depends on the initial conditions of the collapse (Jappsen et al. 2009b).

How does this compare with our results? In Figure 7, we show all star-forming minihalos in the high-resolution region as a function of redshift and distance from the center of the galaxy. Crosses denote minihalos that remain largely unaffected by feedback effects and will likely form Pop III.1 stars. Minihalos that have been severely disturbed by ionizing radiation from a nearby star show elevated electron and H₂ fractions and are denoted by star symbols. These halos also form significant amounts of HD and will most likely cool to the temperature of the CMB before becoming Jeans-unstable, following the canonical pathway of Pop III.2 stars. Finally, triangles denote minihalos that have been disrupted by the SN remnant and are enriched to Z > 10⁻⁶ Z_☉. Two of these halos are even enriched to Z > 10^−3.5 Z_☉, and will almost certainly form Pop II stars.

Zoom In Zoom Out Reset image size

In Figure 8, we show a particularly interesting case that elucidates the importance of timing in regulating the strength of feedback. Here, a minihalo has been partially photo-heated by a nearby star to 10⁴ K, and subsequently enriched with metals. In the process of its collapse, the core has fragmented into two distinct objects: one in which the molecule fraction has been significantly enhanced due to the photoionization, allowing the gas to cool to ≃100 K, and another in which the cooling is unaltered from the standard minihalo pathway. In both cases, the metallicity is not high enough for metal fine-structure cooling to become important, although at higher densities dust-induced cooling may take over and lead to the formation of a small cluster of Pop II stars. This case suggests that the formation of low-mass stars need not await the formation of second-generation halos with M_vir ≳ 10⁸ M_☉.

Zoom In Zoom Out Reset image size

What are the properties of the gas within the newly virialized galaxy? From Figures 5 and 6, it appears that ∼10⁵ M_☉ of cold, dense gas at its center have been enriched to Z ∼ 10⁻³ Z_☉. This is confirmed by Figure 9, where we show its distribution in density and temperature space within the virial radius of the galaxy. It has become highly enriched by accretion from the surrounding IGM and is in a state of collapse. We also find an interesting effect related to photoheating: previously ionized gas forms elevated H₂ and HD abundances and coexists with unheated gas. Similar to the minihalo case discussed in Figure 8, the gas might therefore later fragment and form individual clumps. We note that in the density range resolved here, cooling is dominated by primordial molecules instead of metal fine-structure cooling (Jappsen et al. 2007, 2009a). At somewhat higher densities, we expect that the gas will cool to the temperature of the CMB before the equation of state hardens again. However, a more detailed investigation of the subsequent fragmentation must await high-resolution simulations that resimulate the central few kpc of the galaxy.

Zoom In Zoom Out Reset image size