POPULATION III STARS AND REMNANTS IN HIGH-REDSHIFT GALAXIES

The first generation stars (Population III) form from metal-free gas in dark matter halos with M ∼ 10⁶ M_☉ and have a large characteristic mass (e.g., Abel et al. 2002; Bromm et al. 2002; O'Shea & Norman 2007; Turk et al. 2009; Greif et al. 2012). Due to their high mass, they have short lifetimes (Schaerer 2002) and may go supernova, and enrich their surrounding intergalactic medium (IGM). For these metal-free stars, Type II supernovae (SNe) happen for an initial mass between 10 and 40 M_☉, and much more energetic pair-instability SNe (PISNe) might occur in stars between 140 and 260 M_☉ (Heger et al. 2003). Once the metallicity passes some critical metallicity, ∼10⁻⁶ Z_☉ if dust cooling is efficient (Omukai et al. 2005; Schneider et al. 2006; Clark et al. 2008) or ∼10^−3.5 Z_☉ otherwise (Bromm et al. 2001; Smith et al. 2009), the gas can cool rapidly and lower its Jeans mass. These metal-enriched (Population II) stars have a lower characteristic mass scale and most likely have an initial mass function (IMF) that resembles the present-day one. However at high redshift, heating from the cosmic microwave background may limit the radiative cooling, thus increasing the Jeans mass, resulting in an IMF that also favors massive star formation (Larson 2005; Smith et al. 2009).

The transition from Pop III to Pop II star formation is solely dependent on the metal enrichment from the Pop III SN remnants in the future star-forming halos. Metal enrichment involves complex interactions between SNe blast waves, the IGM, halo mergers, and cosmological accretion. This topic has been extensively studied with semi-analytic models (Scannapieco et al. 2003; Yoshida et al. 2004; Tumlinson 2006; Salvadori et al. 2007; Komiya et al. 2010), post-processing of numerical simulations (Karlsson et al. 2008; Trenti et al. 2009), and direct numerical simulations (Tornatore et al. 2007; Ricotti et al. 2008; Maio et al. 2010; Wise et al. 2012b; Muratov et al. 2012). For example, Trenti et al. (2009) suggested that Pop III stars may still form at the late epoch of z = 6 in the underdense regions of the universe by post-processing of cosmological simulations with blast wave models. Muratov et al. (2012) also showed that Pop III stars continue to form until z = 6 using direct cosmological simulations.

In addition to the metal enrichment from Pop III stars, heating and ionizing effects from their radiation are crucial to modeling early structure formation of the universe (Haiman et al. 2000). Lyman–Werner (LW) photons from Pop III stars photodissociate H₂ by the Solomon process and suppress the formation of Pop III stars in low-mass halos. The higher energy UV radiation from Pop III stars then affects the subsequent structure formation through heating and ionizing the surrounding IGM (Machacek et al. 2001; Yoshida et al. 2003; Wise & Abel 2007b; O'Shea & Norman 2008). Pop III more massive than 260 M_☉ or less massive than 140 M_☉ may directly collapse to form black holes (BHs; Heger et al. 2003). Accretion onto these massive Pop III BHs is a feasible way to form z > 7 quasars (Johnson et al. 2013) and is also an important source of X-ray radiation in high redshifts. X-rays from the accretion onto Pop III BHs (Alvarez et al. 2009) or from Pop III binaries (Turk et al. 2009; Stacy et al. 2010; Stacy & Bromm 2013) may preheat and pre-ionize a large volume of the IGM (Ricotti & Ostriker 2004; Mesinger et al. 2013).

The impact of Pop III feedback on cosmic evolution is dependent on the properties of their host halos, especially their masses (Whalen et al. 2008; Muratov et al. 2012). The sizes of the host halos determine the distance the metals from Pop III SNe reach with their blast waves, and the escape fraction of their UV radiation. This makes a detailed study of the Pop III star and remnant distribution over a wide range of high-redshift galaxies necessary.

It is impossible to observe Pop III stars during their lifetime at high redshift, but they might still possibly be detected directly by looking for their PISN explosions before their death. This idea has been studied and shown to be promising for both LSST (Trenti et al. 2009) and James Webb Space Telescope (Hummel et al. 2012; Whalen et al. 2013). Understanding the population and distribution of Pop III stars and remnants in high-redshift galaxies is helpful in preparing the observation of these events.

In this paper, we focus on the formation of Pop III stars and the population and multiplicity of Pop III stars and remnants in high-redshift galaxies. We have performed a simulation of a survey volume of over 100 comoving Mpc³ that includes a full primordial chemistry network, radiative cooling from metal species, both Pop II and Pop III star formation and their radiative, thermal, mechanical, and chemical feedback. The simulation runs on such a large volume from cosmological initial conditions, so that we can follow the formation and fate of Pop III stars in a statistically significant number of halos with a wide range of masses from a few million solar masses to one billion solar masses. We first describe our simulation model in Section 2. Then in Section 3, we present our results of the Pop III star and remnant distribution over early galaxies. We discuss the findings and possible bias in our simulation in Section 4.

We perform the simulation using the adaptive mesh refinement cosmological hydrodynamics code Enzo (O'Shea et al. 2004). Adaptive ray tracing (Wise & Abel 2011) is used for the radiation transfer of ionizing radiation, which is coupled to the hydrodynamics and chemistry in Enzo. The chemistry, cooling, and star formation and feedback models used in this simulation are the same as in Wise et al. (2012b).

We generate the initial conditions for the simulation using Music (Hahn & Abel 2011) at z = 99 and use the cosmological parameters from the 7 yr Wilkinson Microwave Anisotropy Probe ΛCDM+SZ+LENS best fit (Komatsu et al. 2011): Ω_M = 0.266, Ω_Λ = 0.734, Ω_b = 0.0449, h = 0.71, σ₈ = 0.81, and n = 0.963. We use a comoving simulation box of (40 Mpc)³ that has a 512³ root grid resolution and three levels of static nested grids. We first run a 512³ N-body only simulation to z = 6. Then we select the Lagrangian volume around two ∼3 × 10¹⁰ M_☉ halos at z = 6, and re-initialize the simulation with three more nested grids to have an effective resolution of 4096³ and an effective dark matter mass resolution of 2.9 × 10⁴ M_☉ inside the highest nested grid with a comoving volume of 5.2 × 7.0 × 8.3 Mpc³ (300 Mpc³). During the course of the simulation, we allow a maximum refinement level l = 12, resulting in a maximal resolution of 19 comoving pc. The refinement criteria employed are also the same as in Wise et al. (2012b). The refinements higher than the static nested grids are only allowed in the Lagrangian volume of the two massive halos at z = 6, which contains only high-resolution particles. It has a comoving volume of 3.8 × 5.4 × 6.6 Mpc³ (∼138 Mpc³) at z = 15. This highly refined volume is also the survey volume of this study. At this time, the simulation has a large number (∼1000) of halos with M > 10⁸ M_☉, where new formation of Pop III stars declines rapidly, for statistical analysis. We use results at this redshift for our current study. The simulation has 1.3 billion computational cells and required more than 10 million CPU hours, and there are three >10⁹ M_☉ halos in the refined regions at this time. We will continue this simulation to lower redshift for the study of the Pop II stars and high-redshift galaxies.

Both Pop II and Pop III stars form in the simulation, and we distinguish them by the total metallicity of the densest star-forming cell. Pop III stars are formed if [Z/H] > −4, and Pop II stars are formed otherwise. We use the same star formation models and most of the parameters in Wise et al. (2012b), as well as feedback models. For the initial mass of Pop III stars, we randomly sample from an IMF with a functional form

$$\begin{equation} f(\log M)dM=M^{-1.3}\exp \left[-\left(\frac{M_{\rm char}}{M}\right)^{1.6}\right]\, dM, \end{equation} \tag{ 1 }$$

which behaves as a Salpeter IMF above the characteristic mass, M_char, but is exponentially cut off below that mass (Chabrier 2003; Clark et al. 2009). The only difference between the two simulations is that we here use a characteristic mass of 40 M_☉ for the Pop III IMF, which is more in line with the latest results of Pop III formation simulations (e.g., Turk et al. 2009; Greif et al. 2012), instead of 100 M_☉. Please see Sections 2.2 and 2.3 of Wise et al. (2012b) for the details of the star formation and stellar feedback used in the simulation, respectively.

We show the evolution of the high-resolution region in Figure 1 at three redshifts, 25, 17.91, and 15. Here we show the large-scale structure in the inner 6.6 comoving Mpc in the top row, and in the remaining rows, we focus on the most massive halo at each redshift by showing their density-weighted projections of gas density, temperature, and metallicity.

Zoom In Zoom Out Reset image size

We illustrate the evolution of the number of Pop III stars and remnants, as well as their formation rate from the birth of the first Pop III star to z = 15 in the top panel of Figure 2. Massive Pop III stars have very short lifetimes (Schaerer 2002) and end their life by either directly collapsing to BHs or exploding as SNe (Heger et al. 2003), depending on their initial masses. More specifically, they die as Type II SNe if 11 ⩽ M_⋆/M_☉⩽ 40 and as PISNe if 140 ⩽ M_⋆/M_☉ ⩽ 260, where M_⋆ is the initial stellar mass, or become BHs if their masses are not in these mass ranges. Throughout the paper, we use "Pop III remnants" to refer to all remains of dead Pop III stars, regardless of whether they become BHs or supernova remnants, and we use "Pop III" to refer to both living stars and dead remnants. In the case that remnants have negligible masses after SN events, their star masses are replaced with a very small mass proportional to their initial masses and the star particles are kept in the simulation to follow their remnant kinematic distribution. The first Pop III star forms at redshift z ∼ 29.7. Rates of Pop III star formation steadily increase from ∼10⁻⁸ to above 10⁻⁶ stars per year per comoving Mpc³. The Pop III star formation rate (SFR) shows some signs of saturation at z ∼ 15 in this overdense region of the universe. At z = 15, the entire survey volume of 138 comoving Mpc³ contains 13,123 Pop III stars and remnants and 7677 halos more massive than 5 × 10⁶ M_☉. The number densities of Pop III and halos that host Pop III are 95 and 55 per comoving Mpc³, respectively.

Zoom In Zoom Out Reset image size

The total mass of Pop II stars and the SFR of Pop II and Pop III are shown in the bottom panel of Figure 2. There is ∼3 × 10⁸ M_☉ mass in Pop II stars at z = 15. The star formation histories for both Pop II and Pop III stars are similar to those in Wise et al. (2012b), but shifted to higher redshifts. We will study the details of Pop II star formation in this simulation in a forthcoming paper.

We show the number of Pop III, all halos, and star-hosting halos as functions of halo mass inside the survey volume at z = 17.91 and z = 15 in the left panels in Figure 3, while the numbers of living Pop III stars are stacked over remnants in the same bin. In the right panels in Figure 3 are shown the fraction of halos hosting Pop III and the average Pop III per halo (over halos having Pop III) as functions of halo mass.

Zoom In Zoom Out Reset image size

As our survey volume is chosen over a high-density region, the halo number density is higher than the mean in the universe. The halo mass functions from our simulation are well fit with the fitting function of Warren et al. (2006) boosted by a factor of five for halos heavier than 5 × 10⁶ M_☉ in both redshifts. The number of halos drops for halos <5 × 10⁶ M_☉, due to the lack of dark matter mass resolution. At z = 15, our simulation has a factor of a few more halos than the prediction from the fitting function for halos more massive than 1 × 10⁸ M_☉, showing that the baryon becomes important in massive halos. There are three halos with masses over 10⁹ M_☉ at z = 15, while the first one appears at redshift z ∼ 15.8.

The Pop III distribution over halo mass shows little evolution, except that there are more Pop III appearing in high-mass halos as those halos build up with time. The number of Pop III peaks at ∼3 × 10⁷ M_☉. The fraction of halos hosting Pop III also has little change between redshifts 17.91 and 15 because the chances to host Pop III are very small for halos less massive than 10⁷ M_☉ in this biased region. Then they gradually increase with halo mass and reach 100% for halos ∼10⁸ M_☉. In this particular mass range of 10⁷–10⁸ M_☉, halos can cool and collapse regardless of the LW background and form stars (e.g., Wise & Abel 2007b; O'Shea & Norman 2008; Shang et al. 2010; Wolcott-Green et al. 2011). The exact timing of Pop III star formation depends on the local strength of the LW intensity, causing the gradual increase in the fraction shown in Figure 3. In addition to H₂ suppression, nearby galaxies and Pop III supernovae can chemically enrich halos so that they form metal-enriched stars in their first star formation event without hosting any prior Pop III star formation. Afterward, as halos grow through hierarchical mergers, the number of Pop III per halo increases with both redshift and halo mass. At z = 15, halos between 10⁷ M_☉ and 10⁸M_☉ have an average of 10 Pop III, and the number rises to about 50 for halos ∼10⁹ M_☉.

Shown in Figure 4 are the number distribution of halos as a function of halo mass and number of Pop III stars and remnants per halo. The number of Pop III that is hosted by halos inside a mass bin is very scattered, but there is a clear trend that more massive halos host more Pop III stars and remnants. The same trend is shown in the average number of Pop III per halo in Figure 3. At redshift z = 15, the most massive halo without any Pop III stars or remnants is 7.16 × 10⁷ M_☉. At the low-mass end, the least massive halo hosting Pop III is only 1.82 × 10⁶ M_☉. But this Pop III star forms at about 50 Myr ago, so this star most likely forms in another more massive halo, then it is stripped out with its current host halo from the massive one.

Zoom In Zoom Out Reset image size

In Figure 5, we show the number distribution of Pop III as a function of their host halo mass and their formation time. As old Pop III remnants can be found in halos of all masses, young stars and remnants can only be seen in low-mass halos. No Pop III star forms in halos more massive than 3 × 10⁸ M_☉. The gas in massive halos is enriched by their earlier Pop III, so the further formation of stars transitions to Pop II. The least massive halo having newly formed Pop III stars is ∼4 × 10⁶ M_☉, as the Pop III formation in smaller halos is suppressed by the LW background. But our results for Pop III formation in low-mass halos might be biased. Due to the relatively low dark matter (2.9 × 10⁴ M_☉), we may miss the formation of M ≲ 10⁶ M_☉ halos in our simulation. We discuss this issue in the next section and show that very little Pop III star formation is missed by the current simulation due to resolution effects.

Zoom In Zoom Out Reset image size

We present results of a cosmological simulation representing a high-density region of an unprecedentedly large volume in the early universe to study the Pop III star and remnant population and multiplicity within high-redshift galaxies. The highly resolved simulation volume covers thousand of halos and Pop III stars and remnants, allowing us to make a detailed study of the Pop III star formation history and the Pop III star and remnant distribution over a wide range of halo masses. We observe a continuous Pop III star formation from z ∼ 30 to z ∼ 15, when the calculation was paused for analysis. We expect Pop III star formation to continue to lower redshifts, but this needs to be verified by running the simulation further. We find that the number of Pop III stars and remnants peaks in halos with masses of a few ×10⁷ M_☉ during this time, while the ratio to host Pop III and multiplicity of Pop III increase with the halo mass. The Pop III inside a massive halo are likely coming from mergers of small ones, instead of forming in the big ones. New Pop III stars only form in halos between 4 × 10⁶ and 3 × 10⁸ M_☉, while the formation rate peaks at ∼10⁷ M_☉. So even though there are tens of Pop III remnants in a halo of 10⁹ M_☉, little to no feedback might be expected coming from their still living Pop III stars.

The dark matter mass resolution in this simulation is not high enough to capture all of the Pop III star-forming halos smaller than a few ×10⁶ M_☉, so we could be underestimating the number of Pop III in our simulation. To investigate this, we use the simulation in Wise et al. (2012a) to estimate the missing Pop III stars and remnants in low-mass halos. We show the histogram of halo mass in that 1 comoving Mpc³ box simulation at z = 7.28 in Figure 6. Due to the much higher dark matter mass resolution, there are many more halos at a few 10⁶ M_☉, extending to ∼2 × 10⁵M_☉. Since the simulation box is small, the halo mass function is not well fit by the Warren et al. (2006) function by having too many low-mass halos and too few high-mass halos, but still has the number of halos ∼10⁶ M_☉ close to the cosmic mean. Though there are far fewer halos and Pop III in this smaller simulated volume, it has a similar distribution of Pop III stars and remnants as in our current simulation. Halos with M_halo ≲ 2 × 10⁶ M_☉ do not host any Pop III stars or remnants in this small box simulation. Pop III star formation in these low-mass halos is suppressed when the H₂ is dissociated by LW radiation from nearby stellar sources and the background (Machacek et al. 2001; Wise & Abel 2007b; O'Shea & Norman 2008). The first star forms at z ∼ 20, when the LW background is already at an intensity of ∼10⁻²³ erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹, and suppresses any H₂ formation and thus Pop III star formation in halos below a mass of 10⁶ M_☉. At later times, we also confirm that the minimum star-forming halo mass with this LW background in Wise et al. (2012a) at any time is ≳ 10⁶ M_☉. Thus, we conclude that the Pop III formation rate is small in halos below this mass scale and that we do not miss a significant fraction of Pop III star formation in our simulation because of mass resolution.

Zoom In Zoom Out Reset image size

In this large volume simulation, only the LW flux from point sources is included, while the LW background is ignored. We show the evolution of averaged LW intensity and the fraction distribution of LW intensity at z = 15 over the refined volume in Figure 7. The averaged LW radiation intensity builds up gradually and reaches 10⁻²³ erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹ at redshift z ∼ 17, and maintains at this level for more than 40 Myr. At z = 15, the LW intensity in ∼70% of the volume is higher than 10⁻²³ erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹, while only less than 0.01% of the volume (0.014 comoving Mpc³) is filled with LW radiation stronger than 10⁻²¹ erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹. This LW intensity, locally produced even in a rare high-density region, is still weaker than the expected LW background (Wise & Abel 2005; Wise et al. 2012a). Usually, it is thought that this level of LW radiation may only suppress Pop III formation in ∼10⁶ M_☉ halos but is not strong enough to have a significant impact in larger minihalos (3 × 10⁶–10⁷ M_☉). However, we found that only 20% of halos in this mass range form Pop III stars, delayed by LW feedback. At z = 15, 90% (respectively two-thirds) of the halos with no Pop III stars and M = 5 × 10⁶–3 × 10⁷ M_☉ have LW intensities higher than 10⁻²³ (resp. 2 × 10⁻²³) erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹. The main difference in our simulation to previous studies of LW feedback is environment. In our simulated overdense region, the halo mass accretion rates are high enough so that dynamical heating (Yoshida et al. 2003) becomes important, where mergers and rapid cosmological inflows cause shock heating and virial turbulence (Wise & Abel 2007a; Greif et al. 2008) that can disrupt any dense central cores. Embedding these rapidly growing minihalos in an LW background delays their collapse even further as these disrupted cores must cool once again and re-coalesce. The typical collapse time to form Pop III stars in a 10⁻²³ erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹ LW background is 50 Myr (O'Shea & Norman 2008). Most of the minihalos in our simulation at these very high redshifts are younger than 50 Myr, and thus have not had sufficient time to collapse.

Zoom In Zoom Out Reset image size

An important difference between the simulation in this paper and that of Wise et al. (2012a) is that the characteristic mass of Pop III stars here is much lower at M_char = 40 M_☉, instead of M_char = 100 M_☉. The new choice of characteristic mass results in far fewer Pop III stars having initial mass between 140 and 260 M_☉, and ending their life by PISN. This then significantly reduces the metal generated by Pop III stars. For example, a 40 M_☉ hypernova produces only 8.6 M_☉ metals (Nomoto et al. 2006), compared to 85 M_☉ of metals generated by 180 M_☉ PISN (Heger & Woosley 2002). It is likely due to the lower level of metal generation that Pop III stars continue to form in the larger halos (>10⁸M_☉) in the current simulation. On the contrary, there is no Pop III formed in halos more massive than 5 × 10⁷ M_☉ in Wise et al. (2012b).

Because our simulation covers a rare high-density region, our results are not representative of the cosmic mean. We fit the halo mass function with the fitting function of Warren et al. (2006). It shows that the halo number density at z = 15 is about five times that of the cosmic mean at the same redshift or about the cosmic mean at z = 10, while the halo number density at z = 17.91 is about the cosmic mean at z = 14. The density fluctuations may have effects on local Pop III formation by changing the LW background and metallicity level of the IGM. As discussed before, the strength of LW background during these redshifts might not change the Pop III formation significantly. In addition, because the mean distance between Pop III and halos is still much larger than the direct metal impact distance of a PISN of 5 kpc (Wise et al. 2012b), the metal enrichment by Pop III feedback crossing halos is still weak, except in halo clusters. In this case, the Pop III formation mostly depends on single halo properties, while the multiplicity of Pop III over halos is determined by the clustering of halos. Our results of Pop III distribution and multiplicity over halos should not change significantly in regions of different density fluctuations. The population of Pop III is proportional to the halo density, so the number density from this simulation at z = 15 should be about five times that of the cosmic mean or close to the cosmic mean at z ∼ 10.

One important question about Pop III stars is when and in what environments their formation completely ends. Unfortunately, our simulation does not show the cessation, or even the saturation of Pop III formation yet, though the increase of SFR slows down in the last tens of millions of years. The Pop III star and remnant number density is still low (lower than that of Wise et al. (2012b) 1 Mpc³ simulation), and there still is a significant fraction of halos between 10⁷ and 10⁸ M_☉ that have no Pop III stars and are in a favorable LW radiation field. To show the sufficiency of metal-poor gas inside the survey volume to continue to form Pop III stars, we plot the volume and mass complementary cumulative distributions of the metallicity at z = 15 in Figure 8. There are only about 0.65% of volume and 3.8% of the mass enriched above 10⁻⁶ Z_☉. For comparison, the Wise et al. (2012b) 1 Mpc³ simulation has 0.32% of the volume and 2.6% of the mass at z = 10, and 0.89% of the volume and 7.6% of the mass at z = 7.28, respectively, above 10⁻⁶ Z_☉. So it is reasonable to expect Pop III stars will continue to form inside our survey volume, but at a constant rate, since Wise et al. (2012b) show that the Pop III formation rate is constant from z = 12 to 7. We are continuing the present simulation to lower redshifts to confirm or deny this expectation.

Zoom In Zoom Out Reset image size

We thank an anonymous referee for helpful comments. This research was supported by the National Science Foundation (NSF) grant AST-1109243 to M.L.N. J.H.W. acknowledges support from the NSF grant AST-1211626. The simulation was performed on the Kraken supercomputer operated for the Extreme Science and Engineering Discovery Environment (XSEDE) by the National Institute for Computational Science, ORNL with XRAC allocation MCA-TG98020N. This research has made use of NASA's Astrophysics Data System Bibliographic Services. The majority of the analysis and plots were done with yt (Turk et al. 2011).