DETECTING THE RISE AND FALL OF THE FIRST STARS BY THEIR IMPACT ON COSMIC REIONIZATION

The theory of reionization has not yet advanced to the point of establishing unambiguously its timing and the relative contributions to it from galaxies of different masses. In a cold dark matter ("CDM") universe, these early galactic sources can be categorized by their host halo mass into minihalos ("MHs") and "atomic-cooling" halos ("ACHs"). MHs have masses M ∼ 10⁵–10⁸ M_☉ and virial temperatures T_vir ∼ 10⁴ K, and thus molecular hydrogen (H₂) was necessary to cool the gas below this virial temperature to begin star formation. ACHs have M ≳ 10⁸ M_☉ and T_vir ≳ 10⁴ K for which H-atom radiative line cooling alone was sufficient to support star formation. The ACHs can be split further into low-mass atomic-cooling halos ("LMACHs"; M ∼ 10⁸–10⁹ M_☉), for which the gas pressure of the photoionization-heated intergalactic medium ("IGM") in an ionized patch prevented the halo from capturing the gas it needed to form stars, and high-mass atomic-cooling halos ("HMACHs"; M ≳ 10⁹M_☉), for which gravity was strong enough to overcome this "Jeans-mass filter" and form stars even in the ionized patches.

Once starlight escaped from galactic halos into the IGM to reionize it, the ionized patches ("H ii regions") of the IGM became places in which star formation was suppressed in both MHs and LMACHs. At the same time, UV starlight at energies in the range of 11.2–13.6 eV also escaped from the halos, capable of destroying the H₂ molecules inside MHs through Lyman–Werner band ("LW") dissociation, even in the neutral zones of the IGM. This dissociation eventually prevented further star formation in some of the MHs where the background intensity was high enough. Early estimates, in fact, suggested that this would have made the MH contribution to reionization small (Haiman et al. 2000), and, until now, large-scale simulations of reionization have neglected them altogether.

In this Letter, we report the first radiative transfer (RT) simulations of reionization to include all three of the mass categories of reionization source halos, along with their radiative suppression, in a simulation volume large enough to capture both the global mean ionization history and the observable consequences of its evolving patchiness in a statistically meaningful way.⁶ We overcame the limitation of previous large-volume simulations by applying a newly developed sub-grid treatment to include MH sources (Section 2) and calculating the transfer of LW-band radiation self-consistently with the source population using the scheme of Ahn et al. (2009).

We performed a cosmological N-body simulation of structure formation with 3072³ particles in a 114 h⁻¹ Mpc simulation box, using the WMAP5 background cosmology (Dunkley et al. 2009). For this we used the code CubeP³M (Merz et al. 2005; Iliev et al. 2008; J. Harnois-Deraps et al. 2012, in preparation) in which the gravity is computed by a particle-particle-particle-mesh (P³M) scheme. The simulation was started at redshift z = 300 and run to z = 6. N-body data were recorded at 86 equally spaced times (every 11.53 Myr) from z = 50 to z = 6. Each data time slice was then used to create matter density fields by smoothing the particle data adaptively onto a uniform mesh—or an "RT grid"—of 256³ cells. All cosmological halos with M ⩾ 10⁸ M_☉ (corresponding to 20 particles or more), and thus both LMACHs and HMACHs, were identified on the fly using a spherical overdensity halo finder with overdensity of Δ = 178 with respect to the mean.

Because MHs are too small to be resolved in our simulation box, our RT grid was populated with MHs through a newly developed sub-grid model, as follows. We started with a separate, high-resolution N-body simulation of structure formation in a box with (6.3 h⁻¹ Mpc)³ volume and 1728³ particles, which resolved all MHs with M ⩾ 10⁵M_☉ with 20 particles or more. We then partitioned the box into a uniform grid of 14³ cells, such that each cell is the same size as one of the RT grid cells in our main, 114 h⁻¹ Mpc simulation box, and calculated cell density and the total number of MHs per cell. A strong and tight correlation between the number of MHs located in a cell and its density is observed (Figure 1). The best fit to this correlation at each redshift was then used as the total number of MHs in each grid cell of 114 h⁻¹ Mpc box.

Zoom In Zoom Out Reset image size

Based on these structure formation results for the IGM density field and the source dark matter halos, we then calculated the RT of H-ionizing and H₂-dissociating photons (see Table 1 for the RT simulation parameters). The stars inside ACHs are assumed to produce g_γ ionizing photons per baryon every 10 Myr, where g_γ ≡ f_γ/(t_⋆/10 Myr), and where f_γ ≡ f_ef_⋆N_i, f_e is the escape fraction of ionizing photons, f_⋆ is the star formation efficiency, and N_i is the number of ionizing photons per stellar baryon produced over the star's lifetime t_⋆—we use t_⋆ = 11.53 Myr for both HMACHs and LMACHs, and t_⋆ = 1.92 Myr for MHs. We assign one Pop III star per MH, motivated by numerical simulations of first star formation inside MHs, which find that typically one Pop III star with a mass between 100 and 1000 M_☉ forms per MH in the absence of strong soft UV radiative feedback (Bromm et al. 2002; Abel et al. 2002; Yoshida et al. 2006). At each cell, only those MHs which are newly collapsed every 1.92 Myr are assumed to host Pop III stars to roughly approximate the disruptive radiative and mechanical feedback by the first star and its by-products (such as a supernova) on the halo gas (for this, we create "morphed" density fields every 1.92 Myr by linearly interpolating the N-body density fields in time which are separated by a time interval of 11.53 Myr and finite difference corresponding MH populations on each cell). Star formation in MHs is further suppressed when the local LW background—calculated at each time step in three dimensions using the scheme by Ahn et al. (2009), but now considering both ACHs and MHs and also improved in speed using the fast Fourier transform (FFT) scheme—reaches a certain threshold J_LW, th. At present, the precise value of this threshold is not well determined, but the typical values found by high-resolution simulations of MH star formation are J_LW, th = [0.01–0.1] × 10⁻²¹ erg s⁻¹ cm⁻² sr⁻¹ (Machacek et al. 2001; Yoshida et al. 2003; O'Shea & Norman 2008). We adopted a constant value chosen from this range for each simulation. Even though stars may still form by H₂ cooling, when J_LW > J_LW, th, in MHs in the mass range M ≃ 2 × 10⁶–10⁸ M_☉ and M ≃ 10⁷–10⁸ M_☉ at J_LW, th ≃ 0.01 × 10⁻²¹ erg s⁻¹ cm⁻² sr⁻¹ and J_LW, th ≃ 0.1 × 10⁻²¹ erg s⁻¹ cm⁻² sr⁻¹, respectively (e.g., O'Shea & Norman 2008), we neglect their contribution because they constitute only a small fraction of the whole MH population.

The simulations with ACHs only (and without LW RT) are described in Iliev et al. (2012). Simulation parameters are given in Table 1.

We demonstrate the effects of the first stars by direct comparison of the results from two simulations: a fiducial case which includes all ionizing sources down to the first stars hosted by MHs (case L2M1J1) and a corresponding reference case which includes the larger, atomically cooling halos with exactly the same properties, but no MH sources (case L2, previously presented in Iliev et al. 2012). Our results show that the early reionization history is completely dominated by the first stars, while the late (redshift z ≲ 10) history is driven by the stars inside HMACHs (Figures 2(a) and (b), top panel). The very first stars start to form inside MHs at redshift z ≃ 40 and dominate the reionization process until z ≃ 10 but through self-regulation which slows their contribution to reionization⁷ (Figure 2(b)). Although the abundance of ACHs rises exponentially, they remain relatively rare, and thus sub-dominant, until z ≃ 10. After redshift z ≃ 8, though, the two reionization histories become largely indistinguishable because the same HMACHs then dominate reionization and push J_LW above J_LW, th (at z ≃ 12), halting MH star formation altogether, long before the MHs can complete reionization on their own.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Nonetheless, the MH sources (the first stars) can have quite a dramatic effect on the electron-scattering optical depth τ_es. While intergalactic H ii regions fully overlap (at redshift z_ov, here defined as when the mass-weighted mean ionized fraction in the IGM, x_m, first surpassed 99%) at almost identical redshifts z_ov ≃ 6.8 with (L2M1J1) or without (L2) MHs, the early rise of x_m with MH sources boosts the optical depth by as much as 47 % relative to that without MH sources: τ_es = 0.0861 for L2M1J1, while τ_es = 0.0603 for L2. This satisfies the current observational constraints on reionization: (1) reionization ended no earlier than redshift z = 7 (Fan & et al. 2006; Ota et al. 2010; Mortlock et al. 2011; Bolton et al. 2011; Pentericci et al. 2011) and (2) τ_es = 0.088 ± 0.015 at 68 % confidence level (Larson et al. 2011). Predicted values of τ_es and z_ov are model dependent, and thus we tested the robustness of our conclusions by varying the physical parameters of MHs and ACHs (Figure 3).

The first stars, born inside MHs, imprint a distinctive pattern on the global reionization history. For example, in case L2M1J1, when the LW plateau ended, reionization briefly stalled, since MHs no longer formed the stars which replenished the ionizing background and only ACH sources remained, thereafter; the ACH contribution took a bit more time to climb enough to move reionization forward again. This explains the brief "x_m-plateau" from z ∼ 12 to z ∼ 10 in Figure 2(b) for case L2M1J1, while in case L2, x_m grows continuously without showing such a plateau (Figure 2(b)). This feature is generic (see Figures 3 and 4(b) for different sets of parameters we explored). Reionization histories without MH sources, modeled either by large-scale RT simulations (Ciardi et al. 2003; Iliev et al. 2006, 2007; Trac & Cen 2007; McQuinn et al. 2007; Zahn et al. 2007; Iliev et al. 2012) or semi-analytical calculations (Haiman et al. 2000; Zahn et al. 2007), are all similar in that respect.

Zoom In Zoom Out Reset image size

Reionization histories with MH sources calculated here, however, find an ionization plateau phase. Previous studies that considered MH stars and their impact were not able to settle the issue of their global effect on reionization. This is either because they simulated volumes much too small to represent a fair portion of the universe (Ricotti et al. 2002; Sokasian et al. 2004; Yoshida et al. 2006; Ricotti et al. 2008; Johnson et al. 2012) or else treated reionization by a semi-analytical, 1-zone, homogeneous approximation (either with LW suppression included (Haiman et al. 2000; Furlanetto & Loeb 2005) or without (Shapiro et al. 1994; Wyithe & Loeb 2003; Haiman & Bryan 2006; Wyithe & Cen 2007)), which cannot capture its innate spatially inhomogeneous nature, or made a semi-analytical approximation that accounted statistically for spatial inhomogeneity but without LW suppression (Kramer et al. 2006).

We find that the global reionization history at z ≲ 20 depends on J_LW, th more strongly than M_⋆, III. This is due to the very nature of self-regulation of the first star formation. The larger the J_LW, th is, the weaker the suppression of star formation becomes, thereby temporarily hastening the progress of reionization. If M_⋆, III is smaller (Turk et al. 2009; Stacy et al. 2010; Greif et al. 2011a; Stacy et al. 2012) than those simulated here, those stars produce less ionizing and LW radiation and the resulting suppression is weaker, which partly compensates for the lower emission per star. Similar type of compensation would occur also when the number of MHs with a potential to form the first stars is smaller than our estimate due to the relative offset of baryonic gas from some of the MHs (Tseliakhovich & Hirata 2010; Greif et al. 2011b).

We thus conclude that the first stars hosted by MHs likely made an important contribution to reionization. But how can we probe them observationally? While significant, the effects of the first stars are largely confined to the early stages of reionization, at redshifts z > 10, which puts them beyond the reach of most current instruments. Recent observations by the South Pole Telescope (SPT) have been used to place an upper limit on the kinetic Sunyaev–Zel'dovich (kSZ) effect from the epoch of reionization (Reichardt et al. 2012). While it has been suggested that this restricts the duration of reionization (Zahn et al. 2011), we will show elsewhere (H. Park et al. 2012, in preparation) that the kSZ signal from our fiducial case, L2M1J1, is well below the observed upper bound. The combined effect of the first stars will be reflected in the cosmic microwave background (CMB) polarization anisotropies at large scales. The current best constraints on τ_es by the Wilkinson Microwave Anisotropy Probe satellite are still relatively weak, and thus models with low τ_es values like L2 are still acceptable at the 2σ (95%) confidence level (Figure 2(b), middle panel). However, through the far more precise measurement of the CMB polarization by the Planck mission we should be able to discern the influence of the first stars on reionization.

As we discussed above, current observational constraints suggest that reionization was not complete before z ∼ 7. Imposing this condition as a prior on the allowed reionization histories x_m(z), we predict that the Planck mission will clearly detect the era of first stars (Figure 4). In Figure 4, we show a statistical measure of the Planck sensitivity in detecting the signature of the first stars through the principal component analysis by Hu & Holder (2003) and Mortonson & Hu (2008, who modified COSMOMC (Lewis & Bridle 2002) to allow generic reionization models).⁸ Reionization principal components {S_μ(z)} are eigenvectors of the relevant Fisher matrix (evaluated with an arbitrarily chosen fiducial history x_e, fid(z)),

$$\begin{equation} F_{ij}\equiv {\displaystyle \sum _{l=2}^{l_{{\rm max}}}}\left(l+\frac{1}{2}\right) \frac{\partial ^{2}C_{l}^{{\rm EE}}}{\partial x_{e}(z_{i})\partial x_{e}(z_{j})}, \end{equation} \tag{ 1 }$$

which can be used to describe any generic ionization history x_e(z) with just a small number of modes, such that

$$\begin{equation} x_{e}(z)=x_{e,\,{\rm fid}}(z)+{\displaystyle \sum _{\mu =1}^{N_{{\rm max}}}}m_{\mu }S_{\mu }(z). \end{equation} \tag{ 2 }$$

The mode amplitude m_μ, for a given history x_e(z), becomes

$$\begin{equation} m_{\mu }=\frac{\int _{z_{\rm min}}^{z_{\rm max}} dz \,S_{\mu }(z)[x_{e}(z)-x_{e,\,{\rm fid}}(z)]}{z_{\rm max}-z_{\rm min}}. \end{equation} \tag{ 3 }$$

Based on the Planck data after its full two years of planned operation, the narrow posterior distribution of allowed τ_es values will allow us to distinguish reionization models like L2 and L2M1J1 unambiguously, and thereby strongly constrain the available reionization models. A high measured value of τ_es > 0.085 will be a clear (if indirect) signature of the first stars.

Finally, we note that the presence of MH sources introduces the x_m plateau noted above, which in turn imprints characteristic features on C^EE_l. Hence, Planck might be able to distinguish (albeit at lower statistical significance, of ≳ 2σ or ≳ 95%) reionization models with and without first stars even if they have very similar values of τ_es and z_ov (Figure 4(b)). Full reionization simulations like ours find it hard to satisfy both of these observational constraints without including a significant contribution from the first stars, but some semi-analytical models (Haiman & Bryan 2006; Haardt & Madau 2012; g0.348_67.8 and g2.609C_165.2 in Figure 4(b) which are of unnaturally large gaps in relative efficiencies of LMACH and HMACH) do find such scenarios. However, all such models lack the plateau feature in x_m(z), regardless of the details of the assumed physics, and reside in a narrow window of the m_μ-parameter space adjacent to that occupied by our no-MH cases, as demonstrated in Figure 4(b).

In summary, Planck is capable of distinguishing with high confidence between definitive classes of reionization scenarios allowed by the current constraints, and thereby significantly restricting the available parameter space. Planck will either probe the signature of the first stars or show that the first stars had a negligible impact on reionization. Once these first results confirm the role of the first stars, simulations of the type presented here can be used to study other observable quantities and thus deepen our understanding of the early universe.

K.A. was supported in part by NRF grant funded by the Korean government MEST (Nos. 2009-0068141, 2009-0076868, 2012R1A1A1014646, and 2012M4A2026720). I.T.I. was supported by The Southeast Physics Network (SEPNet) and the Science and Technology Facilities Council grants ST/F002858/1 and ST/I000976/1. This study was supported in part by the Swedish Research Council grant 2009-4088; U.S. NSF grants AST-0708176 and AST-1009799; NASA grants NNX07AH09G, NNG04G177G, and NNX11AE09G; and Chandra grant SAO TM8-9009X. The authors acknowledge the TeraGrid and the Texas Advanced Computing Center (TACC) at The University of Texas at Austin (URL: http://www.tacc.utexas.edu), and the Swedish National Infrastructure for Computing (SNIC) resources at HPC2N (Umeå, Sweden) for providing HPC and visualization resources that have contributed to the results reported within this paper. We acknowledge A. Lewis, A. Liddle, and M. Mortonson for scientific and technical input on Planck forecasts (Lewis, Liddle) and modifying COSMOMC for reionization principal components (Mortonson).