CRITICAL STAR FORMATION RATES FOR REIONIZATION: FULL REIONIZATION OCCURS AT REDSHIFT z ≈ 7

We denote by $\dot{\rho }_{\rm SFR}$ (M_☉ yr⁻¹ Mpc⁻³) the global SFR per comoving volume. Using a simple argument (Madau et al. 1999), balancing photoionization with radiative recombination, we estimate the critical SFR to maintain IGM photoionization at z > 7, assuming that the LyC photons are produced by populations of massive (OB-type) stars. Because the mass in collapsed objects (clusters, groups, galaxies) is still small at high redshift, the IGM contains most of the cosmological baryons, at mean density

$$\begin{equation} \bar{ \rho }_b = \Omega _b \rho _{\rm cr} (1+z)^3 = (4.24 \times 10^{-31} {\rm \,g\, cm}^{-3}) (1+z)^3 \;. \end{equation} \tag{ 1 }$$

For a Hubble constant denoted H₀ = (100 km s⁻¹ Mpc⁻¹) h, we adopt the WMAP-7 (plus BAO + H₀) parameters, Ω_bh² = 0.02255 ± 0.00054 and Ω_mh² = 0.1352 ± 0.0036 (Komatsu et al. 2011) relative to a critical density ρ_cr = 1.8785 × 10⁻²⁹h² g cm⁻³. From the corresponding helium mass fraction Y = 0.2477 ± 0.0029 (Peimbert et al. 2007), we adopt a mean hydrogen number density

$$\begin{equation} \bar{n}_{\rm H} = \frac{ \bar{\rho }_b (1-Y)}{m_{\rm H}} = (1.905 \times 10^{-7} {\rm \,cm}^{-3}) (1+z)^3 \;. \end{equation} \tag{ 2 }$$

In a fully ionized IGM, the CMB optical depth back to z_rei can be written as the integral of n_eσ_Tdℓ, the electron density times the Thomson cross section along proper length,

$$\begin{equation} \tau _{e}(z_{\rm rei}) = \int _{0}^{z_{\rm rei}} n_{e} \sigma _{\rm T} (1+z)^{-1} \; [c/H(z)] \; dz \;, \end{equation} \tag{ 3 }$$

for a standard ΛCDM cosmology (Ω_m + Ω_Λ = 1) with H(z) = H₀[Ω_m(1 + z)³ + Ω_Λ]^1/2. This integral can be done analytically (Shull & Venkatesan 2008):

$$\begin{eqnarray} \tau _{e}(z_{\rm rei}) &=& \left(\frac{c \, \sigma _{\rm T}}{H_0} \right) \left(\frac{2 \Omega _b}{3 \Omega _m} \right) \left[ \frac{\rho _{\rm cr} (1-Y)(1+y) }{m_{\rm H}} \right]\nonumber\\ && \times \lbrace [ \Omega _m (1+z_{\rm rei})^3 + \Omega _{\Lambda } ] ^{1/2} - 1\rbrace . \end{eqnarray} \tag{ 4 }$$

In the high-redshift limit, when Ω_m(1 + z)³ ≫ Ω_Λ, this expression simplifies to

$$\begin{eqnarray} \tau _{e}(z_{\rm rei}) &\approx & \left(\frac{c \, \sigma _{\rm T}}{H_0} \right) \left(\frac{2 \Omega _b}{3 \Omega _m^{1/2}} \right) \left[ \frac{\rho _{\rm cr} (1-Y)(1+y)}{m_{\rm H}} \right]\nonumber\\ && \times (1+z_{\rm rei})^{3/2} \approx (0.0521) \left[ \frac{(1+z_{\rm rei})}{8} \right] ^{3/2} \end{eqnarray} \tag{ 5 }$$

independent of h to lowest order. The helium and electron densities are written n_He = yn_H and n_e = n_H(1 + y) for singly ionized helium, where y = n_He/n_H = (Y/4)/(1 − Y) ≈ 0.0823 by number. To these formulae, we add Δτ_e ≈ 0.002, from electrons donated by He iii reionized at z ⩽ 3 (Shull et al. 2010). Helium therefore contributes ∼8% to τ_e, and a fully ionized IGM produces τ_e = 0.050, 0.060, and 0.070 back to redshifts z_rei = 7, 8, and 9, respectively.

A comoving volume of 1 Mpc³ contains N_H = 5.6 × 10⁶⁶ hydrogen atoms. Our simple ionization criterion requires an SFR density that produces a number of LyC photons equal to N_H over a hydrogen recombination time, t_rec = [n_eα^(B)_HC_H]⁻¹. The hydrogen case-B recombination rate coefficient (Osterbrock & Ferland 2006) is α^(B)_H(T) ≈ (2.59 × 10⁻¹³ cm³ s⁻¹)T₄^−0.845, scaled to an IGM temperature T = (10⁴ K)T₄. For typical IGM ionization histories and photoelectric heating rates, numerical simulations predict that diffuse photoionized filaments of hydrogen have temperatures ranging from 5000 K to 20,000 K (Davé et al. 2001; Smith et al. 2011). These are consistent with temperatures inferred from observations of the Lyα forest at z < 5 (Becker et al. 2011).

Owing to gravitational instability, a realistic IGM is inhomogeneous and filamentary. Semi-analytical models of the reionization of the universe often adopt a "clumping factor," C_H ≡ 〈n²_e〉/〈n_e〉², to account for inhomogeneity in estimates of the enhanced recombination rate in denser IGM filaments. The clumping factor therefore plays an important role in computing the critical SFR density needed to maintain the reionization of the universe. The clumping factor is also used in numerical simulations to implement "sub-grid physics," in which changes in the density field occur on scales below the resolution of the simulation and are also approximated by the factor C_H (Gnedin & Ostriker 1997; Madau et al. 1999; Miralda-Escudé et al. 2000; Miralda-Escudé 2003; Kohler et al. 2007).

We correct the recombination time for density variations scaled to a fiducial C_H ≈ 3, found in the simulations described below. At z ≈ 7, the IGM filaments have electron density n_e ≈ (10⁻⁴ cm⁻³)[(1 + z)/8]³C_H, and the characteristic times for hydrogen recombination and Hubble expansion are

$$\begin{eqnarray} t_{\rm rec} &\approx & (386\, {\rm Myr})(3/C_{\rm H}) T_4^{0.845} \left[ \frac{(1+z)}{8} \right]^{-3} \;, \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} t_{\rm H} &\approx & \big[H_0 \Omega _m^{1/2} (1+z)^{3/2}\big]^{-1} \approx (1.18\, {\rm Gyr}) \left[ \frac{(1+z)}{8} \right]^{-3/2}\!\!{.}\quad \end{eqnarray} \tag{ 7 }$$

In our calculations, we express the reionization criterion as $N_{\rm H} ({\rm Mpc}^{-3}) = \dot{\rho }_{\rm crit} t_{\rm rec} Q_{\rm LyC} f_{\rm esc}$, where $\dot{\rho }_{\rm crit}$ is the critical SFR density (M_☉ yr⁻¹ Mpc⁻³) and Q_LyC is the conversion factor from $\dot{\rho }_{\rm SFR}$ to the LyC production rate (see Section 2.2). We define f_esc as the fraction of LyC photons that escape from their galactic sources into the IGM (Dove & Shull 1994). Recent statistical estimates (Nestor et al. 2011) suggest that f_esc ≈ 0.1 for an ensemble of 26 Lyman-break galaxies and 130 Lyα emitters at z ≈ 3.09 ± 0.03, and it could be higher for the lower-mass galaxies that likely dominate the escaping LyC at z > 6 (Fernandez & Shull 2011). The LyC production efficiency, Q_LyC, is expressed in units 10⁶³ photons M⁻¹_☉ of star formation, since typical massive stars emit (1–10) × 10⁶³ LyC photons over their lifetime. To evaluate Q_LyC, we convert the SFR (by mass) into numbers of OB stars and compute the total number of ionizing photons produced by a star of mass m over its lifetime. We then integrate over an IMF, Ψ(m) = Km^−α, with a range m_min < m < m_max. The standard mass range is 0.1–100 M_☉, but changes in the mass range and IMF slope will affect the LyC production substantially. For example, differences in the IMF have been associated with higher Jeans masses in low-metallicity gas in the high-redshift IGM (Abel et al. 2002; Bromm & Loeb 2003), and Tumlinson (2007) and Smith et al. (2009) noted the potential influence of CMB temperature on modes of high-redshift star formation. Consequently, most cosmological simulations or calculations include a metallicity-induced IMF transition (Trenti & Shull 2010) between low-metallicity Population III star formation and metal-enhanced Population II.

The number of LyC photons produced per total mass in star formation depends on the IMF of the stellar population and is given by the conversion factor Q_LyC. We calculate this conversion by integrating the total number of LyC photons produced over the entire mass in star formation:

$$\begin{equation} Q_{\rm LyC} \equiv \frac{N_{\rm LyC}}{\dot{\rho }_{\rm SFR}t_{\rm rec}} = \frac{\int _{m_{{\rm OB}}}^{m_{{\rm max}}} \Psi (m) Q(m) \,dm}{\int _{m_{{\rm min}}}^{m_{{\rm max}}} \Psi (m) m \, dm} \;. \end{equation} \tag{ 8 }$$

Here, Q(m) is the lifetime-integrated number of LyC photons as a function of mass calculated from stellar atmosphere models and evolutionary tracks. The IMF, Ψ(m) = Km^−α, is integrated over the mass range m_min < m < m_max, where m is expressed in solar units. In the SFR simulator discussed in Section 3.3, the user can choose between a normal and broken IMF power law. The lower integration limit in the numerator, m_OB, is the mass at which stars no longer produce significant amounts of LyC photons.

In our calculator, we use and compare two models that calculate Q(m). First, using stellar atmospheres and evolutionary tracks (R. S. Sutherland & J. M. Shull, unpublished), we find that, over its main-sequence and post-main-sequence lifetime, an OB star of mass m produces a total number N_LyC = Q₆₃(m) × 10⁶³ of ionizing photons, where Q₆₃ ≈ 1–10 over the mass range m = 30–100 and for metallicities Z = 0.02–2.0 Z_☉. We have fitted our results to the form Q₆₃(m) ≈ Am − B, for m ⩾ m_OB = B/A (the mass m_OB defines the effective lower limit for stars that produce significant numbers of LyC photons). For metallicities in the range 0.002 ⩽ Z ⩽ 0.04 (where Z = Z_☉ = 0.02 is the solar metallicity) and for masses 15 < m < 100, the fitted coefficients are B = 1.578 and A(Z) = 0.0950 − 1.059Z. Thus, for Z = Z_☉ ≈ 0.02, we have A = 0.0738 and m_OB = 21.4 M_☉, while for Z = 0.1 Z_☉ = 0.002 we have A = 0.0929 and m_OB = 17 M_☉. Inserting these values for Q(m) into Equation (8), we integrate these LyC photon yields over the IMF to derive the conversion coefficient, Q_LyC, from total mass in star formation to total number of LyC photons produced (in units of 10⁶³ photons M⁻¹_☉).

The integrals in Equation (8) can be done analytically as functions of the IMF parameters and metallicity. For a Salpeter IMF (α = 2.35, 0.1 ⩽ m ⩽ 100) we find Q_LyC = 0.00236 (for Z = Z_☉), Q_LyC = 0.00366 (for Z = 0.2 Z_☉), and Q_LyC = 0.00384 (for Z = 0.1 Z_☉). Adopting a typical (low-Z) value Q_LyC = 0.004 (4 × 10⁶⁰ photons M⁻¹_☉), we can solve for the critical SFR for reionization, scaled to clumping factor C_H = 3, escape fraction f_esc = 0.2, and gas temperature T₄ = 1:

$$\begin{eqnarray} \dot{\rho }_{\rm crit} &=& (0.018 \,M_{\odot } \,{\rm yr}^{-1}\, {\rm Mpc}^{-3}) \left[ \frac{(1+z)}{8} \right]^{3}\nonumber\\ && \times \left( \frac{C_{\rm H}/3}{f_{\rm esc}/0.2} \right) \left( \frac{0.004}{Q_{\rm LyC}} \right) T_4^{-0.845} \;. \end{eqnarray} \tag{ 9 }$$

Our chosen value of C_H ≈ 3 is consistent with recent downward revisions (Pawlik et al. 2009) and with our numerical simulations discussed in Section 2.3. Escape fractions f_esc ≈ 0.1–0.2 have been inferred from observations at z ≈ 3 (Shapley et al. 2006; Nestor et al. 2011) and theoretical expectations (Fernandez & Shull 2011). Equation (9) agrees with prior estimates (Madau et al. 1999) when adjusted for our new scaling factors, in particular the ratio C_H/f_esc = 15. Their earlier paper assumed C_H = 30, f_esc = 1, Ω_bh² = 0.020, Q_LyC = 0.005, and z = 5. Our fiducial redshift has increased from z = 5 to z = 7, appropriate for the new discoveries of high-redshift galaxy candidates for reionization. Expressed with the same coefficients in Equation (9), their coefficient would be nearly the same, 0.020 M_☉ yr⁻¹ Mpc⁻³ at z = 7. One of the improvements in our formulation is to better identify the dependences on the physical parameters (C_H, f_esc, T_e) and the SFR-to-LyC conversion factor (Q_LyC), which can change with different IMFs and atmospheres.

A related calculation is the production rate of ionizing (LyC) photons per unit volume, needed to balance hydrogen recombinations. With the same assumptions as above, this is

$$\begin{eqnarray} \left(\frac{d {\cal N}_{\rm LyC}}{dt} \right)_{\rm crit} &=& n_{e} \, n_{\rm H\,\mathsc{ii}} \, \alpha _{\rm H}^{(B)} = (4.6 \times 10^{50} {\rm \,s}^{-1} {\rm \,Mpc}^{-3})\nonumber\\ && \times \left[ \frac{(1+z)}{8} \right]^3 T_{4}^{-0.845} \left(\frac{C_{\rm H}}{3} \right) \;. \end{eqnarray} \tag{ 10 }$$

For standard high-mass stars (O7V, solar metallicity) each with LyC photon luminosity 10⁴⁹ s⁻¹, this rate corresponds to an O-star space density n_O ≈ 50 Mpc⁻³ at z = 7. Schaerer (2002, 2003) also computed models of Population III and low-metallicity stars based on non-LTE model atmospheres and new stellar evolution tracks and evolutionary synthesis models. The lifetime total number of LyC photons produced per star of mass m is $Q(m) = \bar{Q}(H) \times t_\ast$. For stars of mass parameter x = log (m/M_☉), the number of ionizing photons, $\bar{Q}$, emitted per second per star, is given by

$$\begin{equation} \hbox{\fontsize{7}{8}\selectfont ${\rm log}_{10}[\bar{Q}_{\rm H}/{\rm s}^{-1}]= \left\lbrace \begin{array}{@{}lll}43.61 + 4.90x - 0.83x^2&\!\! Z = 0,&\!\! 9\hbox{--}500 \,M_\odot, \\ \ms 39.29 + 8.55x&\!\! Z=0,& \!\! 5\hbox{--}9 \,M_\odot,\\ \ms 27.80+30.68x-14.80x^2+2.50x^3&\!\! Z = 0.02 \,Z_\odot, &\!\! 7\hbox{--}150 \,M_\odot, \\ \ms 27.89+27.75x-11.87x^2+1.73x^3&\!\! Z= Z_\odot,&\!\! 7\hbox{--}120 \,M_\odot \\ \ms \end{array} \right.$} \end{equation} \tag{ 11 }$$

and the star's lifetime t_* is given by

$$\begin{equation} \hbox{\fontsize{7}{8}\selectfont ${\rm log}_{10}\left[t_\ast \right]= \left\lbrace \begin{array}{@{}ll}9.785-3.759x+1.413x^2-0.186x^3&Z = 0, \\ \ms 9.59-2.79x+0.63x^2&Z = 0.02 \,Z_\odot, \\ \ms 9.986-3.497x+0.894x^2&Z = Z_\odot. \end{array} \right.$} \end{equation} \tag{ 12 }$$

Sample values of Q_LyC and $\dot{\rho }_{\rm crit}$ for different IMFs are shown in Table 1, for both model atmospheres. One can obtain different values of Q_LyC for power-law IMFs by varying their high-mass slope β and the minimum and maximum masses. Increases in Q_LyC translate into decreases in critical SFR. For example, at 0.1 Z_☉, if the minimum mass is raised to 1 M_☉ with a Salpeter slope (β = 2.35), the LyC production factor rises to Q_LyC = 0.0098, some 2.5 times higher than the equivalent model, Q_LyC = 0.00384 for m_min = 0.1 M_☉. If the IMF is flatter, β = 2 instead of 2.35, with minimum mass fixed at m_min = 0.1 M_☉, Q_LyC = 0.0177 (4.6 times higher). At the upper end of the IMF, if one increases m_max from 100 to 200 M_☉, keeping β = 2.35 and m_min = 0.1, one finds Q_LyC = 0.0054 (1.4 times higher). All of these variations can be explored with our reionization/SFR simulator, described in Section 3.3.

Table 1. Sample Values of Q_LyC and $\dot{\rho }_{\rm crit}$ for Different IMFs

Mmin	Mmax	α	Z	QLyCa	QLyCb	$\dot{\rho }_{\rm crit}$a	$\dot{\rho }_{\rm crit}$b
(M☉)	(M☉)
0.1	100	2.35	Z☉	0.00236	0.00286	0.0306	0.0253
0.1	100	2.35	0.02 Z☉	0.00397	0.00558	0.0181	0.0129
0.1	100	2.35	0	0.00401	0.00752	0.0180	0.0096
0.1	100	2.35	0.2 Z☉	0.00365	⋅⋅⋅	0.0197	⋅⋅⋅
0.1	100	2.35	0.1 Z☉	0.00383	⋅⋅⋅	0.0188	⋅⋅⋅
1.0	100	2.35	0.1 Z☉	0.00976	⋅⋅⋅	0.0074	⋅⋅⋅
0.1	100	2.00	0.1 Z☉	0.01267	⋅⋅⋅	0.0057	⋅⋅⋅
0.1	200	2.35	0.1 Z☉	0.00543	⋅⋅⋅	0.0133	⋅⋅⋅

Notes. Production efficiency of ionizing (LyC) radiation, Q_LyC, in units of 10⁶³ photons M⁻¹_☉. Critical SFR, $\dot{\rho }_{\rm crit}$, is given in units of M_☉ yr⁻¹ Mpc⁻³ (comoving) assuming C_H = 3, f_esc = 0.2, and T₄ = 1. ^aSutherland & Shull unpublished model atmospheres. ^bSchaerer model atmospheres.

Download table as:  ASCIITypeset image

The critical SFR of 0.018 M_☉ yr⁻¹ Mpc⁻³ can be understood from simple arithmetic. Over the 386 Myr recombination time (Equation (6)) with C_H = 3, f_esc = 0.2, and z = 7, approximately 20 million stars are formed per comoving Mpc³ in a Salpeter IMF (0.1–100 M_☉) with mean stellar mass 0.352 M_☉. Of these stars, a small fraction, f_OB ≈ 7 × 10⁻⁴, are massive OB stars (m > 20 M_☉). At an average of 2 × 10⁶³ LyC photons and escape fraction f_esc = 0.2, these ∼14,000 OB stars will produce a net (escaping) ∼5 × 10⁶⁶ LyC photons over their lifetime. These photons are sufficient to ionize N_H = 5.6 × 10⁶⁶ hydrogen atoms Mpc⁻³.

The simulations used in this study were performed with Enzo, an Eulerian adaptive mesh-refinement (AMR), hydrodynamical + N-body code (Bryan & Norman 1999; O'Shea et al. 2004, 2005). Smith et al. (2011) enhanced this code by adding new modules for star formation, primordial chemistry, and cooling rates consistent with ionizing radiation, metal transport, and feedback. The ionizing background is spatially constant and optically thin, but variable in redshift. At this stage, we have not implemented radiative transfer or spectral filtering by the IGM.

For the clumping factor calculation, we ran a simulation on a 50 h⁻¹ Mpc static grid ("unigrid") cube with 1024³ cells, denoted as run 50_1024_2 in Table 1 of Smith et al. (2011). The radiative heating from the ionizing background plays an important role in determining the properties of the filamentary structure. As a filament is ionized and heated, its density drops and its temperature rises; both effects reduce the recombination rate. To study these effects, we ran four moderate-resolution (50 h⁻¹ Mpc unigrid cube with 512³ cells) simulations with varying ionizing backgrounds, summarized in Table 2. After initial submission of this paper, we ran a 1536³ simulation to check convergence and assess the "cosmic variance" among eight sub-volumes of the 1024³ and 1536³ simulations. The standard UV background was taken from Haardt & Madau (2001), although we also explore new computations of high-z SFRs by Trenti et al. (2010) and Haardt & Madau (2012). These four simulations were (1) no photoionizing background, (2) UV background ramped up from z = 7 to z = 6 (run 50_512_2 from Smith et al. 2011), (3) UV background ramped up from z = 9 to z = 8, and (4) UV background ramped up from z = 9 to z = 8, but twice as strong as in (3). Post-processing of the simulations was performed using the data analysis and visualization package yt,³ documented by Turk et al. (2011).

Table 2. Parameters of Simulations

Run	l	N1/3cells	zUV	Relative Strength
	(h−1 Mpc)			of Background
$50\_1024\_7$a	50	1024	7	1
$50\_512\_7$b	50	512	7	1
$50\_512\_9$	50	512	9	1
$50\_512\_9\_2$	50	512	9	2
$50\_512\_0$	50	512	N/A	0
$50\_1536\_9$	50	1536	9	1

Notes. z_UV is the redshift at which the ionizing background radiation (Haardt & Madau 2001) is turned on, ramping up to a constant value by z = z_UV − 1. ^a,bSimulations $50\_1024\_2$ and $50\_512\_2$ from Smith et al. (2011).

Download table as:  ASCIITypeset image

For regions of ionized hydrogen of density $n_{\rm H\,\mathsc{ii}}$, we calculate the clumping factor, C_H, using two different methods. The first calculation, which has been used in some earlier studies, uses density weighting. In this "density field" (DF) method we define

$$\begin{equation} C_{\rm DF} = \frac{ \big\langle n_{\rm H\,\mathsc{ii}}^2 \big\rangle }{ \langle n_{\rm H\,\mathsc{ii}} \rangle ^2} \; \end{equation} \tag{ 13 }$$

and average the density fields (quantities x_i, denoting either $n_{\rm H\,\mathsc{ii}}$ or $n_{\rm H\,\mathsc{ii}}^2$), where parameter averages are computed by summing over grid cells (j), subject to various "cuts" on the IGM overdensity and gas temperature and weighted by factors, w_j,

$$\begin{equation} \left<x_i\right> =\displaystyle \sum \limits _j x_j w_j \big/ \displaystyle \sum \limits _j w_j \;. \end{equation} \tag{ 14 }$$

In the second calculation, we compare the local recombination rate to the global average recombination rate, averaged over density and temperature in cells,

$$\begin{equation} C_{\rm RR} = \frac{\big\langle n_{e} \, n_{\rm H\,\mathsc{ii}} \; \alpha _{\rm H}^{(B)} (T) \big\rangle }{ \langle n_{e} \rangle \langle n_{\rm H\,\mathsc{ii}} \rangle \, \big\langle \alpha _{\rm H}^{(B)} (T) \big\rangle } \;, \end{equation} \tag{ 15 }$$

where again 〈x〉 denotes a weighted average and α^(B)_H(T) is the case-B radiative recombination rate coefficient for hydrogen, as tabulated by Osterbrock & Ferland (2006). We refer to this as the "recombination rate" (RR) method.

Because the purpose of the clumping factor is to correct for an enhanced recombination rate, we believe C_RR to be a more appropriate representation. The recombinations that are important to removing ionizing photons occur in the filamentary structure of the IGM, and the clumping factor should only be calculated in these structures. To assess the critical SFR necessary to maintain an ionized medium, we focus on grid cells that are significantly ionized. Because a negligible amount of recombination occurs in cells containing mostly neutral gas, these cells should not contribute to the clumping factor. If we do not exclude neutral gas, the clumping factor is extremely high (C_H ∼ 100) before the ionizing background is turned on. We find large density gradients between the neutral gas (uniformly distributed over the simulation) and the ionized gas. Because the clumping factor is a measure of the inhomogeneity of the medium, a large density gradient yields a large value of C_H. If low-density voids are included in the calculation, the larger density gradient leads to an overestimate of the clumping. By setting both upper and lower density thresholds, we can exclude collapsed halos and low-density voids from our calculations.

Previous studies (Miralda-Escudé et al. 2000; Miralda-Escudé 2003; Pawlik et al. 2009) addressed these issues by setting a density threshold that excludes collapsed halos, but they did not set a lower limit to exclude voids from their calculations. To explore these effects in a filamentary IGM, we make various cuts of our data in baryon overdensity ($\Delta _b \equiv \rho _b/\bar{\rho }_b$) and in temperature, metallicity, and hydrogen ionization fraction ($x \equiv n_{\rm H\,\mathsc{ii}}/n_{\rm H}$). In our standard formulation, we include only those cells that meet the following criteria: 1 < Δ_b < 100, 300 K < T < 10⁵ K, Z < 10⁻⁶ Z_☉, and x > 0.05. We believe this "data cut" adequately represents unenriched IGM lying on an adiabat (see Figure 19 of Smith et al. 2011). We also explore the clumping factor with no lower density threshold (total range Δ_b < 100). Our results, presented in Section 3.1, show a small increase in C_H(z) from the wider range in densities by including low-density cells with Δ_b < 1. However, these low-density voids do not contribute substantially to the recombination rate.

In summary, our prescriptions for calculating the clumping factor yield a more physical representation of the enhanced recombination rate that is an important component to many reionization models. By not assuming a fully ionized medium and by specifically following $n_{\rm H\,\mathsc{ii}}$, we are able to exclude denser neutral gas that does not contribute appreciably to the recombination rate. We make simple assumptions on the reionization process and reionization history (Section 3.2), turning on the ionizing radiation field at redshifts z = 7 or z = 9 and following the thermal history arising from photoelectric heating, radiative cooling, and pressure smoothing (Pawlik et al. 2009). The metal-line and molecular cooling, metal transport, and feedback included in our simulations also allow us to accurately represent the thermodynamics of the gas, which has been shown to have a significant effect on the evolution of the clumping factor. Future models will include discrete sources and radiative transfer, accounting for temperature increases arising from photoheating with spectral hardening (Abel & Haehnelt 1999). Because our current simulations employ a spatially constant ionizing radiation field, we anticipate carrying out these more realistic situations.

Early studies of IGM clumping adopted high values, C_H > 40 at z < 5 (Gnedin & Ostriker 1997). As discussed earlier, we believe these values are too high for the ionized IGM filaments, which have expanded as a result of the heat deposited by LyC photons. The differences between observations and inferred critical SFR densities can largely be attributed to this high clumping factor (Sawicki & Thompson 2006; Bouwens et al. 2007). More recent studies have trended toward less clumping. Raičević & Theuns (2011) argue that using a global clumping factor overestimates the recombination rate, and that local values should be used instead. In this study, we calculate the clumping factor for a series of high-resolution cosmological simulations for ionized hydrogen and helium that explore how the photoheating of an ionizing background affects the IGM thermodynamics (density and temperature). We impose criteria in overdensity, temperature, metallicity, and ionization fraction to constrain our calculations to IGM filaments where recombinations and ionization fronts are most important. The clumping in simulation $50\_1024\_7$ is calculated as a function of redshift from z = 15 to z = 0 for both C_DF and C_RR, weighted by overdensity and by volume. Note that the unigrid simulations automatically yield volume weighting, when averaged over cells. For the remainder of our simulations, we find a power-law overdensity distribution, with an average form f(Δ_b) = 10⁶ Δ^−1.5_b. Weighting these cells by overdensity gives undue emphasis on small numbers of high-density cells. This method is also biased, since it does not calculate the clumping factor where most of the recombinations are occurring.

Figure 1 compares the DF and RR methods and shows the difference between the different weights used in our calculations. In both methods, weighting by overdensity results in higher clumping factors. The RR method (Equation (15)) results in a lower clumping factor because it accounts for both density and temperature effects from the ionizing background on the clumping averages. Photoelectric heating during photoionization causes the filaments to expand, lowering the density. The increased electron temperature also reduces the radiative recombination rate coefficient. The DF method (Equation (13)) includes no dependence on the recombination rate coefficient and therefore only counts the effect of photoheating on the filament density. See Section 3.2 for further discussion.

Zoom In Zoom Out Reset image size

In our current simulations, we turn on a spatially constant UV background at z = 7 or z = 9. The effects of different "density cuts" in the summation over cells (see Section 2.3) are shown in Figure 2. We find little difference at z > 9 (turn-on of ionizing radiation) and see a small increase in C_H at z < 9 when we include the lower density cells with Δ_b < 1. In future, more refined simulations with radiative transfer, we expect photoionization to be initiated in high-density regions, where the stars are located. Once the UV background is turned on, lower-density regions are easily photoionized. The photoelectric heating of the UV background causes the filaments to expand and become diffuse, resulting in a lower clumping factor. The elevated temperature also reduces the recombination rate coefficient, which in turn lowers C_H when calculated by the RR method. When calculating a mean recombination rate, only the recombining ionized gas is relevant. By z ∼ 6, the IGM is nearly completely ionized and only small regions of H i remain. As noted earlier, in more refined simulations, the remaining H i would be largely self-shielded from the UV background and likely to reside in clumps of high density yet unreached by ionization fronts.

Zoom In Zoom Out Reset image size

Figure 3 explores the convergence and cosmic variance among different simulations. In two panels, we compare results from our 512³ simulations with larger simulations with 1024³ and 1536³ cells. In each panel, we show the average clumping factors, together with those in eight sub-volumes. The average values of C_H agree in the 1024³ and 1536³ simulations, and the small (±10%) variance among these sub-volumes suggests that the 1024³ and 1536³ simulations are converged. Over the redshift range 5 < z < 9 in the 1536³ simulation, the clumping factor is well fitted by a power law,

$$\begin{equation} C_{\rm H}(z) = (2.9) \left[ \frac{(1+z)}{6} \right] ^{-1.1} \;, \end{equation} \tag{ 16 }$$

representing a slow rise in clumping to lower redshift, after the turn-on of ionizing radiation.

Zoom In Zoom Out Reset image size

To study the effect of radiative heating from the ionizing background on filamentary structure, we ran a suite of moderate-resolution simulations with different ionizing backgrounds (see Table 2 for details). The redshift at which the UV background is turned on affects the clumping factor, causing it to drop and then recover at lower redshift. Pawlik et al. (2009) attribute this effect to Jeans filtering, where the photoheating raises the cosmological Jeans mass, preventing further accretion onto low-mass halos and smoothing out small-scale density fluctuations. The photoheating also heats the filaments we are focusing on, which lowers the hydrogen recombination rate coefficient, α^(B)_H∝T^−0.845, and causes the filaments to expand; both of these effects reduce the clumping factor. Without a UV background, the clumping continues to rise as the filaments gravitationally collapse. Pawlik et al. (2009) also claim that the clumping factor at z = 6 is insensitive to when the background is turned on, as long as it is turned on at z > 9. When an ionizing background is introduced, photoheating acts as a positive feedback to reionization by lowering the clumping factor and making it easier to stay ionized. This same photoheating mechanism also suppresses star formation in low-mass halos, which in turn lowers the ionizing photon production rate by star-forming galaxies and acts as a negative feedback to reionization (Pawlik et al. 2009). These thermodynamic processes affect clumping and structure and emphasize the importance of carefully modeling the strength of feedback processes and their effects on Jeans mass.

We have explored what happens when the ionizing radiation turned on earlier, especially during the interval z = 6–8 marking the transition from a neutral to fully ionized IGM. In a series of 512³ simulations, we explore the influence of turn-on of photoionizing radiation between redshifts z = 7 and z = 9. Figure 4 shows C_H, computed for the moderate-resolution simulations via the DF method and weighted by volume. We do not plot simulation $50\_512\_9\_2$, since it is identical to simulation $50\_512\_9$. As in the results of Pawlik et al. (2009), we find that photoheating from the ionizing background results in a decrease in the clumping factor. The clumping factors fall along two tracks: a higher track at redshifts above the turn-on of the UV background, and a lower one after turn-on. By z = 5, we find nearly identical clumping factors for both backgrounds. After a substantial recovery time, the redshift when the background is turned on is not important. Before this recovery, the clumping factor of the earlier background is lower by a factor of ∼2, resulting in earlier reionization.

Zoom In Zoom Out Reset image size

One can compare the strengths of feedback to reionization, where "positive feedback" lowers the clumping factor and "negative feedback" suppresses star formation. At z = 5 the total stellar mass (SFR density) of simulation $50\_512\_7$ is 1.19 (1.13) times lower than that of simulation $50\_512\_0$, while the clumping factor is 1.64 times lower. For simulation $50\_512\_9$, the total stellar mass (SFR density) and clumping factor are, respectively, 1.52 (1.49) and 1.66 times lower than those of simulation $50\_512\_0$. This suggests that the positive feedback introduced by the background is greater than the negative feedback. However, the stellar mass (SFR density) does not recover in the same manner as the clumping factor. Once photoheating suppresses the formation of small-mass halos, the Hubble flow takes over and prevents them from collapsing and forming stars. Therefore, the redshift at which the background is turned on determines whether the positive or negative feedback dominates and whether the ionizing background will cause reionization to be accelerated or delayed.

In connection with this project, we have developed a user interface for calculating the critical SFR density $\left(\dot{\rho }_{\rm crit}\right)$ needed to maintain the IGM ionization at a given redshift. The software computes the effects of variations in the stellar IMF (slope, mass range) and model atmospheres, as well as the redshift evolution of metallicity and gas thermodynamics (density, temperature, coupling to the CMB). The clumping factor and LyC escape fraction are free parameters in the calculator. Our simulations find ranges of C_H ≈ 1–10 depending on redshift, overdensity, and thermal phase of the (photoionized or shock-heated) IGM. For the new 1536³ simulation, the global mean clumping factor is 〈C_H〉 ≈ 3, with a power-law fit C_H(z) = (2.9)[(1 + z)/6]^−1.1 for redshifts between 5 < z < 9.

This calculator is a useful tool for determining the population of galaxies responsible for reionization. The critical SFR per comoving volume (Equation (9)) is obtained by balancing the production rate of LyC photons with the number of hydrogen recombinations. Here, C_H is the clumping factor of ionized hydrogen, f_esc is the escape fraction of LyC photons from their host galaxies, and T₄ is the temperature scaled to 10⁴ K. The conversion factor, Q_LyC, from stellar mass to total number of LyC photons produced, is regulated by the IMF and model atmospheres. The user has the option of controlling these parameters to determine the resulting critical SFR density, subject to several observational constraints. This simulator can be accessed at http://casa.colorado.edu/~harnessa/SFRcalculator with an html interface for easy use.

Two standard tests of the SFR use the simulator to calculate the ionization histories of H ii and He iii and compare them to the ionized volume filling factor, $Q_{\rm H\,\mathsc{ii}}(z)$, and the CMB optical depth, τ_e(z). The calculator derives τ_e(z) by integrating the differential form of Equation (3) for various SFR histories and parameters. The average evolution of $Q_{\rm H\,\mathsc{ii}}$ is found by numerical integration of the rate equation (Madau et al. 1999) expressing the sources and sinks of ionized zones,

$$\begin{equation} \frac{dQ_{\rm H\,\mathsc{ii}}}{dt} = \frac{\dot{n}_{\rm LyC}}{\langle n_{\rm H} \rangle } - \frac{Q_{\rm H\,\mathsc{ii}}}{\langle t_{\rm rec}(C_{\rm H}) \rangle } \;. \end{equation} \tag{ 17 }$$

The source term, $\dot{n}_{\rm LyC} = \dot{\rho }_{\rm SFR} f_{\rm esc} Q_{\rm LyC}$, represents the net ionizing photon production rate, computed from selected models of the SFR density, $\dot{\rho }_{\rm SFR}$. Here, 〈n_H〉 is the mean hydrogen number density, and 〈t_rec〉 is the hydrogen recombination timescale, which depends on C_H as shown in Equation (6). The IGM is assumed to be fully ionized when $Q_{\rm H\,\mathsc{ii}} = 1$. We integrate a similar equation for $Q_{\rm He\,\mathsc{iii}}$ to follow He ii photoionization in the QSO 4 ryd continuum. We use the QSO emissivities at 1 ryd (Haardt & Madau 2012), extrapolated to 4 ryd assuming a spectrum with specific flux F_ν∝ν^−1.8.

Figure 5 shows ionization histories, quantified by $Q_{\rm H\,\mathsc{ii}}(z)$ and τ_e(z), for various values of clumping factor, escape fraction, IGM temperature, and LyC-production efficiencies (Q_LyC). We adopt the SFR history from Trenti et al. (2010) with an evolving luminosity function. Our standard model adopts C_H = 3, f_esc = 0.2, T₄ = 2, and Q_LyC = 0.004 (4 × 10⁶⁰ photons/M_☉). Different curves show the effects of changing C_H and f_esc, including two models in which these parameters evolve with redshift. With these parameters, we are typically able to complete reionization by z ≈ 7, consistent with observations of Gunn–Peterson troughs, redshift evolution in Lyα emitters, and IGM neutral fraction. The model with a constant f_esc = 0.05 does not complete the ionization until z ≈ 4, which is far too late. Thus, we be believe that LyC escape fractions must be considerably larger (⩾20%), perhaps evolving to higher values at z > 6 shown by green and magenta curves.

Zoom In Zoom Out Reset image size

Figure 6 compares the dependence of $Q_{\rm H\,\mathsc{ii}}(z)$ and τ_e(z) on SFR histories, computed with an evolving luminosity function (Trenti et al. 2010) and new calculations (Haardt & Madau 2012) of the SFR density, $\dot{\rho }_{\rm SFR}(z)$. Figure 7 compares the effects of two choices of model atmospheres, with a fixed SFR history from Trenti et al. (2010). The online calculator provides ionization fractions, $Q_{\rm H\,\mathsc{ii}}(z)$ and $Q_{\rm He\,\mathsc{iii}}(z)$, together with CMB optical depth, τ_e(z). In the simulator, users can select other IGM parameters and SFR histories. The main difference between the two SFR models lies in the assumptions at z > 8. Haardt & Madau (2012) rely on an empirical extrapolation of the SFR as a function of redshift, while Trenti et al. (2010) adopt a physically motivated model based on the evolution of the dark matter halo mass function. The two approaches are similar at z ≲ 8 but differ significantly at higher redshift, where an empirical extrapolation does not capture the sharp drop in the number density of galaxies observed at z ∼ 10 (see Figure 8 in Oesch et al. 2012). As a consequence, the reionization history from the Haardt & Madau (2012) model is more extended at high z, especially when the efficiency of reionization is increased because of evolving clumping factor and escape fraction (our preferred models, green and magenta curves in Figure 5). The two models yield quite different predictions for the duration of reionization, defined as the redshift interval, Δz, over which $Q_{\rm H\,\mathsc{ii}}$ evolves from 20% to 80% ionized. Our preferred SFR models (green and magenta lines in Figure 5) have Δz ≈ 3 (from z ≈ 10.5 to z ≈ 7.5), whereas the Haardt–Madau SFR histories exhibit a more extended interval, Δz ≈ 6 (from z ≈ 13 to z ≈ 7). This difference in Δz could be tested by upcoming 21 cm experiments; see Bowman & Rogers (2010) for an initial constraint, Δz > 0.06.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 8 illustrates a third constraint on the reionization epoch (Pritchard et al. 2010; Lidz et al. 2011) comparing SFR histories with estimates of the ionizing background at z = 5.5 ± 0.5 (Fan et al. 2006; Bolton & Haehnelt 2007; Haardt & Madau 2012). The LyC comoving emissivity (in photons s⁻¹ Mpc⁻³), defined as $\dot{n}_{\rm LyC} = \dot{\rho }_{\rm SFR} f_{\rm esc} Q_{\rm LyC}$, can be related to the ionizing background at z = 5–6, using recent estimates of the hydrogen photoionization rate, $\Gamma _{\rm H\,\mathsc{i}}(z)$, from Haardt & Madau (2012), and the LyC mean free path, $\lambda _{\rm H\,\mathsc{i}}$, from Songaila & Cowie (2010). The LyC emissivity is proportional to the SFR density, computed from our halo mass function model (Trenti et al. 2010), integrated down to absolute magnitudes M_AB = −18 (Bouwens et al. 2011a) or to M_AB = −10, the faint limit suggested by Trenti et al. (2010). We adopt a fiducial LyC production parameter Q_LyC = 0.004, corresponding to 4 × 10⁶⁰ LyC photons produced per M_☉ of star formation, and we use two different models for LyC escape fraction, f_esc (constant at 20% and varying with redshift). For quantitative values, we assume an ionizing background with specific intensity, J_ν = J₀(ν/ν₀)^−α (in units of erg cm⁻² s⁻¹ sr⁻¹ Hz⁻¹) with power-law index α ≈ 2 at energies above hν₀ = 1 ryd. The hydrogen photoionization rate is $\Gamma _{\rm H\,\mathsc{i}} = [4 \pi J_0 \sigma _0/ h (\alpha + 3)]$ for a hydrogen photoionization cross section σ_ν ≈ σ₀(ν/ν₀)⁻³ with σ₀ = 6.3 × 10⁻¹⁸ cm². For this spectrum, the frequency-integrated ionizing intensity is J_tot = J₀ν₀/(α − 1), and the LyC photon flux (photons cm⁻² s⁻¹) integrated over all solid angles is Φ_LyC = (4πJ₀/hα). Because we analyze the photon mean-free path, we normalize J₀ to Φ_LyC and $\Gamma _{\rm H\,\mathsc{i}}$,

$$\begin{equation} J_0 = \frac{h \alpha \Phi _{\rm LyC}}{4 \pi } = \frac{h (\alpha + 3)}{4 \pi \sigma _0} \; \Gamma _{\rm H\,\mathsc{i}} \;. \end{equation} \tag{ 18 }$$

By approximating Φ_LyC as the product of the LyC emissivity and mean free path, we arrive at the calibrations:

$$\begin{eqnarray} J_{\rm tot} &=& \frac{h \nu _0}{4 \pi \sigma _0} \frac{(\alpha +3)}{(\alpha - 1)} \; \Gamma _{\rm H\,\mathsc{i}} = (6.88\times 10^{-7} {\rm \,erg\, cm}^{-2} {\rm \,s}^{-1} {\rm \,sr}^{-1})\nonumber\\ && \times \left[ \frac{\Gamma _{\rm H\,\mathsc{i}}}{5 \times 10^{-13} {\rm \,s}^{-1}} \right], \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} \Phi _{\rm LyC} &=& \frac{(\alpha +3)}{\alpha } \frac{\Gamma _{\rm H\,\mathsc{i}}}{ \sigma _0} = (1.98 \times 10^5 {\rm \,photons \,cm}^{-2} {\rm \,s}^{-1})\nonumber\\ && \times \left[ \frac{\Gamma _{\rm H\,\mathsc{i}}}{5 \times 10^{-13} {\rm \,s}^{-1}} \right] \;. \end{eqnarray} \tag{ 20 }$$

Finally, we relate Φ_LyC to the star formation rate $\dot{\rho }_{\rm SFR}$ and ionization rate $\Gamma _{\rm H\,\mathsc{i}}$,

$$\begin{eqnarray} \dot{\rho }_{\rm SFR} &=& (0.056 \,M_{\odot } \,{\rm yr}^{-1}\, {\rm Mpc}^{-3}) \left[ \frac{\Gamma _{\rm H\,\mathsc{i}}}{5 \times 10^{-13} {\rm \,s}^{-1}} \right] \left[ \frac{0.1}{f_{\rm esc}} \right]\nonumber\\ && \times \left[ \frac{0.004}{Q_{\rm LyC}} \right] \left[ \frac{9.8 {\rm \,pMpc}}{\lambda _{\rm H\,\mathsc{i}}} \right] \left[ \frac{6.5}{1+z} \right]^{3} \;. \end{eqnarray} \tag{ 21 }$$

Here, we have scaled the parameters to the same values assumed by Lidz et al. (2011), namely, $\Gamma _{\rm H\,\mathsc{i}} = 5 \times 10^{-13} {\rm \,s}^{-1}$, f_esc = 0.1, an ionizing spectral index α = 2, and redshift z = 5.5, at which Songaila & Cowie (2010) fit $\lambda _{\rm H\,\mathsc{i}} \approx 9.8$ proper Mpc. Our parameter Q_LyC = 0.004 corresponds to the Lidz et al. (2011) LyC photon production calibration, 10^53.1 s⁻¹ per M_☉ yr⁻¹ of star formation. However, our coefficient, 0.056 M_☉ yr⁻¹ Mpc⁻³, is slightly larger than their value, 0.039 M_☉ yr⁻¹ Mpc⁻³, an effect that may arise from our different method of relating $\Gamma _{\rm H\,\mathsc{i}}$ to the ionizing radiation field and SFR.

Zoom In Zoom Out Reset image size

More careful examination of these parameters suggests that, at z ≈ 5.5, the hydrogen ionization rate $\Gamma _{\rm H\,\mathsc{i}} \approx 3.6 \times 10^{-13} {\rm \,s}^{-1}$, given in Table 3 of Haardt & Madau (2012), and the observed mean free path, $\lambda _{\rm H\,\mathsc{i}}$, at z = 5.5 may be closer to 6 proper Mpc (see Figure 10 of Songaila & Cowie 2010). Rescaling to those two parameters, we find a similar coefficient of 0.066 M_☉ yr⁻¹ Mpc⁻³. This SFR density at z = 5.5 ± 0.5 (see also Figure 8) is comparable to the critical value needed to maintain reionization at z = 7, but the observed rates and mean free paths are declining rapidly with redshift. All three constraints on SFR suggest that full reionization is more likely to occur at redshift z_rei ≈ 7 than at z_rei = 10.