CHEMICAL ENRICHMENT OF DAMPED Lyα SYSTEMS AS A DIRECT CONSTRAINT ON POPULATION III STAR FORMATION

The mean mass growth rate of dark matter halos in the ΛCDM cosmology can be parameterized by the fitting function (Fakhouri et al. 2010)

$$\begin{eqnarray} \dot{M}(z) &=& 46.1\,M_\odot \,\mathrm{yr}^{-1}\,(1+1.1z) \nonumber \\ &&\times \sqrt{\Omega _m(1+z)^3+\Omega _\Lambda }\,\left(\frac{M(z)}{10^{12}\,M_\odot }\right)^{1.1}, \end{eqnarray} \tag{ 1 }$$

where M(z) denotes the mean mass in the main branch of the merger tree of a halo at redshift z. This growth rate corresponds to the mean rate of growth of halo mass. This relation can be understood in terms of the extended Press–Schechter picture of structure formation (Neistein & Dekel 2008). By comparing it with cosmological simulations, Fakhouri et al. (2010) have shown that it is accurate over a mass range spanning at least five orders of magnitude (10¹⁰–10¹⁵ M_☉ at z = 0) and a wide range of redshifts (0 ⩽ z ≲ 15). Integrating Equation (1) with the boundary condition, M(z = 0) = M_h gives the mean mass of a halo at any redshift. Note that this describes the evolution of the dark matter mass of the halo.

Baryonic content of galaxies is influenced by gas inflows, gas outflows, and star formation. Hydrodynamical simulations suggest that these processes tend to be in equilibrium (Davé et al. 2012). We now describe these processes and calculate the mean evolution of gas mass (M_g), stellar mass (M_*), and metal mass (M_Z) in each halo mass bin. The metal mass term describes the mass of the metals in the gas phase (ISM); it excludes the metals that are locked up inside stars.

2.2.1. Gas

Gas inflows replenish the gas content of halos as they evolve. We assume that the mean gas inflow rate is proportional to the mean halo growth rate

$$\begin{equation} \dot{M}_{g,\mathrm{in}}(z) = f_{g,\mathrm{in}}\left(\frac{\Omega _b}{\Omega _m}\right)\dot{M}(z). \end{equation} \tag{ 2 }$$

We follow Bouché et al. (2010) and take the constant of proportionality to be f_{g, in} = 0.7, except when the halo mass is below the filtering mass M_min (described below), in which case there is no gas inflow and we set f_{g, in} = 0. The choice of these values is motivated by the comparison of measured SFR in massive high-redshift galaxies with predicted gas accretion rates (Genel et al. 2008; Dekel et al. 2009; Bauermeister et al. 2010).

Once inside the halo, some of the accreted gas is lost to star formation. The amount of gas lost is

$$\begin{eqnarray} \dot{M}_{g, \mathrm{sf}}(z) &=& {-}\psi (z) +\int _{m_l}^{m_u}dm \,\phi (m)\cdot \psi [t(z)-\tau (m)] \nonumber\\ && \cdot [m-m_r(m)], \end{eqnarray} \tag{ 3 }$$

where ψ is the SFR (which we describe below), and the second term accounts for mass loss from evolving stars via stellar winds and supernova explosions. In the second term, ϕ(m) is the stellar IMF, m_r(m) is the remnant mass left by a star with initial mass m, τ(m) is its lifetime, and t(z) is the cosmic time. The lower limit m_l of the integral is such that τ(m_l) = t(z). (The upper limit m_u corresponds to the high-mass limit of the stellar IMF, which we discuss below.) The instantaneous recycling approximation (IRA) is often used in calculating the mass loss from stars (e.g., Krumholz & Dekel 2012). In this approximation, stars with initial mass above a certain value m₀ are assumed to die instantaneously (τ = 0), while those with initial mass less than m₀ are assumed to live forever (τ = ∞). However, although this is a reasonable assumption to make while studying, say, the solar neighborhood, it is not a good approximation at high redshift, when stellar lifetimes are comparable to the Hubble time. Therefore, we do not use IRA in Equation (3).

Gas is also lost due to outflows resulting out of supernova-driven winds. This effect is expected to be proportional to the SFR and inversely proportional to the depth of the halo potential well (Daigne et al. 2006; Erb 2008; Bouché et al. 2010; Krumholz & Dekel 2012). We write the mass loss due to outflows as

$$\begin{eqnarray} \dot{M}_{g, \mathrm{out}}(z) &=& {-} \frac{2\epsilon }{v^2_\mathrm{esc}}\int _{m_l}^{m_u}dm\,\left\lbrace \phi (m)\right.\nonumber \\ && \times\left. \psi [t(z)-\tau (m)]\cdot E_\mathrm{kin}(m)\right\rbrace, \end{eqnarray} \tag{ 4 }$$

where $v^2_\mathrm{esc}=2GM/R_\mathrm{vir}$ is the escape velocity of the halo, and E_kin(m) is the kinetic energy released by a star of mass m. The parameter is fixed by calculating the total baryon fraction in structures, f_{b, struct}, at z = 0 (Fukugita & Peebles 2004).

Combining the contribution of the above three processes of star formation, inflows, and outflows, the evolution of the gas mass of a halo can be written as

$$\begin{equation} \dot{M}_g(z) = \dot{M}_{g,\mathrm{in}}(z) + \dot{M}_{g, \mathrm{sf}}(z) + \dot{M}_{g, \mathrm{out}}(z), \end{equation} \tag{ 5 }$$

where each term is redshift and halo mass dependent.

2.2.2. Stars

We assume that the SFR, ψ, in a halo tracks the total amount of cold gas, M_cool, inside that halo. Thus,

$$\begin{equation} \psi = f_*\left(\frac{M_\mathrm{cool}}{t_\mathrm{dyn}}\right), \end{equation} \tag{ 6 }$$

where t_dyn is the halo dynamical time and f_* is a free parameter. The cold gas mass is calculated by defining a local cooling rate within the halo,

$$\begin{equation} t_\mathrm{cool}(r)=\frac{3k_BT\rho _g(r)}{2\mu n_{{\rm H}}^2(r)\Lambda (T)}, \end{equation} \tag{ 7 }$$

where ρ_g is the (spherically symmetric) gas density profile, n_H is the hydrogen number density, μ is the molecular weight, and Λ is the cooling function. A cooling radius r_cool can then be defined by

$$\begin{equation} t_\mathrm{cool}(r_\mathrm{cool})=t_\mathrm{dyn}, \end{equation} \tag{ 8 }$$

which gives the cooling rate as

$$\begin{equation} \frac{{dM}_\mathrm{cool}}{dt} = 4\pi \rho _g(r_\mathrm{cool})r^2_\mathrm{cool}\frac{dr_\mathrm{cool}}{dt}. \end{equation} \tag{ 9 }$$

Initially all gas entering the halo is shock-heated to the virial temperature. Integrating Equation (9) then gives the amount of cool gas at any time. We use a gas profile close to isothermal, which fits the results of Springel & Hernquist (2003) to better than 10%. Equation (6) follows from the empirical Kennicutt–Schmidt (KS) relation (Krumholz & Thompson 2007) in the limit of marginally unstable disks. The KS relation does not exhibit any evolution till z ∼ 2 (Daddi et al. 2010).

The galaxy dynamical time can be written as

$$\begin{equation} t_\mathrm{dyn} = 2\times 10^7 \,\mathrm{yr} \left(\frac{R_{1/2}}{4\,\mathrm{kpc}}\right)\left(\frac{V_c}{200\,\mathrm{km}\,\mathrm{s}^{-1}}\right)^{-1}, \end{equation} \tag{ 10 }$$

where R_1/2 is galaxy disk scale length and V_c is the halo circular velocity (Bouché et al. 2010). We take the disk scale length to be proportional to the halo virial radius

$$\begin{equation} R_{1/2} = 0.05 \left(\frac{\lambda }{0.1}\right)R_\mathrm{vir}, \end{equation} \tag{ 11 }$$

where the halo spin parameter λ is taken to be 0.07 (Krumholz & Dekel 2012). The halo circular velocity is given by

$$\begin{eqnarray} V_\mathrm{c}&=&23.4\left(\frac{M}{10^8\,h^{-1}\,M_\odot }\right)^{1/3}\left[\frac{\Omega _m}{\Omega _m^z}\frac{\Delta _c}{18\pi ^2}\right]^{1/6}\nonumber \\ &&\times \left(\frac{1+z}{10}\right)^{1/2} \,\mathrm{km}\,\mathrm{s}^{-1}, \end{eqnarray} \tag{ 12 }$$

where

$$\begin{equation} \Omega _m^z=\frac{\Omega _m(1+z)^3}{\Omega _m(1+z)^3+\Omega _\Lambda +\Omega _k(1+z)^2}, \end{equation} \tag{ 13 }$$

and Δ_c is the halo overdensity relative to the cosmological critical density, given by

$$\begin{equation} \Delta _c=18\pi ^2+82d-39d^2, \end{equation} \tag{ 14 }$$

for the ΛCDM cosmology, where $d=\Omega _m^z-1$ (Barkana & Loeb 2001). The virial radius is given by

$$\begin{eqnarray} R_\mathrm{vir}&=&0.784\left(\frac{M}{10^8 \,h^{-1}\,M_\odot }\right)^{1/3}\left[\frac{\Omega _m}{\Omega _m^z}\frac{\Delta _c}{18\pi ^2}\right]^{-1/3}\nonumber \\ &&\times \left(\frac{1+z}{10}\right)^{-1} \,h^{-1}\,\mathrm{kpc}. \end{eqnarray} \tag{ 15 }$$

The total stellar mass in the halo evolves as

$$\begin{eqnarray} \dot{M}_*(z) &=& {-} \dot{M}_{g, \mathrm{sf}}(z) = \psi (z)-\int _{m_l}^{m_u}dm\,\phi (m)\nonumber \\ &&\cdot \psi [t(z)-\tau (m)]\cdot [m-m_r(m)], \nonumber\\ \end{eqnarray} \tag{ 16 }$$

where, as in Equation (3), the second term on the right-hand side accounts for mass loss from stars.

2.2.3. Metals

The initial mass of metals in a halo is zero. Metals are produced in stars in the halo and are mixed in the halo gas when stars explode as supernovae, and due to stellar winds. The total mass of metals M_Z in the halo is also affected by inflows from the IGM and outflows into it. The inflow term is given by

$$\begin{equation} \dot{M}_{Z,\mathrm{in}}(z) = Z_\mathrm{IGM}(z)\dot{M}_{g,\mathrm{in}}, \end{equation} \tag{ 17 }$$

where Z_IGM is the metal abundance in the IGM, which we describe below.

Metal outflows can be described in similar fashion, as

$$\begin{equation} \dot{M}_{Z,\mathrm{out}}(z) = -Z_g(z)\dot{M}_{g,\mathrm{out}}, \end{equation} \tag{ 18 }$$

where Z_g ≡ M_Z/M_g is the gas metal abundance in the halo.

Lastly, the effect of star formation can be expressed as (Krumholz & Dekel 2012)

$$\begin{eqnarray} \dot{M}_{Z,\mathrm{sf}}(z) &=& \int _{m_l}^{m_u}dm\,\phi (m)\cdot \psi [t(z)-\tau (m)]\nonumber \\ &&\times mp_Z(m)(1-\zeta). \end{eqnarray} \tag{ 19 }$$

Here, p_Z(m) is the mass fraction of a star of initial mass m that is converted to metals and ejected. The factor 1 − ζ takes into account the fact that not all of the newly enriched material is going to mix in the ISM. Some of it will be directly ejected out of the halo due to supernova explosions. We follow Krumholz & Dekel (2012) and write

$$\begin{equation} \zeta = \zeta _l\,\exp {(-M_{h,12}/M_\mathrm{ret})}, \end{equation} \tag{ 20 }$$

where M_{h, 12} is the halo mass in units of 10¹² M_☉ and we set the parameters M_ret = 0.3 and ζ_l = 0.9. These values are expected to depend on the stellar IMF and the problem geometry. Our choice of these values is motivated by the simulations of Mac Low & Ferrara (1999) and the findings of Krumholz & Dekel (2012). These authors find that changing these quantities within reasonable limits has only a modest effect on the star formation and metallicity evolution.

The evolution of the metal mass of a halo can now be written as

$$\begin{equation} \dot{M}_Z(z) = \dot{M}_{Z,\mathrm{sf}}(z) + \dot{M}_{Z,\mathrm{out}}(z) + \dot{M}_{Z,\mathrm{in}}(z). \end{equation} \tag{ 21 }$$

We always take the stellar IMF to have the Salpeter form, which is given by

$$\begin{equation} \phi (m)=\phi _0\,m^{-2.3}, \end{equation} \tag{ 22 }$$

where the constant ϕ₀ is chosen such that

$$\begin{equation} \int _{m_0}^{m_1}dm\,m\,\phi (m)=1\, M_\odot. \end{equation} \tag{ 23 }$$

We take m₀ = 0.1 M_☉ and m₁ = 100 M_☉ for Population I/II stars. Population III stars have different IMF in each of our models, as we want to discuss the effect of Population III IMF on various quantities. These are shown in Table 1. The stellar IMF in model 3 covers the PISN range, whereas models 1 and 2 explore the effect of asymptotic giant branch stars and core-collapse supernovae, respectively. Each model is calibrated separately to reproduce the cosmic SFR history, the IGM Thomson scattering optical depth to the last scattering surface (τ_e = 0.089 ± 0.014; Hinshaw et al. 2012), and the hydrogen photoionization rate in the IGM as measured from Lyα forest observations.

Our model requires stellar lifetimes because we do not assume instantaneous recycling. We take lifetimes for low- and intermediate-mass stars (0.1 M_☉ < M < 100 M_☉) from Maeder & Meynet (1989). Lifetimes for stars with higher mass (all of which are Population III) are taken from Schaerer (2002). The model also requires chemical and UV photon yields of stars with different IMFs. Woosley & Weaver (1995) have calculated the chemical yields of massive stars (12 M_☉ < M < 40 M_☉) with different metallicities (Z = 0, 10⁻⁴, 0.01, 0.1, and 1 times solar). We use their results by interpolating between different metallicities and extrapolating beyond the mass range for lower and higher stellar masses. Chemical yields of Population III stars (which are metal-free and high mass) are taken from the calculations of Heger & Woosley (2002). Stellar spectra of Population I/II stars are calculated using starburst99 (Leitherer et al. 1999; Vázquez & Leitherer 2005) with respective metallicities. Synthetic spectra of Population III stars are taken from Schaerer (2002).

The transition from Population III star formation to Population II star formation in any halo is implemented via a critical metallicity, Z_crit, with our fiducial value as Z_crit = 10⁻⁴ (Bromm et al. 2001; Bromm & Loeb 2003; Frebel et al. 2007). When the ISM metallicity in a halo crosses Z_crit, new stars are formed according to a Population II IMF. We consider the effect of changing the value of Z_crit below.

Note that in the chemical evolution model described above, we have only accounted for core-collapse supernovae; the contribution of Type Ia supernovae (SNe Ia) has been ignored. As implicitly assumed in Equation (19), the core-collapse supernova rate traces star formation activity. The SN Ia rate has a more complex dependence on the stellar IMF and SFR.

Classically, progenitors of SNe Ia are thought to be intermediate-mass stars. An unknown delay is expected between the death of the progenitor star and the SN Ia explosion. By using a similar model of cosmic star formation as that presented here, Daigne et al. (2006) argued that the observed rate of SN Ia occurrence suggests that the typical delay time is long (∼3–3.5 Gyr). Using a different model, Gal-Yam & Maoz (2004) reached a similar conclusion, arguing that observations favor long delay times (∼1 Gyr). This is also in agreement with recent DLA enrichment measurements at z = 4–5 by Rafelski et al. (2012), who report α-enhancement of about 0.3 in these objects, suggesting that these objects have not yet been enriched by SNe Ia. (As we will discuss below in Section 4, the effect of dust depletion is expected to be low in these systems due to their low metallicity.) If the SN Ia delay times are indeed sufficiently long, the results of this paper will not be affected, as we are interested in DLA abundances at the highest redshifts (z ≳ 4).

However, Mannucci et al. (2006) have argued that the observed dependence of SN Ia rate on the colors of parent galaxies and observed SN Ia rates in radio-loud early-type galaxies favor a bimodal delay time distribution, in which about half of SNe Ia explode soon after the birth of their progenitor star (short delay time of ∼100 Myr), while the remaining half have long delay time of the order of 3 Gyr. In the context of DLAs, Calura et al. (2003) have also shown that when dust depletion is taken into account, most of the claimed DLA α-enhancements vanish, which suggests a role of SNe Ia (cf. Malaney & Chaboyer 1996). While dust is not expected to play a major role in z ∼ 5 DLAs, these prompt SNe Ia with short delay times can have an effect on the results of our model. Nevertheless, in this paper, we focus on the effect of Population III star formation on the abundance ratios in this paper. We will consider the effects of dust depletion and prompt SNe Ia in a future work (G. Kulkarni et al., in preparation). We conjecture here, however, that while such prompt SNe Ia will affect our results, a strong degeneracy with the effect of Population III IMF is unlikely as the latter leave a sufficiently large signature in ratios between α-element abundances (∼0.5–1 dex at z ∼ 6).

Cosmological inflows bring material from the IGM into the ISM of galaxies. The extent to which this dilutes the metal content of the ISM depends on the metallicity of the IGM. We follow the model of Daigne et al. (2006) to calculate the average IGM metallicity at any redshift. This model divides cosmic baryons into three reservoirs: (1) intergalactic medium (subscript "IGM"), (2) interstellar medium ("ISM"), and (3) stars ("str"). We denote the mass densities of these as M_IGM, M_ISM, and M_str, respectively. Note that M_ISM and M_str are mass function weighted integrals of M_g and M_* defined for each halo mass bin above.

The IGM mass density evolves according to

$$\begin{equation} \frac{{dM}_\mathrm{IGM}}{dt}=-a_b(t)+o(t), \end{equation} \tag{ 24 }$$

where a_b is the total rate of accretion of baryons from the IGM onto halos of all masses, and o(t) is the total rate of outflow of baryons from halos into the IGM. The stellar mass density evolves as

$$\begin{equation} \frac{{dM}_\mathrm{str}}{dt}=\Psi (t)-e(t), \end{equation} \tag{ 25 }$$

where Ψ is the total cosmic SFR and e is rate of ejection of material of stars into the ISM via winds and supernovae, which depends on the stellar IMF. The total cosmic SFR is given by

$$\begin{equation} \Psi (t)=\int ^{M_h^\mathrm{high}}_{M_h^\mathrm{low}} {dM} \dot{M}_*(M,t) N(M,t), \end{equation} \tag{ 26 }$$

where $M_h^\mathrm{high}$ and $M_h^\mathrm{low}$ are the halo mass limits considered in our calculation (as discussed above), $\dot{M}_*(M,t)$ is SFR within a halo of mass M at cosmic time t (given by Equation (16)), and N(M, t) is the halo mass function. Lastly, the ISM mass density evolves as

$$\begin{equation} \frac{{dM}_\mathrm{ISM}}{dt} = -\frac{{dM}_\mathrm{IGM}}{dt} - \frac{{dM}_\mathrm{str}}{dt}. \end{equation} \tag{ 27 }$$

To calculate the IGM metallicity, in addition to the total masses, the mass fraction of individual elements in each reservoir is also evolved. For an element i, the mass fraction in the IGM is $X^\mathrm{IGM}_i=M^\mathrm{IGM}_i/M_\mathrm{IGM}$ and that in the ISM is $X^\mathrm{ISM}_i=M^\mathrm{ISM}_i/M_\mathrm{ISM}$. It can then be shown that mass fractions in the IGM evolve as

$$\begin{equation} \frac{{dX}^\mathrm{IGM}_i}{dt}=o\left(X_i^\mathrm{ISM}-X_i^\mathrm{IGM}\right)/M_\mathrm{IGM}, \end{equation} \tag{ 28 }$$

and those in the ISM evolve as

$$\begin{eqnarray} \frac{{dX}^\mathrm{ISM}_i}{dt}&=&\left[a_b\left(X_i^\mathrm{IGM}-X_i^\mathrm{ISM}\right)+\left(e_i-eX_i^\mathrm{ISM}\right)\right]/M_\mathrm{ISM}.\nonumber \\ \end{eqnarray} \tag{ 29 }$$

Here, o is the outflow rate from the ISM and e_i is the rate at which metal i is ejected by stars. In order to completely specify the model, all that is left to specify is (1) initial conditions for above equations; and (2) method of calculation for the three quantities o, a_b, and e.

The total baryon accretion rate is given by

$$\begin{eqnarray} a_b(t) &=& \Omega _b\left(\frac{3 {\rm H}_0^2}{8\pi G}\right)\left(\frac{dt}{dz}\right)^{-1}\left|\frac{df_{\mathrm{struct}}}{dz}\right| \nonumber \\ &=& 1.2\,h^3\,M_\odot \,\mathrm{yr}^{-1}\,\mathrm{Mpc}^{-3}\left(\frac{\Omega _b}{0.044}\right)\nonumber \\ & &\times (1+z)\sqrt{\Omega _\Lambda +\Omega _m(1+z)^3}\left|\frac{df_{\mathrm{struct}}}{dz}\right|, \end{eqnarray} \tag{ 30 }$$

where f_struct is the fraction of total mass in star-forming halos. It is given by

$$\begin{equation} f_{\mathrm{struct}}(z)=\frac{\int ^\infty _{M_\mathrm{min}} {dM} M f_\mathrm{PS}(M,z)}{\int ^\infty _0 {dM} M f_\mathrm{PS}(M,z)}, \end{equation} \tag{ 31 }$$

and f_PS is the Press–Schechter mass function.

The outflows from the ISM into the IGM are described using the following relation:

$$\begin{equation} o(t)=\frac{2\epsilon }{v^2_\mathrm{esc}(z)}\int ^{m_u}_{m_l} dm\,\phi (m)\cdot \Psi [t(z)-\tau (m)]\cdot E_\mathrm{kin}(m), \end{equation} \tag{ 32 }$$

where, as before, ϕ is the stellar IMF, and E_kin is the kinetic energy released by the explosion of a star of mass m. The quantity v_esc is the "escape velocity," which is calculated according to

$$\begin{equation} v^2_\mathrm{esc}(z)=\frac{\int ^\infty _{M_\mathrm{min}}{dM} f_\mathrm{PS}(M,z) M (2GM/R_\mathrm{vir})}{\int ^\infty _{M_\mathrm{min}}{dM} f_\mathrm{PS}(M,z) M}. \end{equation} \tag{ 33 }$$

The ejection of enriched gas from stars into the ISM is given by

$$\begin{equation} e(t)=\int _{m_l}^{m_u}dm\,\phi (m)\cdot \Psi [t(z)-\tau (m)]\cdot [m-m_r(m)], \end{equation} \tag{ 34 }$$

as above, where m_u is the upper limit on mass of stars that explode and produce supernova remnants, and m_r is the remnant mass.

Our model for reionization and thermal history of the average IGM is essentially that developed by Choudhury & Ferrara (2005), with the main difference being in implementation of sources. In this model, which matches a wide range of observational constraints, the IGM is gradually reionized and reheated by star-forming galaxies between z ∼ 20 and z ∼ 6. We summarize the main elements of the model here and refer the reader to that paper for details.

The model accounts for IGM inhomogeneities by adopting a lognormal density distribution with the evolution of volume filling factor of ionized hydrogen (H ii) regions, $Q_{{\rm H\,\scriptsize{II}}}(z)$ being calculated according to the method outlined by Miralda-Escudé et al. (2000). The volume filling factor of H ii regions evolves as

$$\begin{equation} \frac{d[Q_{{\rm H\,\scriptsize{II}}}F_M(\Delta)]}{dt} = \frac{\dot{n}_\nu }{n_{{\rm H}}} - Q_{{\rm H\,\scriptsize{II}}}\alpha _R(T)n_eR(\Delta)(1+z)^3, \end{equation} \tag{ 35 }$$

where F_M(Δ) is the mass fraction of the IGM occupied by regions with density less than Δ, $\dot{n}_\nu$ is the rate at which ionizing photons are introduced in the IGM by galaxies, n_e is the electron number density, α_R(T) is the temperature-dependent recombination rate, and R(Δ) is the clumping factor of the IGM, which is related to the second moment of the IGM density distribution. The IGM mass fraction F_M(Δ) is given by

$$\begin{equation} F_M(\Delta _\mathrm{crit}) = \int _0^{\Delta _\mathrm{crit}} d\Delta \Delta P(\Delta), \end{equation} \tag{ 36 }$$

and the clumping factor is given by

$$\begin{equation} R(\Delta _\mathrm{crit}) = \int _0^{\Delta _\mathrm{crit}} d\Delta \Delta ^2 P(\Delta). \end{equation} \tag{ 37 }$$

Reionization is said to be complete once all low-density regions (say, with overdensities Δ < Δ_crit ∼ 60) are ionized. We refer to the redshift at which this happens as the redshift of reionization, z_reion. We follow the ionization and thermal histories of neutral and H ii regions simultaneously and self-consistently, treating the IGM as a multi-phase medium. We assume that helium is singly ionized together with hydrogen and contributes to the photoionization heating. We do not consider the double ionization of helium, which is thought to have occurred at z ∼ 3 (e.g., Worseck et al. 2011).

The rate of ionizing photons injected by the galaxies into the IGM per unit time per unit volume at redshift z is denoted by $\dot{n}_\mathrm{ph}(z)$, which is determined by the SFR. Using our star formation model described above, we can write the rate of emission of ionizing photons per unit time per unit volume per unit frequency range, $\dot{n}_\nu (z)$, as

$$\begin{equation} \dot{n}_\nu (z)=f_\mathrm{esc}\left[N^\mathrm{II}_\gamma (\nu)\Psi ^\mathrm{II}(z)+N^\mathrm{III}_\gamma (\nu)\Psi ^\mathrm{III}(z)\right], \end{equation} \tag{ 38 }$$

where f_esc is the escape fraction of UV photons, and $N^\mathrm{II}_\gamma (\nu)$ and $N^\mathrm{III}_\gamma (\nu)$ are the total number of ionizing photons emitted per unit frequency range per unit stellar mass from Population I/II and Population III star-forming halos, respectively. These quantities are calculated by integrating the stellar spectra over an appropriate range. Ψ^II(z) and Ψ^III(z) are the cosmic Population I/II and Population III SFRs. At any redshift,

$$\begin{equation} \Psi (z) = \Psi ^\mathrm{II}(z) + \Psi ^\mathrm{III}(z). \end{equation} \tag{ 39 }$$

The total rate of emission of ionizing photons per unit time per unit volume, $\dot{n}_\mathrm{ph}(z)$, is obtained by simply integrating Equation (38) over the correct frequency range. The hydrogen photoionization rate is then given by

$$\begin{equation} \Gamma _{\mathrm{H\,\scriptsize{I}}}(z)=(1+z)^3\int _{\nu _0}^\infty d\nu \lambda (z,\nu) \dot{n}_\nu (z) \sigma (\nu), \end{equation} \tag{ 40 }$$

where σ(ν) is the hydrogen photoionization cross-section and λ(z, ν) is the redshift-dependent mean free path of UV photons. The mean free path is governed by the distribution of high-density neutral clumps and reproduces the observed density of Lyman limit systems at low redshift. The photoionization rate is used with the average radiative transfer equation to calculate the evolution of $Q_{{\rm H\,\scriptsize{II}}}(z)$ in the pre-reionization universe and that of Δ_crit(z) in the post-reionization universe. Note that this approach assumes ionization equilibrium and local absorption of UV photons.

The mean free path is defined as the average distance between high-density regions of the inhomogeneous IGM. It can be written as (Choudhury & Ferrara 2005)

$$\begin{equation} \lambda (z) = \frac{\lambda _0}{[1-F_V(\Delta _\mathrm{crit})]^{2/3}}, \end{equation} \tag{ 41 }$$

where F_V(Δ_crit) is the volume fraction of the IGM occupied by regions with density less than Δ_crit, given by

$$\begin{equation} F_V(\Delta _\mathrm{crit}) = \int _0^{\Delta _\mathrm{crit}} d\Delta P(\Delta), \end{equation} \tag{ 42 }$$

where P(Δ) is the probability distribution function of IGM overdensities, taken to have a lognormal form (Miralda-Escudé et al. 2000; Choudhury & Ferrara 2005). We take λ₀ to be proportional to the Jeans length (Equation (47)),

$$\begin{equation} \lambda _0=\lambda _{\mathrm{mfp,0}}r_J, \end{equation} \tag{ 43 }$$

and constrain the constant of proportionality, λ_{mfp, 0}, using the number of Lyman limit systems at low redshift, which is related to the mean free path by

$$\begin{equation} \frac{dN_\mathrm{LLS}}{dz}=\frac{c}{\sqrt{(}\pi)\lambda (z)H(z)(1+z)}. \end{equation} \tag{ 44 }$$

We assume that the mean free path is independent of the photon energy up to the ionization threshold of He ii and is proportional to $(\nu /\nu _{\mathrm{He\,\scriptsize{II}}})^{1.5}$ at higher energies (Choudhury & Ferrara 2005).

The IGM temperature evolves as

$$\begin{eqnarray} \frac{dT}{dt}&=&{-}2H(z)T-\frac{T}{1+X_{\mathrm{H\,\scriptsize{I}}}}\frac{{dX}_{\mathrm{H\,\scriptsize{I}}}}{dt}\nonumber \\ &&+\frac{2}{3k_Bn_b(1+z)^3}\frac{dE}{dt}, \end{eqnarray} \tag{ 45 }$$

where dE/dt is the total photoheating rate and $X_{\mathrm{H\,\scriptsize{I}}}$ is the H i fraction. The photoheating term accounts for heating by the UV background and cooling due to recombinations and Compton scattering of cosmic microwave background photons. The IGM temperature governs the thermal feedback on galaxy formation by stopping gas inflow on halos with mass less than M_min(z). We take this to be the Jeans mass, given by

$$\begin{equation} M_\mathrm{min}=\frac{4}{3}\pi r_J^3, \end{equation} \tag{ 46 }$$

where r_J is the comoving Jeans length

$$\begin{equation} r_J = \left[\frac{2k_BT\gamma }{8\pi G \rho _c m_p \Omega _m (1+z)}\right]^{1/2}, \end{equation} \tag{ 47 }$$

where γ is the adiabatic index and m_p is the mass of the proton. Note that before reionization is complete, the H ii regions are distinct. During this epoch, the temperature in H ii regions is different from that in H i regions. The average temperature is then given by the mean of these two temperatures, weighted by $Q_{{\rm H\,\scriptsize{II}}}(z)$. Similarly, both regions have different M_min. This is relevant in understanding the role that M_min plays in our star formation prescription above.

Note that in all of the above, we evolve two distinct sets of halos corresponding to halos forming in H i and H ii regions, respectively. Each set of halos is affected by photoionization feedback differently as the temperatures of H i and H ii regions are different. Any average property, such as the cosmic SFR, at a given redshift is an average of these two sets weighted by the volume filling fraction of H i and H ii regions. This distinction is unimportant at redshift z < z_reion, when there are no H i regions anymore.

Given the above model, we obtain best-fit parameters by comparing with various observations, such as the redshift evolution of photoionization rate obtained from the Lyα forest (Bolton & Haehnelt 2007), the electron scattering optical depth (Hinshaw et al. 2012), and the total baryon fraction in structures at z = 0 (Fukugita & Peebles 2004). These parameters are summarized in Table 2 along with their fiducial values and the constraints used to obtain them.

Table 2. Free Parameters in Our Model

Parameter	Fiducial Value	Constraint
fesc	0.04	τe and $\Gamma _{\mathrm{H\,\scriptsize{I}}}$
f* (Pop. II)	0.02	τe and $\Gamma _{\mathrm{H\,\scriptsize{I}}}$
f* (Pop. III)	0.04	τe and $\Gamma _{\mathrm{H\,\scriptsize{I}}}$
λmfp, 0	1.7 × 10−3	dNLLS/dz
	1.0 × 10−6	fb, struct(z = 0)
fgas, in	0.7	Bouché et al. (2010)
ζ	Equation (20)	Krumholz & Dekel (2012)
Zcrit	10−4 Z☉	⋅⋅⋅

Download table as:  ASCIITypeset image

Before moving on to the results of our model, we compare our galaxy formation model with conventional semianalytical models of galaxy formation used in the literature (e.g., Kauffmann et al. 1999; Cole et al. 2000; De Lucia et al. 2004; Baugh 2006). The basic approach in these models is two-step: first, halo merger trees are constructed from N-body cosmological simulations; next, baryonic physics is calculated in each individual halo using prescriptions similar to those described above. When halos merge, the stellar, gas, and metal masses are suitably combined, possibly with some extra prescriptions for events, such as starbursts or active galactic nuclei (AGNs). Our approach in modeling galaxy formation in this paper is also two-step. However, instead of drawing halo merger trees from N-body simulations, we only draw the average mass assembly histories of halos. This is a significant reduction in information, as the average mass at redshift z of a halo with mass M at z = 0 is the mean of the masses of all progenitors at redshift z of all such halos in the simulation box (see Fakhouri et al. 2010 for details of how this is done). In particular, by considering only the mean assembly histories, we lose any information regarding the intrinsic scatter in any property (such as halo mass, stellar mass, gas mass, or chemical abundance), even for fixed halo mass. Nonetheless, this simplification affords us the ability to calculate the IGM thermal and ionization histories and the evolution of H i and H ii regions in the IGM by solving the average radiative transfer equation (Equation (35)). It is not possible to do this in conventional merger-tree-based semianalytic models.

Although our approach is different from conventional semianalytic models, we expect our results to correspond to averaged results of those models. Thus, e.g., the stellar mass at a given redshift of a given halo of mass, say, M in our model corresponds to the mean stellar mass of all halos at that redshift in a merger-tree-based model. This can be seen from the agreement between our model predictions and observational data for average quantities, such as cosmic SFR, gas-phase metallicity, and mass–metallicity relation, as we discuss below.

Our model for galaxy and IGM evolution is constrained to produce the observed SFR density evolution, the fraction of total baryon density in halos at low redshift, and the reionization history of the IGM as measured by the Thomson scattering optical depth to the last scattering surface (τ_e) and the hydrogen photoionization rate ($\Gamma _{\mathrm{H\,\scriptsize{I}}}$). Figure 1 shows the evolution of the cosmic SFR density. The solid curves show the prediction of our model, and the red data points show observational measurements from a compilation by Hopkins & Beacom (2006). All three models reproduce the observed rise in the SFR between z = 0 and z ∼ 3. The total SFR predictions of the three models are virtually identical. As we discuss in detail below, this is because the contribution from Population III stars is very small. However, given the error bars on the measurements and the simplicity of our model, we consider the match to be adequate. Note that recent measurements of the SFR density at z ∼ 7–10 by Bouwens et al. (2010) and Ellis et al. (2013) are smaller than our model prediction at these redshifts. These measurements report values of ≲ 10⁻² M_☉ yr⁻¹ Mpc⁻³ at z ∼ 8 and ≲ 10⁻³ M_☉ yr⁻¹ Mpc⁻³ at z ∼ 10. Our values at both these redshifts are ≳ 10⁻² M_☉ yr⁻¹ Mpc⁻³. However, this discrepancy is because the SFR estimates of Bouwens et al. (2010) and Ellis et al. (2013) do not count the large contribution to the SFR by faint galaxies (which are more numerous). These estimates are obtained by integrating the UV luminosity functions down to $0.05L^*_{z=3}$. This limit corresponds to an AB magnitude of about −18 at 1500 Å, while we find galaxies down to ≳ − 17 in our model. In a forthcoming paper (G. Kulkarni et al., in preparation), we analyze this issue in detail and show that our model predictions are in good agreement with high-redshift SFR estimates when the faintness limit of Bouwens et al. (2010) and Ellis et al. (2013) is applied.

Zoom In Zoom Out Reset image size

Figure 2 shows the thermal evolution of the IGM in our models in comparison with some observational measurements (Schaye et al. 2000; Lidz et al. 2010; Becker et al. 2011). The solid curve in the left panel shows the evolution of the temperature of the IGM at mean density (T₀), which is the temperature usually measured from absorption spectra of high-redshift quasars. At redshift z ≳ 20, T₀ decreases adiabatically ∝(1 + z)². However, as the H ii regions around first galaxies form, this temperature gradually rises to about 2 × 10⁴ K at z = z_reion ∼ 6. Once reionization is complete, only the highest density regions in the IGM are left to be ionized. Thus, the mean density regions evolve unaffected and T₀ then decreases to ∼10⁴ K. At redshift higher than z_reion, T₀ is the mean of temperatures in H i and H ii regions, weighted by the volume filling factor. As a result, at these redshifts temperatures in these regions are more relevant. In Figure 2 (left panel), the red dashed curve shows the temperature evolution in H ii regions, and the blue dashed curve shows the same in H i regions. The average H ii region temperature climbs to a high value of a few times 10⁴ K at z ∼ 20 when the first sources turn on. However, this does not affect T₀ at these redshifts as the H ii region volume filling factor is small. The average temperature in H i regions continues to decrease with time due to adiabatic cooling. Note that as we do not model He ii reionization in this paper, our ability in matching the observational measurements is restricted. All observational measurements are made in the redshift range z = 2–4, which is believed to be the epoch of He ii reionization. We also note in passing that, while the only source of heat in our model is photoionizations, the temperature measurements in the redshift range z = 2–4 are not precise enough to rule out other, exotic, sources (Bolton et al. 2004; McQuinn et al. 2009; Calura et al. 2012; Meiksin & Tittley 2012; Rudie et al. 2012; Garzilli et al. 2012).

Zoom In Zoom Out Reset image size

Recall that here the volume filling factor of H ii regions is calculated according to Equation (35). Present observational constraints on this quantity are not very strong. However, several upper limits have been suggested in the literature. These include: (1) a lower limit of $Q_{\mathrm{H\,\scriptsize{II}}}\gtrsim 0.6$ derived by Ouchi et al. (2010) at z = 6.6 from the evolution of the number density of Lyα emitting galaxies, (2) an upper limit on $Q_{\mathrm{H\,\scriptsize{I}}}$ at z = 5.5 and 6 measured by Mesinger (2010) by counting dark pixels in quasar spectra, (3) a lower limit of $Q_{\mathrm{H\,\scriptsize{II}}}\gtrsim 0.5$ inferred by McQuinn et al. (2007) at z = 6.6 from the lack of an increase in the clustering of Lyα emitting galaxies, (4) an upper limit of $Q_{\mathrm{H\,\scriptsize{II}}}\lesssim 0.1$ inferred from the proximity zone of a z = 7.1 quasar by Bolton et al. (2011), and (5) an upper limit of $Q_{\mathrm{H\,\scriptsize{I}}}\lesssim 0.5$ inferred by McQuinn et al. (2008) from the red damping wing of the Lyα absorption line in the spectrum of a z = 6.3 GRB. Our model is consistent with all of these constraints as in our model reionization is complete by redshift z = 8.

The thermal evolution of the IGM dictates the evolution of M_min, the minimum mass of star-forming halos, which is shown in the right panel of Figure 2. This quantity implements photoionization feedback in our model. As described in the previous section, we calculate M_min using a Jeans criterion. At a given redshift, the minimum mass of star-forming halos is either given by the atomic cooling threshold of 10⁴ K, or by the local Jeans mass, whichever is higher. Other models of M_min are present in the literature, which take into account the full thermal history of the IGM instead of the instantaneous temperature (Gnedin 2000; Okamoto et al. 2008). We ignore these improvements and focus on the Jeans method for simplicity; this approximation does not affect our results strongly. The solid curve in Figure 2 (left panel) shows the average M_min evolution. It increases monotonically with decreasing redshift and is ∼2 × 10⁸M_☉ at z = z_reion. It is ∼6 × 10⁸M_☉ at redshift 1. It is important to note that this is only the average M_min; the minimum mass is different in H i and H ii regions. These are shown by the blue and red curves in Figure 2. Thus, galaxies forming in these regions are affected differently. The minimum mass in H ii region grows to about 10⁸ M_☉ as soon as the first sources are turned on and the temperature of these regions is boosted to ∼10⁴ K. On the other hand, the minimum mass in H i regions increases much more slowly with decreasing redshift as it is always determined by the atomic cooling threshold alone. The relative contribution of the two regions to the solid black curve in Figure 2 is determined by the filling factor of H ii regions. Of course, this distinction is not important in the post-reionization phase, where there are no H i regions.

As seen in the thermal evolution in Figure 2, IGM reionization in our model is gradual. This gradual change in the ionization state of the IGM helps us simultaneously reproduce the two most robust observational constraints on reionization: (1) the observed Thomson scattering optical depth to the last scattering surface τ_e = 0.089 ± 0.014 (Hinshaw et al. 2012), and (2) the measurements of the hydrogen photoionization rate, $\Gamma _{\mathrm{H\,\scriptsize{I}}}$, from observations of the Lyα forest (Meiksin & White 2004; Bolton & Haehnelt 2007; Faucher-Giguère et al. 2008). The hydrogen photoionization rate in our model 3 is shown in Figure 3; our other models (1 and 2) give very similar results. The model prediction matches very well with the observational measurements. The photoionization rate increases rapidly as UV photon sources build up at z > z_reion. There is a sudden jump at z = z_reion when different H ii regions overlap. This is because at this redshift H ii regions overlap and a given point in the IGM starts "seeing" multiple sources, which rapidly enhances the UV photon mean free path, thereby affecting the photoionization rate.

Zoom In Zoom Out Reset image size

We now turn to the mass–metallicity relation prediction in our models, which is shown in Figure 4. The mass–metallicity relation can be easily predicted in our model because for each halo mass, we calculate the stellar mass and gas-phase abundances of various metals at various redshifts, as explained above in Section 2. In particular, in Figure 4, we compare the stellar mass with the oxygen abundance. We find that the model predictions match quite well with observations at z = 2.27 (Erb et al. 2006) and z = 3.7 (Maiolino et al. 2008). The oxygen abundance increases with stellar mass for galaxies in the intermediate-mass range. At the high-mass end (M_* > 10¹¹ M_☉), we find a slow decline in the oxygen abundance. This is because the gas accretion rate in these halos is high, which results in dilution of their ISMs. Also, it is difficult for outflows to remove gas from the deeper potential wells of these halos. Note that the measurements of the mass–metallicity relationship are quite uncertain (Kewley & Ellison 2008; Krumholz & Dekel 2012). Therefore, to compare our estimate of metallicity to the observations, we force the zero point of the predictions to match with that of the observations at z = 2.27. The model mass–metallicity relationship has a slightly different slope than the observational points at z = 2.27. However, given the simplicity of our model, and our focus on ratios of abundances in this work, we take the level of agreement to be sufficient. We will also see below that when the average ISM metallicity of galaxies is averaged it agrees perfectly with the mean ISM metallicity inferred from observations of DLAs.

Zoom In Zoom Out Reset image size

In summary, our average galaxy formation and chemical evolution model is consistent with several observations, including the observed Thomson scattering optical depth to the last scattering surface, the observed hydrogen photoionization rate inferred from the Lyα forest, the slope of the mass–metallicity relation of galaxies at various redshifts, and the evolution of cosmic SFR density. It also agrees well with the low-mass slope, characteristic mass, and normalization of the stellar-to-halo mass relation, and, as we will show below, the evolution of the mean ISM metallicity inferred from observations of DLAs.

To study how the Population III stellar IMF affects above results, we now consider three different Population III stellar IMFs shown in Table 1. Several early studies predicted that Population III stars have a characteristic mass of a few hundred M_☉ (e.g., Abel et al. 2002). In recent years, this prediction has come down in the range of 20–100 M_☉ (e.g., Turk et al. 2009). Recently, Dopcke et al. (2013) argued that in the absence of metal-line cooling, dynamical effects can still lead to fragmentation in protostellar gas clouds. In their simulations, this resulted in Population III stars with masses as low as 0.1 M_☉. (This picture predicts the existence of many Population III stars surviving in the Galaxy till the present day.) Enrichment beyond Z_crit merely results in a change in the slope of the IMF, not in its mass values. Given this lack of certainty in our knowledge of the Population III IMF, in this work we choose to work with three different IMFs: 1–100 M_☉ Salpeter, 35–100 M_☉ Salpeter, and 100–260 M_☉ Salpeter. In each case, our Population II IMF is fixed to 0.1–100 M_☉ Salpeter. Each of these three models is calibrated independently, as described above. Our goal is to understand how the difference in these Population III scenarios is reflected in chemical evolution of galaxies. The 1–100 M_☉ and 100–260 M_☉ IMFs are selected to represent two extreme possibilities: if metal-free gas can fragment due to dynamical effects, the IMF will favor small stellar masses, but in the absence of fragmentation Population III stars could have masses of the order of ∼100 M_☉. The intermediate (35–100 M_☉) IMF is chosen to consider the effect of core-collapse supernovae. If the distinct chemical signatures of these three IMFs turn out to be observable in high-redshift galaxies, then such observations could act as a probe of the Population III IMF.

The contribution of Population III stars to the total cosmic SFR density is shown in Figure 1, in which the bottom panel shows the evolution of the Population III SFR (Ψ^III(z)) and the top panel shows the fractional contribution to the total SFR (Ψ^III(z)/Ψ(z)). Evidently, Population III stars contribute very little in terms of SFR. Furthermore, any Population III contribution is essentially zero by redshift z ∼ 7. The initial burst of Population III stars is sufficient to enhance the ISM metallicity beyond Z_crit in any halo mass bin. When this happens, the corresponding galaxy stops forming Population III stars. Since most of the mass in the universe is contained in M_* halos, the globally averaged Population III SFR starts declining as soon as these halos cross the Z_crit threshold. Since the IMF in model 1 has a lower metal yield than that in model 3, halos take longer to cross the Z_crit threshold. Therefore, the Population III SFR contribution is larger in model 1. These results are in good agreement with the hydrodynamic simulations of Wise et al. (2012), who use only very massive Population III stars. We find that while the amplitude of the Population III SFR is sensitive to this effect, the lowest redshift with nonzero cosmic Population III SFR density is more dependent on the evolution of the minimum mass of star-forming halos, M_min. Further, the small contribution of Population III stars to the global SFR is reflected in the fact that the effect of changing the Population III stellar IMF on the reionization and thermal history of the IGM is negligible in our model. As a result, the temperature evolution, as well as the photoionization rate evolution, is practically independent of the Population III IMF.

However, we find that although the effect of Population III on the SFR and the IGM reionization and thermal history is small, it has a significant effect on the chemical properties of low-mass galaxies. This is seen in Figure 5 (left column), which shows abundance ratios of various metal species as a function of halo mass at z = 6 for our three models. We focus on [C/Fe], [O/Si], [Zn/Fe], and [N/O]. It is seen that in halos with mass greater than a few times 10⁹ M_☉, all four abundance ratios are constant, and equal to the Population II values. Changing the Population III IMF has no effect on the abundance ratios in these halos. This is because in these halos, Population III star formation has ceased at a much earlier time and the subsequent Population II star formation has wiped out any chemical signature of Population III star formation. The constancy of the abundance ratios is due to the fact that for a time-independent IMF, abundance ratios are equal to the ratio of corresponding chemical yields (Tinsley 1980; Pagel 2009). In low-mass halos, on the other hand, chemical signatures of Population III star formation have not yet been wiped out. As a result, abundance ratios in these halos depend on the Population III IMF and are different in each of our Population III models. Thus, the chemical signature of Population III stars slowly emerges as we move toward low halo masses. Our model suggests that this is probably the best hope of constraining the Population III SFR density on cosmic timescale and the Population III IMF.

Zoom In Zoom Out Reset image size

When we change the Population III model, the abundance ratios in low-mass halos change significantly. For example, the [O/Si] ratio shows a change of more than 1 dex when we change from model 2 to model 3. This can be seen in Table 3, which shows the Population III yields of O and Si in all three models. Here, the yield of species i is defined as

$$\begin{equation} y_i=\frac{\int dm\cdot \phi (m)\cdot m\cdot p_i(m)}{\int dm\cdot m\cdot \phi (m)}, \end{equation} \tag{ 48 }$$

where p_i(m) is the fraction of initial stellar mass m that is converted to species i. (Table 3 shows log (y_i).) The ratio of oxygen yield to silicon yield is comparable to solar in model 1.³ However, it is much lower than solar in model 3, and much higher than solar in model 2. This is because in model 3 high-mass Population III stars produce large amounts of silicon during the O-burning process. This is reflected in the [O/Si] ratio shown in Figure 5. Similar considerations explain the differences seen in other abundance ratios in low-mass halos in the three models. As we move toward lower redshift, we would expect to see the abundance-ratio values to move away from their Population III values toward the Population II values even in small-mass halos. This is seen in Figure 6 (left column), which shows the same abundance ratios as Figure 5 at z = 2. The abundance-ratio values in low-mass halos are much closer to the Population II values.

Zoom In Zoom Out Reset image size

As described above, our model predicts that Population III stars will leave an imprint on the ratio of abundances of metals in low-mass halos, while any such imprint is wiped out in high-mass halos by subsequent Population II star formation. In this section, we discuss the reason behind this contrast between low-mass and high-mass halos. In a nutshell, it is caused by a difference in Population III-to-II transition redshifts of these halos, which in turn is a result of photoionization feedback.

The solid curves in Figure 7 (top panel) show the average growth of three different halos from H ii regions in one of our models (model 3). These halos have masses of 10^9.5M_☉, 10^10.5 M_☉, and 10^11.5 M_☉, respectively, at z = 1. Their growth is described by Equation (1). The dashed line in this panel shows the evolution of the filtering mass M_min in H ii regions in this model. It can be seen that the three halos satisfy the star formation criterion M > M_min at three different redshifts, approximately z = 10, 20, and 25, respectively. Thus, each halo forms its first stars, which are Population III, at very different times. This has an effect on the chemical evolution of these three mass bins, which is shown in the bottom panel of Figure 7. In this panel, the three solid curves show the evolution of [O/Si] in the same halo mass bins as the top panel. (We still focus on model 3.) Since the Population III IMF in model 3 (100–260 M_☉) produces negative values of [O/Si], we expect a negative value of [O/Si] in each halo mass bin when Population III stars form. However, as time progresses, Population II star formation takes over in each mass bin, because of which the value of [O/Si] approaches +0.2. This is seen to happen in all three halo mass bins. However, this happens at different redshifts for different halo masses because they have different redshifts of Population III star formation (and therefore also the redshifts of Population III-to-II transition). The redshift at which Population III star formation first occurs is the redshift at which the condition M > M_min is first met, and is thus governed by the thermal feedback. As a result, at a given redshift, low-mass halos retain memory of their Population III enrichment history. This is highlighted by the vertical line marking z = 6 on the bottom panel of Figure 7. It can be seen that at this redshift, the two higher halo mass bins have already converged to [O/Si] = 0.2, but the lower mass bin still has [O/Si] = −0.3. It is this effect that causes the trends seen in Figure 6.

Zoom In Zoom Out Reset image size

We now see why changing the Population III stellar IMF will selectively affect low-mass halos. Figure 8 shows the evolution of [O/Si] for the same three halos as Figure 7 in each of our three Population III models. As expected from the yields in Table 3, [O/Si] for Population III stars in models 1 and 2 is greater than 0.2, but is less than 0.0 for model 3. The curves for models 1 and 2 show an enhancement in [O/Si] during the first star formation episode in each halo mass bin. But this occurs at different redshifts for different halo masses, and explains the variation in the trends seen in Figure 6. An alternate way of understanding this is that in a closed system with a non-evolving IMF, all abundance ratios get locked at the value given by the ratio of their respective yields (Tinsley 1980; Pagel 2009)

$$\begin{equation} \frac{Z_\mathrm{O}}{Z_\mathrm{Si}}=\frac{y_\mathrm{O}}{y_\mathrm{Si}}, \end{equation} \tag{ 49 }$$

where the yields are defined in Equation (48). While this is strictly valid only under the IRA, it holds for real galaxies to a very good approximation (Tinsley 1980). As a result, abundance ratios will change only when there is a change in the IMF. This is exactly what happens when the stellar IMF in a galaxy changes from Population III to Population II. For low-mass galaxies, this happens at a lower redshift. Therefore, they retain the Population III signatures at these redshifts. Note that the gas-phase metallicity in these galaxies is greater than Z_crit, so they have already stopped forming Population III stars. Still, their relative abundances have not yet settled on the Population II values, and therefore carry a dependence on the Population III IMF.

It is often mentioned in the literature that an "odd–even" pattern in the abundances is a signature of Population III stars. This refers to the fact that in metal-free stars, the production of species with odd nuclear charge is preferentially suppressed relative to their solar values. As we have seen above, the low-mass halos in our models retain memory of Population III star formation and therefore also show the odd–even pattern in elemental abundances. However, it is important to understand that this odd–even pattern is not unique to Population III and also occurs in metal-poor Population II stars (Cooke et al. 2011). Hence, the crucial point here is that one really wants to search for relative abundance trends which differ from the Population II predictions, and this information is encoded in low-mass halos.

Zoom In Zoom Out Reset image size

In all of the above, the critical metallicity Z_crit, at which a Population III IMF changes over to a Population II IMF, plays a crucial role. However, the value of Z_crit is currently a topic of debate (Klessen et al. 2012; Dopcke et al. 2013). When metallicity Z is greater than Z_crit, the availability of metal-line cooling and/or dust cooling is expected to lead to fragmentation in the ISM that could help form small-mass stars with a Population II type IMF. In this paper, we have assumed Z_crit = 10⁻⁴ Z_☉, a value motivated by most studies of fragmentation in metal-poor gas (Bromm et al. 2001; Bromm & Loeb 2003; Frebel et al. 2007). However, some studies have argued for a much smaller value of Z_crit = 10⁻⁶ Z_☉ (Klessen et al. 2012). Since this number is not very well constrained, we comment on the effect on our results of variation in its value.

We find that our results remain unchanged for any value of Z_crit below 10⁻⁴ Z_☉. This is because in all three of our models, the metallicity of ISM in all halos crosses Z_crit in a very short time after the first burst of Population III star formation. We explicitly checked this in our calculation by changing Z_crit to 10⁻⁵ Z_☉ and 10⁻⁶ Z_☉. In our models 2 and 3, in which the Population III metal yield is higher, there is no effect of a change in the critical metallicity. In model 1, lowering the critical metallicity to 10⁻⁶ Z_☉ has an effect of decreasing the average Population III SFR density, but only for z > 20. When Z_crit is lowered, halos in this model take shorter time to stop forming Population III stars. These results are in qualitative agreement with those of Salvadori et al. (2007), who also found that the influence of Z_crit depends on the Population III stellar IMF. Given the negligible effect of changing Z_crit, we keep Z_crit = 10⁻⁴ Z_☉ in this paper.

As we saw above, low-mass galaxies (halo masses ∼10⁹ M_☉) at high redshift (z ∼ 6) are likely to be useful in constraining Population III IMF and SFR. However, these galaxies are difficult to observe. Even otherwise, measurements of their relative abundances are unlikely to be accurate, as measurements from line emissions involve highly uncertain modeling of the observed H ii regions. As a result, we need to look for an alternative method of observing these galaxies and measuring their metal abundances. We now argue that DLAs in high-redshift quasar spectra provide just such an alternative.

DLAs represent H i reservoirs at high redshift. Metal abundances of DLAs have been measured up to z = 6. Using these data, an "age–metallicity relation" can be constructed for DLAs. Moreover, chemical abundance measurements in DLAs are very accurate. This is because in the cold H i reservoirs of DLAs, most metal species are in their singly ionized state, which can be probed in the optical from ground-based telescopes. Errors in DLA H i column density measurements are typically quite low, around 0.05 dex. Therefore, errors in the corresponding metal abundance measurements ([M/H]) are also low, around 0.1 dex. Thus, observations of DLAs can be used to measure gas-phase metallicities at large cosmological look-back times with high precision (Wolfe et al. 2005). Furthermore, in DLAs, relative abundances can still be measured accurately deep into the reionization epoch (z > 6) using metal-line transitions redward of Lyα, even though Gunn–Peterson absorption precludes measurement of neutral hydrogen (cf. Becker et al. 2012).

In order to predict the properties of DLAs in our model, we assign a mass-dependent, physical, neutral hydrogen cross-section ("size") to each halo, at every redshift. This assignment is performed using the fitting function

$$\begin{equation} \Sigma (M)=\Sigma _0\left(\frac{M}{M_0}\right)^2\left(1+\frac{M}{M_0}\right)^{\alpha -2}, \end{equation} \tag{ 50 }$$

where the constants take the values of α = 0.2, M₀ = 10^9.5M_☉, and Σ₀ = 40 kpc² at z = 3 (Pontzen et al. 2008; Font-Ribera et al. 2012). Values at other redshifts are calculated by mapping halos at these redshifts to halos z = 3 according to circular velocity (Font-Ribera et al. 2012). Our choice of the fitting function in Equation (50) is inspired by Pontzen et al. (2008), who find that it provides a good fit to absorption systems in their hydrodynamical simulation at z = 3. As another check on Equation (50), we look at the average evolution of ISM metallicities. The average gas-phase metallicity at a given redshift is

$$\begin{equation} Z(z) = \frac{\int dm\cdot N(m,z)\cdot Z_\mathrm{ISM}(m,z)\cdot \Sigma (m,z)}{\int dm\cdot N(m,z)\cdot \Sigma (m,z)}. \end{equation} \tag{ 51 }$$

Figure 9 shows this quantity in comparison with metallicity estimates for 100 DLAs with z ≲ 4 by Prochaska et al. (2003). Our models are in good agreement with the mean metallicities of the observed DLAs. It is encouraging that our results are also comparable to those of the recent hydrodynamical simulations (Cen 2012; Fumagalli et al. 2011; Davé et al. 2011).

Zoom In Zoom Out Reset image size

For any measurable property p (such as metallicity, or abundance ratio) of DLAs, we can calculate the number of systems with different values of p in a sample of DLAs. This is called the line density distribution, and with Equation (50) in hand, it can be written as (e.g., Wolfe et al. 2005)

$$\begin{equation} \frac{d^2N}{dXdp}=N(M)\cdot \Sigma (M)\cdot \frac{dl}{{dX}}\frac{{dM}}{{dp}}\cdot (1+z)^3. \end{equation} \tag{ 52 }$$

Here, X is an absorption length element given by

$$\begin{equation} \frac{dl}{{dX}}=\frac{c}{{\rm H}_0(1+z)^3}, \end{equation} \tag{ 53 }$$

dl = cdt is a length element, and p is the property in consideration. The halo mass is denoted by M, N(M) is the comoving number density of halos (i.e., the halo mass function), and Σ(M) is the halo cross-section given by Equation (50). The quantity dM/dp in Equation (52) can be easily calculated in our model, as properties like metallicity and relative abundances are known for all halo masses. The integral of Equation (52) over all values of p is just the total line density of DLAs, dN/dX. The evolution of dN/dX in our model is shown in Figure 10. Our results are in reasonable agreement with the values obtained from the Sloan Digital Sky Survey (SDSS) DR5 by Prochaska et al. (2005) and from SDSS DR9 by Noterdaeme et al. (2012). (Note that we have assumed a slightly smaller value of α compared to Pontzen et al. (2008) to get a good match to observed dN/dX.)

Zoom In Zoom Out Reset image size

Using Equation (52), we plot the line density distribution of DLAs at z = 6 as a function of various metal abundance ratios in Figure 5 (right column). We focus on [O/Si] for this discussion, which can be easily generalized for other ratios in Figure 5. The panel corresponding to [O/Si] shows the DLA line density distribution for the three Population III IMFs considered in this paper (Table 1). The distribution corresponding to model 2 (35–100 M_☉; blue curve) is peaked at [O/Si] = 0.2 and has a tail spreading out to higher values up to [O/Si] = 0.8. The distribution corresponding to model 3 (100–260 M_☉; red curve) is also peaked at [O/Si] = 0.2, but has a tail spreading out to lower values down to [O/Si] = −0.4. Finally, model 1 (1–100 M_☉) is simply a delta function at [O/Si] = 0.2. The line density distribution in this model has no tail; all DLAs have the same abundance ratio of [O/Si] = 0.2. This suggests that the distribution of relative abundance values in a sample of DLAs depends on the Population III IMF, at least at sufficiently high redshifts.

This dependence of the line density distribution on Population III IMF can be understood in terms of the dependence of halo relative abundances on the Population III IMF. We discussed the latter in Sections 3.2 and 3.3 above, and in the left-hand column of Figure 5. Focusing again on [O/Si], recall that most high-mass halos had [O/Si] = 0.2, while low-mass halos had different values of [O/Si], depending on the Population III IMF. This same effect is reflected in the line density distribution in the right-hand column of Figure 5, since all of these halos contribute to the line density distribution via Equation (52). Thus, e.g., as low-mass halos in model 2 (35–100 M_☉; blue curve) have [O/Si]  >  0.2, DLAs corresponding to these halos form the tail that spreads toward these values. Similarly, in model 3 (100–260 M_☉; red curve), low-mass halos have negative [O/Si], which results in a tail in the DLA line density distribution that extends toward these values. Finally, in model 1 (1–100 M_☉), all halos have the same value of [O/Si]. As a result, the DLA line density is a delta function centered on this value (0.2) and has no tail. This also suggests that with decreasing redshift, as more and more halos move to their Population II-producing phase, the DLA line density distribution will reduce its spread and move toward the Population II value. This is exactly what is seen as we move from z = 6 to z = 2, as shown in Figure 6 (right column).

In sum, Figure 5 displays the general result that the distribution of DLAs in the abundance-ratio space at sufficiently high redshift is sensitive to the Population III IMF. In the absence of Population III stars causing any change in stellar yields over time, this distribution will be a delta function in the abundance-ratio space (neglecting corrections for the impact of dust on relative abundances, which we discuss below). However, the addition of a new Population III stellar population with different chemical yields spreads DLAs out in the abundance-ratio space. The shape of this spread can constrain the IMF of Population III stars, while its evolution can constrain the Population III SFR history. (Note that as discussed before in Section 3.3 above, the gas-phase metallicity of these DLAs is higher than Z_crit, similar to a vast majority of observed DLAs. Still their relative abundances have not yet settled on the Population II values. This results in a dependence on the Population III IMF.)

Further, to check whether this effect is detectable and to show that it can indeed be used to probe the Population III IMF, in Figure 11 we show a simulated data set of 10 [O/Si] measurements at z = 6 taken from the DLA line density distribution corresponding to model 2 (35–100 M_☉; shown in the top panel of Figure 11 and also in Figure 5). This data set is produced by random sampling the predicted distribution of [O/Si] values in the model, and adding a Gaussian error of 0.1 dex to each sample point. An error of 0.1 dex corresponds to the typical accuracy with which abundance ratios in DLAs are measured (e.g., Becker et al. 2012). Using a KS test, we find that this sample rejects the Population II distribution at 4σ (D = 0.9 and p = 0.00017). A set of 100 samples (each of size 10 and a Gaussian error of 0.1 dex) preferred the Population III IMF at 3.8σ on average (〈D〉 = 0.8, 〈p〉 = 0.007). This suggests that the effects of the Population III IMF on DLA abundance ratio are significant enough to be detectable with just 10 measurements accurate to about 0.1 dex.⁴ A larger sample of 20 measurements rejects the Population II distribution at an even higher significance of about 5σ (〈D〉 = 0.85 and 〈p〉 = 10⁻⁵; again assuming an accuracy of 0.1 dex). Note that this test used a single abundance ratio ([O/Si]). In practice, using multiple ratios as shown in Figure 5 can further improve the significance of the constraints. Also note that the metallicity of these low-mass halos (M ∼ 10⁹ M_☉) with Population III signatures is high enough to produce detectable lines at z ∼ 6. For example, we can estimate the O i column density by Salvadori & Ferrara (2012)

$$\begin{equation} N_\mathrm{OI}=\frac{3}{2\pi }\frac{M_\mathrm{O}/\mu _\mathrm{O} m_p}{\alpha ^2r_\mathrm{vir}^2}, \end{equation} \tag{ 54 }$$

where μ_O is the atomic weight of oxygen and r = αr_vir is the gas radius, where usually α ∼ 0.8 (Salvadori & Ferrara 2012). The oxygen mass M_O is calculated in each halo in our chemical evolution model. This gives N_OI > 10¹⁵ cm⁻² for α = 0.8. This value of the column density is already higher than that measured in several z = 6 DLAs by Becker et al. (2012), which span the range N_OI = 10^13.49–10^14.47 cm⁻².

Zoom In Zoom Out Reset image size

It is important to understand that the distribution of relative abundance values at high redshifts is used here, not the mean value at those redshifts. In fact, as we have seen, all line density distributions in Figure 5 are anchored on the Population II values of the respective abundance ratios. This is because even at high redshift, high-mass halos have already moved to their Population II star-forming stage; only low-mass halos are expected to carry the memory of their Population III enrichment. Therefore, we expect the mean value of any abundance ratio to be close to the Population II value. It will have very little dependence on the Population III IMF. This can be understood by explicitly calculating the mean value of different abundance ratios, given by

$$\begin{equation} \langle p\rangle =\int \frac{d^2N}{dXdp}\cdot p\cdot {dp}\cdot {dX}, \end{equation} \tag{ 55 }$$

where p = [M₁/M₂] is the abundance ratio of two species M₁ and M₂, and dX is defined in Equation (53). Figure 12 shows the mean value of [C/Fe] and [O/Si] calculated in this fashion, corresponding to our three Population III models. Observational measurements of several DLAs from Becker et al. (2012) are also shown. It is seen that the average values predicted by our model are consistent with the observations, regardless of the Population III IMF used. (Note that the observational data shown here are uncorrected for effects of dust depletion.) The most striking feature in Figure 12 is the lack of evolution over the large redshift range from z ≳ 6 to z ∼ 2. This is explained in our model by the fact that all large halos show constant abundance ratios, corresponding to Population II. A second feature of Figure 12 is that varying the Population III IMF results in only a small change in the average value of the abundance ratio, as expected. It is only when we look at the spread of these values at different redshift that the effect of Population III IMFs becomes visible.

Zoom In Zoom Out Reset image size