THE ROLE OF DUST IN THE EARLY UNIVERSE. I. PROTOGALAXY EVOLUTION

SNe II are believed to be the dominant sources of dust at a high redshift of z > 5 because of the short lifetimes (<10⁷ yr) of their massive progenitors (e.g., Dwek et al. 2007; Gallerani et al. 2010; Gall et al. 2011). Dust formation in the ejecta of primordial SNe II has been investigated theoretically (Todini & Ferrara 2001; Nozawa et al. 2003; Cherchneff & Dwek 2010). The amount and the size distributions of dust grains injected into the ISM have been investigated by considering the destruction in SN remnants (SNRs; Bianchi & Schneider 2007; Nozawa et al. 2007; Nath et al. 2008; Silvia et al. 2010). Although the amount of dust that forms in the ejecta is still debated (Kozasa et al. 2009, for review), the recent observations of Cas-A SNR revealed the presence of ∼0.07 M_☉ dust condensed in the ejecta (Barlow et al. 2010; Sibthorpe et al. 2010), which is consistent with the dust mass predicted by the theoretical model taking into account formation and destruction processes of dust in an SN IIb (Nozawa et al. 2010). Valiante et al. (2009) and Dwek & Cherchneff (2011) have proposed that the contribution from asymptotic giant branch (AGB) stars in the high-redshift quasar J1148+5251 cannot be neglected for the total dust budget even at z ∼ 6. However, the size distribution of dust formed in the mass-loss wind of AGB stars has not yet been fully studied (Ferrarotti & Gail 2006; Zhukovska et al. 2008). If SNe Ia could occur in such an early epoch, they are unlikely to be efficient sources of dust (Nozawa et al. 2011). Therefore, in order to follow the evolution of dust size distribution and reveal the resulting influence on galaxy evolution, we consider SNe II as the source of dust in the early universe.

The basic quantity for governing the production and destruction history of dust by SNe II is the rate of SN II explosions, γ_SN(t), given by

$$\begin{equation} \gamma _{\rm SN}(t)=\int _{m_{\rm SN}^{l}}^{m_{\rm SN}^{u}}{d}m\Psi (t-\tau (m))\phi (m), \end{equation} \tag{ 1 }$$

where Ψ(t) is the SFR at time t, ϕ(m) is the stellar IMF, τ(m) is the lifetime of a star whose mass is m, and m^u_SN and m^l_SN are the upper and lower mass limits of SN II progenitors, respectively. In this paper we adopt the Salpeter IMF (ϕ(m)∝m^−2.35; Salpeter 1955) with the stellar mass range between 0.1 M_☉ and 60 M_☉; we assume m^l_SN = 8 M_☉ and m^u_SN = 40 M_☉ (Heger et al. 2003). For τ(m), we adopt the model of zero-metallicity stars without mass loss (Schaerer 2002).

In this paper, we do not consider Population III stars, for simplicity. In the forthcoming paper, we will consider the possible contribution of Population III stars. Population III stars formed out of the primordial gas are considered to be much more massive than Population I/II stars (Yoshida et al. 2008; Bromm et al. 2009 for reviews); thus, the primordial IMF might be biased toward a higher mass (≳ 10 M_☉) than that in the present universe. Furthermore, Population III stars as massive as 140–260 M_☉ are predicted to end their lives as pair-instability SNe (PISNe; Heger & Woosley 2002) and to produce a large amount of metals and dust (Nozawa et al. 2003; Schneider et al. 2004). However, Joggerst et al. (2010a, 2010b) address the growing nucleosynthetic "forensic" evidence that the majority of primordial stars may have been 15–40 M_☉ objects. On the other hand, once the gas is enriched up to a critical metallicity of Z ≃ 10⁻⁶–10⁻⁵ Z_☉, formation of low-mass stars is triggered, leading to the transition of the star formation mode from massive Population III stars to low-mass Population I/II stars, if dust is present (Omukai et al. 2005; Schneider et al. 2006; Schneider & Omukai 2010; Omukai et al. 2010). If there is no dust, the transition of the star formation mode is expected to occur at 10^−3.5 Z_☉ (Mackey et al. 2003). In this case, the cosmic microwave background (CMB) limits the lower masses of stars to a few tens of M_☉ (Smith et al. 2008; Schneider & Omukai 2010). Formation history of galaxies considering the time-dependent IMF from the top-heavy to the Salpeter-like IMF and taking into account the production and destruction of dust by PISNe will be explored in the forthcoming paper (D. Yamasawa et al. 2011, in preparation).

Throughout this paper, we adopt the models by Nozawa et al. (2003, 2007) for dust formation and destruction. Nozawa et al. (2003) investigated the dust production in the ejecta of primordial SNe II as well as PISNe, applying a theory of non-steady-state nucleation and grain growth. They revealed the grain species formed in the ejecta and their size distributions for the unmixed and mixed elemental compositions within the He core. In what follows, we apply the results of calculation for the unmixed ejecta of SNe II with the progenitor mass m = 13, 20, 25, and 30 M_☉ and the explosion energy 10⁵¹ erg, and extrapolate the data to the mass range from 8 to 40 M_☉.

The solid line in Figure 1 shows the IMF-averaged mass distribution of dust formed in the ejecta of SNe II, $\overline{{\cal M}_{{d}}^{0}}(a)$, which is weighted by the Salpeter IMF and is summed up over all the grain species as

$$\begin{eqnarray} \overline{{\cal M}_{d}^{0}}(a)&=&\sum _{j}\overline{{\cal M}_{{d},j}^{0}}(a) =\frac{\sum _{j}\int _{{\rm m}_{\rm SN}^{l}}^{m_{\rm SN}^{u}}{d}m\ {\cal M}_{{d},j}^{0}(a,m)\phi (m)}{\sum _{j}\int _{m_{\rm SN}^{l}}^{m_{\rm SN}^{u}}{d}m\ \phi (m)},\nonumber\\ \end{eqnarray} \tag{ 2 }$$

where ${\cal M}_{{d},j}^{0}(a,m){d}a$ is the mass of the jth dust species produced in an SN II with radii between a and a + da as a function of progenitor mass m, and superscript 0 indicates the case with no destruction by a reverse shock. In Figure 1, we plot $a\overline{{\cal M}_{{d}}^{0}}(a)$ in the vertical axis to make clear the mass fraction in each logarithmic bin. We can see that the grain radii range from a few Å up to a few μm and that the size spectrum of dust in mass has a peak at a ∼ 0.1 μm.

Zoom In Zoom Out Reset image size

In the course of their injection into ISM, dust grains formed in the ejecta are destroyed due to sputtering in the hot gas between the reverse and forward shocks, which is hereafter referred to as destruction by the reverse shock. Nozawa et al. (2007) investigated the survival of the newly formed dust in the shocked gas within the SNRs expanding into the uniform ISM with hydrogen number densities of n_SN = 0.1, 1.0, and 10.0 cm⁻³ and showed that the destruction efficiency of newly formed dust is not only sensitive to the initial size distribution but also strongly depends on n_SN. To investigate the dependence of destruction of dust on the ISM densities, we extend their models to six cases of n_SN = 0.03, 0.1, 0.3, 1.0, 3.0, and 10.0 cm⁻³.

Figure 1 shows the IMF-averaged mass distribution of dust, $\overline{{\cal M}_{{d}}^{n_{\rm SN}}}(a)$, injected into the ISM after destruction by the reverse shock for n_SN = 0.03, 0.1, 0.3, 1.0, 3.0, and 10.0 cm⁻³, which is weighted by the Salpeter IMF and is summed up over all the grain species as in Equation (2). The mass distribution of the jth dust species for the case of n_SN, $\overline{{\cal M}_{{d},j}^{n_{\rm SN}}}(a)$, is the IMF-averaged one after the destruction by the reverse shock. We can see that the change in the dust mass distribution through processing in SNRs becomes more (less) prominent for higher (lower) gas density; small-size grains get deficient with increasing n_SN. The dust grains with radii below 0.01 μm are preferentially destroyed by sputtering for n_SN > 0.03 cm⁻³, while dust with radii larger than ∼1 μm is almost intact for n_SN ⩽ 3.0 cm⁻³. As a result, the mass of dust injected into the ISM is dominated by grains with radii above ∼0.1 μm.

The total geometrical cross section of dust per metal mass is an important quantity for H₂ formation on the grain surface at a certain metallicity level. The total geometrical cross section of dust injected into the ISM from an SN depends on the gas density around the SN progenitor, n_SN. In many papers, the total geometrical cross section of dust is scaled to metal mass under the following two assumptions: (1) the depletion factor, which is defined as the dust mass per metal mass, is identical to that in the Milky Way (MW); and (2) the dust size distribution, which determines dust area per unit dust volume 〈a²〉/〈a³〉, is the same as that in the MW. However, the depletion factor and the size distribution of the SN II dust are quite different from those in the MW.

Figure 2 shows the depletion factor, $\overline{M_{\rm SN,d}^{n_{\rm SN}}}/\overline{m_{m}}$, after the destruction by the reverse shock as a function of ISM density nSN, where MnSNSN, d is the IMF-averaged total dust mass ejected into the ISM by an SN,

$$\begin{equation} \overline{M_{\rm SN,d}^{n_{\rm SN}}}=\int _{0}^{\infty }{d}a\overline{{\cal M}_{d}^{n_{\rm SN}}}(a), \end{equation} \tag{ 3 }$$

and $\overline{m_{m}}$ is the IMF-averaged total metal mass ejected into the ISM,

$$\begin{eqnarray} \overline{m_{m}}&=&\sum _{i}\overline{m_{{\rm m},i}} =\frac{\sum _{i}\int _{m_{\rm SN}^{l}}^{m_{\rm SN}^{u}}m_{{\rm m},i}(m)\phi (m){d}m}{\sum _{i}\int _{m_{\rm SN}^{l}}^{m_{\rm SN}^{u}}\phi (m){d}m}. \end{eqnarray} \tag{ 4 }$$

Here, m_{m, i}(m) is the mass of the ith element of metal ejected from an SN with progenitor m, taken from Umeda & Nomoto (2002). The depletion factor in the case of n_SN = 1.0 cm⁻³ is a factor of six smaller than that in the MW.

Zoom In Zoom Out Reset image size

Figure 3 shows the ratio of the total geometric cross section of dust to the total metal mass, $\overline{D_{\rm area}^{n_{\rm SN}}}/\overline{m_{\rm m}}$, after the destruction by the reverse shock as a function of ISM density n_SN. The total geometrical cross section of dust weighted by the Salpeter IMF, $\overline{D_{\rm area}^{n_{\rm SN}}}$, is written as

$$\begin{equation} \overline{D_{\rm area}^{n_{\rm SN}}}=\sum _{j}\frac{3}{4\rho _{j}}\int _{0}^{\infty }{d}a\overline{{\cal M}_{{d},j}^{n_{\rm SN}}}(a)/a, \end{equation} \tag{ 5 }$$

where ρ_j is the bulk density of the jth dust species. The ratio, $\overline{D_{\rm area}^{n_{\rm SN}}}/\overline{m_{\rm m}}$, is smaller for larger n_SN because small dust grains are efficiently destroyed by sputtering in the SNR under large n_SN; note that the surface area per dust mass is larger for smaller size grains. Also with increasing n_SN, the reverse shock becomes stronger and destroys dust by sputtering more effectively (see Nozawa et al. 2007). In addition, we plot the ratio for the case without reverse shock together with the typical value in the MW for comparison. $\overline{D_{\rm area}^{n_{\rm SN}}}/\overline{m_{\rm m}}$ in the model without reverse shock is a factor of four smaller than that in the MW. $\overline{D_{\rm area}^{n_{\rm SN}}}/\overline{m_{\rm m}}$ in the case with n_SN = 1.0 cm⁻³ is 40 times smaller than that in the MW. Thus, the rescaling of the cross section of dust by the metal mass using the MW value results in significant overestimate for H₂ formation in high-redshift galaxies.

Zoom In Zoom Out Reset image size

Dust grains injected into the ISM are subjected to destruction by the blast waves (the high-velocity interstellar shocks) driven by the ambient SNe (e.g., Jones et al. 1994). Nozawa et al. (2006) investigated the processing of interstellar dust by sputtering in the hot gas swept up by the SN forward shock. Adopting the dust model by Nozawa et al. (2003) as the size distribution of interstellar dust, Nozawa et al. (2006) have shown that the destruction efficiency of dust depends on the ISM density and the explosion energy of SNe as well as the initial size distribution of dust. It should be noted that the size distribution and the destruction efficiency change as functions of time because interstellar dust is supplied from SNe and processed in ISM successively according to star formation activity. Thus, we must deal with the destruction process in a way that is applicable to any dust size distribution to explore the global evolution of dust size distribution.

In order to evaluate the destruction efficiency of interstellar dust for any initial size distribution, here we introduce the conversion efficiency as defined below. Consider that the jth dust species residing in the ISM, whose size distribution is given by the number of dust grains with radii between a and a + da, f_j(a)da, is processed by sputtering in hot plasma produced through a single passage of an SN shock. The conversion efficiency η_j(a, a') is defined as the number fraction of dust grains with radii between a' and a' + da' that are converted to grains with radii between a and a + da by sputtering through the passage of an SN shock. The number of dust grains with radii between a and a + da produced by the sputtering is given as η_j(a, a')f_j(a')da'. Note that η_j(a, a') = 0 for a > a'. Then the change in the number of dust grains with radii between a and a + da caused by shock processing is given by

$$\begin{eqnarray} {d}N_{j}(a)&=&\sum _{a^{\prime }>a}^{\infty }\eta _{j}(a,a^{\prime })f_{j}(a^{\prime }){d}a^{\prime } -[1-\eta _{j}(a,a)]f_{j}(a){d}a \nonumber \\ &=&\int _{0}^{\infty }\eta _{j}(a,a^{\prime })f_{j}(a^{\prime }){d}a^{\prime }-f_{j}(a){d}a. \end{eqnarray} \tag{ 6 }$$

The corresponding change of the mass is given by

$$\begin{eqnarray} {d}M_{{d},j}(a)&=&\frac{4\pi }{3}a^{3}\rho _{j}\int _{0}^{\infty }\eta _{j}(a,a^{\prime })f_{j}(a^{\prime }){d}a^{\prime } -{\cal M}_{{d},j}(a){d}a,\nonumber\\ \end{eqnarray} \tag{ 7 }$$

where ${\cal M}_{{d},j}(a){d}a$ is the mass of the pre-shocked dust. The size distribution function after the shock processing f'_j(a) is given by f'_j(a) = f_j(a) + dN_j/da.

The conversion efficiency η(a, a') and the mass of ISM gas swept up by shock M_swept depend on the progenitor mass, expanding energy, and type of the SN as well as the structure, number density, and metallicity of the ISM gas. For these parameters of SNe and ambient ISM, once M_swept and η(a, a') for each dust species are calculated, the time evolution of dust mass and size distribution can be followed in a consistent way with the star formation activity in galaxies as described in Section 3.

The calculations of η(a, a') and M_swept are performed using the method developed by Nozawa et al. (2006) as follows: the efficiency of dust destruction increases with increasing explosion energy and/or increasing n_SN but is almost independent of the progenitor mass as long as the explosion energy is the same (Nozawa et al. 2006). We assume that SNe driving high-velocity shock in ISM are represented by SN II with the progenitor mass of 20 M_☉ and the explosion energy of 10⁵¹ erg. The ISM surrounding the SN is considered to be uniform with hydrogen number densities n_SN = 0.03, 0.1, 0.3, 1.0, 3.0, and 10 cm⁻³. By distributing dust grains with radius a' uniformly in the ISM, the conversion efficiency η(a, a') is evaluated for each grain species by calculating the erosion of dust by sputtering until the truncation time t_tr, which is defined as a time when the shock velocity is decelerated below 100 km s⁻¹ (see Nozawa et al. 2006, for the details). In the calculations, the radii of grains in the ISM range from 0.00013 to 6.3 μm for each grain species.

In Figure 4, we present the changes in the dust size distributions due to the shock processing for different ISM densities, where the initial mass distribution of dust is set to be $a{\cal M}_{{d},j}(a)=1$ for clarity. As seen in the figure, small-size grains are destroyed due to the erosion by sputtering, and more dust grains are processed for a higher ISM density.

Zoom In Zoom Out Reset image size

The mass of gas swept up by the forward shock until the truncation time t_tr, M_swept, depends not only on the ISM density but also on the initial metallicity of the gas in the ISM. As the line cooling by heavy elements becomes more efficient for a higher gas metallicity, the forward shock is decelerated more quickly, resulting in a smaller M_swept. By fitting M_swept calculated for different n_SN and Z, we derive the following approximation formula,

$$\begin{equation} M_{\rm swept}/M_{\odot }=1535n_{\rm SN}^{-0.202}\left[\left(Z/Z_{\odot }\right)+0.039\right]^{-0.298}, \end{equation} \tag{ 8 }$$

whose fitting accuracy is within 16% for 0.03 cm⁻³ ⩽ n_SN ⩽ 30 cm⁻³ and for 10⁻⁴ ⩽ Z/Z_☉ ⩽ 1.0.

In terms of the conversion efficiency describing the processing of dust by sputtering, here we formulate the time evolution of the mass of jth dust grains with radii between a and a + Δa in our model galaxies, $\Delta M_{{d},j}(a,t)=\frac{4\pi }{3}a^{3}\rho _{j}f_{j}(a,t)\Delta a$, as

$$\begin{eqnarray} \frac{{d}\Delta M_{{d},j}(a,t)}{{d}t}&=&\overline{\Delta M_{{{\rm SN},d},j}}(a)\gamma _{\rm SN}(t)\nonumber\\ &&-\frac{M_{\rm swept}}{M_{\rm ISM}(t)}\gamma _{\rm SN}(t) \times \bigg\lbrace \Delta M_{{d},j}(a,t)\nonumber\\ && -\int _{0}^{\infty }{d}a^{\prime }\eta _{j}(a,a^{\prime })f_{j}(a^{\prime },t)\rho _{j}\frac{4\pi }{3}a^{3}\bigg\rbrace \nonumber \\ && -\Psi (t)\frac{\Delta M_{{d},j}(a,t)}{M_{\rm ISM}(t)}, \end{eqnarray} \tag{ 9 }$$

where M_ISM(t) is the total mass of gas and dust, and f_j(a, t) is the size distribution function of dust species j in the ISM at a time t. We note that (M_sweptγ_SN(t)/M_ISM(t))⁻¹ is a timescale of sweeping the whole ISM by SNe. The IMF-averaged mass of dust species j with radii between a and a + Δa injected from SNe II into the ISM with number density n_{ISM, SN} is defined by $\overline{\Delta M_{{\rm SN},d,j}}(a)=\overline{{\cal M}_{d,j}^{n_{\rm SN}}}(a)\Delta a$. The first term on the right-hand side is the injection rate of dust from SNe II. The second term is the destruction rate of interstellar dust by SN blast waves, and the third term is the rate at which the interstellar dust is incorporated into stars.

We quantify the properties of dark matter halos, assuming a dynamically equilibrium state. The radius of dark matter halo, r_vir, is estimated in terms of the mass of dark halo, M_vir, and the redshift of virialization, z_vir, as

$$\begin{equation} \frac{4}{3}\pi r_{\rm vir}^{3}\left\lbrace 1+\delta _{\rm c}(z_{\rm vir})\right\rbrace \rho _{{\rm c}0}\Omega _{M}(1+z_{\rm vir})^{3}=M_{\rm vir}, \end{equation} \tag{ 10 }$$

where ρ_c0 ≡ 3H²₀/8πG is the critical density of the universe at z = 0, δ_c(z_vir) is the overdensity of a dark matter halo virialized at z_vir, and G is the gravitational constant.

We assume dark matter halos as singular isothermal spheres and rotating uniform gas disks in their gravitational potentials. Cosmological N-body simulations show that structures of dark matter halos are described well by the Navarro–Frenk–White (NFW) profile (Navarro et al. 1996). Mo et al. (1998) studied a simple disk model in the gravitational potential of the singular isothermal sphere and a more realistic disk model in the gravitational potential of the NFW halo profile. We adopt a radius of the disk, r_disk ≃ 0.18r_vir (Ferrara et al. 2000; Hirashita & Ferrara 2002), by considering the conservation of angular momentum and assuming a typical value for the spin parameter λ = 0.04 from the paper by Ferrara et al. (2000), who estimate the radius of the disk as r_disk = 4.5λr_vir in a modified isothermal halo.

In our one-zone model, we need a virial temperature for the initial gas temperature and a dynamical timescale of gas in the disk. Gas collapsed at z_vir in the dark matter halo of M_vir has a virial temperature, T_vir, defined as

$$\begin{equation} T_{\rm vir}\equiv \frac{G\mu \rm m_{\rm H}M_{\rm vir}}{3k_{\rm B}r_{\rm vir}}, \end{equation} \tag{ 11 }$$

where k_B is the Boltzmann constant, μ is the mean molecular weight, and m_H is the mass of a hydrogen atom. The initial value for the temperature of gas, T, is assumed to be T_vir. A circular velocity, v_c, is defined as

$$\begin{equation} v_{\rm c}\equiv \left(\frac{{\it GM}_{\rm vir}}{r_{\rm vir}}\right)^\frac{1}{2}, \end{equation} \tag{ 12 }$$

and we also define a rotation timescale, t_cir, as

$$\begin{equation} t_{\rm cir}\equiv \frac{2\pi r_{\rm disk}}{v_{\rm c}}. \end{equation} \tag{ 13 }$$

Note that the rotation timescale of the gas disk, t_cir, depends only on virialization redshift, z_vir, as

$$\begin{equation} t_{\rm cir}=9.3\times 10^{7}\ {\rm yr} \left(\frac{11}{1+z_{\rm vir}}\right)^{\frac{3}{2}}\left(\frac{18\pi ^{2}}{\delta _{c}(z_{\rm vir})}\right)^{\frac{1}{2}}. \end{equation} \tag{ 14 }$$

We must estimate the number density of the hydrogen gas, n_H, because it affects both the chemical reaction rate and the cooling rate. The cooling time of halo gas is much shorter than the Hubble timescale for the objects of interest in this paper (T_vir ≳ 10000 K; e.g., Hutchings et al. 2002). It is widely understood that most z ∼ 10 galaxies are not clear disk galaxies; numerical simulations that proceed from cosmological initial conditions (z ∼ 100–200) clearly reveal that they possess highly irregular structures whose SFRs are not easily quantifiable, that filamentary accretion and frequent mergers are still churning the halo at this epoch, and that turbulent flows that arise in the center of the halo prevent coherent disks from forming on the spatial scales of galaxies (Johnson et al. 2008; Greif et al. 2010; Wise et al. 2010). We assume for simplicity that a significant fraction of baryons finally collapses to a disk in the dark matter halo potential. In semi-analytic models, it is assumed that the cooled halo gas settles into the disk (e.g., Cole et al. 2000). We make a similar assumption in our model, but a more detailed treatment for H₂ formation, dust evolution, and star formation in the gas disk. The radius of disk, r_disk, is determined following Hirashita & Ferrara (2002). We estimate the typical scale height, H, from hydrostatic equilibrium (Shakura & Sunyaev 1988)

$$\begin{eqnarray} H &=& \sqrt{2}\frac{v_{\rm s}}{v_{\rm c}}r_{\rm disk} = \left(\frac{2T}{3T_{\rm vir}}\right)^{\frac{1}{2}}r_{\rm disk} \end{eqnarray} \tag{ 15 }$$

for H/r_disk ⩽ 0.1; otherwise, we assume H/r_disk = 0.1, where $v_{s}=(k_{B}T/\mu \rm m_{\rm H})^{\frac{1}{2}}$ is the isothermal sound velocity. Therefore, the initial hydrogen density of disk, n_H, is estimated as

$$\begin{equation} n_{\rm H}=\frac{M_{\rm H}}{\pi r_{\rm disk}^{2}2Hm_{\rm H}}. \end{equation} \tag{ 16 }$$

The initial mass of hydrogen in the galaxy, M_H, is written as

$$\begin{eqnarray} M_{\rm H}&=&M_{\rm gas}-M_{\rm He} =M_{\rm gas}\frac{m_{\rm H}}{(m_{\rm H}+m_{\rm He}y_{\rm He})} \nonumber \\ &=&\frac{M_{\rm vir}\Omega _{b}}{\Omega _{M}}\frac{m_{\rm H}}{(m_{\rm H}+m_{\rm He}y_{\rm He})}, \end{eqnarray} \tag{ 17 }$$

where M_He is the initial mass of helium in the galaxy, M_gas = M_H + M_He, m_He is the mass of a helium atom, and y_He is the helium abundance. We assume y_He = 0.0972 (Omukai 2000). Note that since massive stars ionize surrounding gas and form an expanding H ii region (e.g., Whalen et al. 2004; Kitayama et al. 2004), we adopt gas density around the SN progenitor, n_SN, as being different from the hydrogen gas density in our one-zone galaxy model, n_H. This is because the gas density around the SN progenitor, n_SN, is closely related to the dust destruction process by reverse shocks driven by SNe, as shown in Section 2.2.

We expect that the SFR, Ψ(t), is roughly proportional to t⁻¹_cir, since a representative timescale of the dynamics of the gas disk is t_cir. We assume that

$$\begin{equation} \Psi (t)=\frac{f_{{\rm H}_{2}}(t)M_{\rm H}(t)}{t_{\rm cir}(z_{\rm vir})}, \end{equation} \tag{ 18 }$$

where $f_{\rm H_{2}}$ is the mass fraction of molecular hydrogen to the total gas. We should note that observationally Bigiel et al. (2008) find that H₂ is converted into stars at a constant efficiency in nearby spirals, and Gnedin et al. (2009) show that the star formation recipe in a galaxy formation simulation in which star formation occurs only in the molecular gas can reproduce the observational correlations between SFR and the total gas density.

We calculate the time evolutions of the masses of hydrogen and helium gases, M_gas, stars, M_star, and metal of element i, M_{m, i}, in the galaxy using the following equations:

$$\begin{eqnarray} \frac{{d}M_{\rm gas}(t)}{{d}t}&=&{-}\Psi (t)\frac{M_{\rm gas}(t)}{M_{\rm ISM}(t)}+\overline{m_{\rm gas}}\gamma _{\rm SN}(t) \nonumber \\ \frac{{d}M_{\rm star}(t)}{{d}t}&=&\Psi (t)-\overline{m_{\rm ejecta}}\gamma _{\rm SN}(t) \nonumber \\ \frac{{d}M_{{\rm m},i}}{{d}t}&=&{-}\Psi (t)\frac{M_{{\rm m},i}(t)}{M_{\rm ISM}(t)}+\overline{m_{{\rm m},i}}\gamma _{\rm SN}(t), \end{eqnarray} \tag{ 19 }$$

where M_ISM(t) = M_gas(t) + ∑_iM_{m, i}(t), $\overline{m_{\rm ejecta}} = \overline{m_{\rm gas}}+ \sum _{i}\overline{m_{{\rm m},i}}$ and $\overline{m_{\rm gas}}$ and $\overline{m_{{\rm m},i}}$ are the gas mass of hydrogen and helium and the metal mass of element i in SN ejecta, respectively. The mass returning to the ISM, m_ejecta(m), m_gas(m), and m_{m, i}(m) through an SN with progenitor mass, m, is taken from Umeda & Nomoto (2002) in the case of m = 13, 20, 25 and 30 M_☉. $\overline{m_{\rm gas}}$, $\overline{m_{\rm ejecta}}$, and $\overline{m_{m,i}}$ are IMF-averaged m_gas(m), m_ejecta(m), and m_{m, i}(m), respectively. Note that metals consist not only of heavy elements in gas phase but also those in dust grains.

We follow the time evolution of molecular mass fraction, $f_{{\rm H}_{2}}$, ionization degree, x, and gas temperature, T. We define the molecular fraction of hydrogen as

$$\begin{equation} f_{{\rm H}_{2}}\equiv \frac{2n_{{\rm H}_{2}}}{n_{\rm H}}, \end{equation} \tag{ 20 }$$

where $n_{{\rm H}_{2}}$ and n_H are the number densities of molecular hydrogen and hydrogen nuclei, respectively. The molecular fraction is very important in our models, because it determines the final cooling rate of low-metallicity gas. The metal-free gas evolution with chemical reactions and cooling is studied using the model by Tegmark et al. (1997), Hutchings et al. (2002), and Hirashita & Ferrara (2002). In Table 1, we summarize chemical reactions considered in this paper and their rate coefficients (R_n; n = 1, ..., 11). The equations are based on Hirashita & Ferrara (2002), but we include the effect of the dust size distribution on H₂ formation and the metal-line cooling process.

Table 1. Reaction Rates Needed to Calculate the Abundance of H₂

No.	Reaction	Rate	Ref.
		(cm3 s−1)
1	${{\rm H}+e^{-}\longrightarrow {\rm H}^++2e^-}$	exp [ − 32.71 + 13.54ln (T(eV))	1
		−5.739(ln (T(eV)))2 + 1.563(ln (T(eV)))3
		−0.2877(ln (T(eV)))4 + 3.483 × 10−2(ln (T(eV)))5
		−2.632 × 10−3(ln (T(eV)))6
		+1.120 × 10−4(ln (T(eV)))7
		−2.039 × 10−6(ln (T(eV)))8]
2	${{\rm H}^++e^-\longrightarrow {\rm H}+\gamma }$	exp [ − 28.61 − 0.7241(ln (T(eV)))	1
		−2.026 × 10−2(ln (T(eV)))2
		−2.381 × 10−3(ln (T(eV)))3
		−3.213 × 10−4(ln (T(eV)))4
		−1.422 × 10−5(ln (T(eV)))5
		+4.989 × 10−6(ln (T(eV)))6
		+5.756 × 10−7(ln (T(eV)))7
		−1.857 × 10−8(ln (T(eV)))8
		−3.071 × 10−9(ln (T(eV)))9]
3	${{\rm H}+e^{-}\longrightarrow {\rm H}^-+\gamma }$	1.4 × 10−18T0.928exp (− T/1.62 × 104)	1
4	${{\rm H}^-+{\rm H}\longrightarrow {\rm H}_2+e^-}$	4.0 × 10−9T−0.17(T > 300);	1
		1.5 × 10−9(T < 300)
5	${\rm H^-+H^+\longrightarrow 2H}$	5.7 × 10−6T−1/2 + 6.3 × 10−8	1
		−9.2 × 10−11T1/2 + 4.4 × 10−13T
6	${\rm H+H^+\longrightarrow H_2^++\gamma }$	dex[ − 19.38 − 1.523log10T	1
		+1.118(log10T)2 − 0.1269(log10T)3]
7	${\rm H_2^++H\longrightarrow H_2+H^+}$	6.4 × 10−10	1
8	${{\rm H}_2^++e^-\longrightarrow 2{\rm H}}$	2.0 × 10−7T−1/2	1
9	${\rm H_2+H^+\longrightarrow H_2^++H}$	3.0 × 10−10exp (− 21050/T) (T < 104)	2
		1.5 × 10−10exp (− 14000/T) (T > 104)
10	${\rm H_2+H\longrightarrow 3H}$	k1 − aHkLa	1
		kL = 1.12 × 10−10exp (− 7.035 × 104/T)
		kH = 6.5 × 10−7T−1/2
		× exp (− 5.2 × 104/T)[1 − exp (− 6000/T)]
		a = 4.0 − 0.416log10(T/104) − 0.327(log10(T/104))2
11	${{\rm H}_2+e^-\longrightarrow 2{\rm H}+e^-}$	4.4 × 10−10T0.35exp (− 1.02 × 105/T)	1
Dust	${\rm H+H+\mbox{grain}\longrightarrow H_2+\mbox{grain}}$	See Section 3.5

Note. The unit of the gas temperature T is K unless otherwise stated. References. (1) Omukai 2000; (2) Galli & Palla 1998.

Download table as:  ASCIITypeset image

The time evolution of the ionization degree is described as

$$\begin{equation} \frac{{d}x}{{d}t}=xf_{0}R_{1}n_{\rm H}-x^{2}R_{2}n_{\rm H}+\Gamma _{12}f_{0}, \end{equation} \tag{ 21 }$$

where $f_{0}=1-x-f_{{\rm H}_{2}}$ is the neutral fraction of hydrogen. The terms on the right-hand side are the rates of collisional ionization, recombination, and photoionization. Next, the time evolution of the molecular fraction is written as

$$\begin{eqnarray} \frac{{d}f_{{\rm H}_{2}}}{{d}t}&=&\left[\frac{{d}f_{{\rm H}_{2}}}{{d}t}\right]_{\rm gas}+\left[\frac{{d}f_{{\rm H}_{2}}}{{d}t}\right]_{\rm dust}+\left[\frac{{d}f_{{\rm H}_{2}}}{{d}t}\right]_{\rm dest} \nonumber \\ &&+\left[\frac{{d}f_{{\rm H}_{2}}}{{d}t}\right]_{\rm UV}+\left[\frac{{d}f_{{\rm H}_{2}}}{{d}t}\right]_{\rm star}, \end{eqnarray} \tag{ 22 }$$

where the terms on the right-hand side are the H₂ formation rate in the gas phase, the H₂ formation rate on dust grains, the destruction rate in the gas phase, the destruction rate by UV photons, and the decreasing rate by star formation, respectively. These terms are given by

$$\begin{eqnarray*} \left[\frac{{d}f_{{\rm H}_{2}}}{{d}t}\right]_{\rm gas}&=&2f_{0}^{2}xn_{\rm H}(R_{{\rm eff},1}+R_{{\rm eff},2}), \nonumber \\ \left[\frac{{d}f_{{\rm H}_{2}}}{{d}t}\right]_{\rm dust}&=&2R_{\rm dust}{\cal D}n_{\rm H}f_{0}, \nonumber \\ \left[\frac{{d}f_{{\rm H}_{2}}}{{d}t}\right]_{\rm dest}&=&{-}f_{{\rm H}_{2}}n_{\rm H}(x^{2}R_{{\rm eff},3}+f_{0}R_{10}+xR_{11}), \nonumber \\ \left[\frac{{d}f_{{\rm H}_{2}}}{{d}t}\right]_{\rm UV}&=&{-}\Gamma _{13}f_{{\rm H}_{2}}, \nonumber \end{eqnarray*}$$

and

$$\begin{eqnarray} \left[\frac{{d}f_{{\rm H}_{2}}}{{d}t}\right]_{\rm star}&=&{-}(1-f_{{\rm H}_{2}})\Psi (t)\frac{1}{M_{\rm ISM}(t)}. \end{eqnarray} \tag{ 23 }$$

The effective formation rates of H₂ including the effect of the destruction rate of H⁻ and H⁺₂ are

$$\begin{equation} R_{{\rm eff},1}\equiv \frac{R_{3}R_{4}}{f_{0}R_{4}+xR_{5}+\Gamma _{14}/n_{\rm H}} \end{equation} \tag{ 24 }$$

and

$$\begin{equation} R_{{\rm eff},2}\equiv \frac{R_{6}R_{7}}{f_{0}R_{7}+xR_{8}+\Gamma _{15}/n_{\rm H}}, \end{equation} \tag{ 25 }$$

respectively, and the destruction of H⁺₂ due to H⁻ collision is

$$\begin{equation} R_{{\rm eff},3}\equiv \frac{R_{8}R_{9}}{f_{0}R_{7}+xR_{8}+\Gamma _{15}/n_{\rm H}}. \end{equation} \tag{ 26 }$$

We will give the dust-to-gas mass ratio, ${\cal D}$, and the production rate of molecular hydrogen via dust surface reaction, R_dust, in Section 3.5 and reaction rates of the photo-process, Γ_n(n = 12, ..., 15), in Section 3.6.

At temperature <10⁴ K, the main coolant is molecular hydrogen in low-metallicity gas. The cooling rate for molecular hydrogen, $\Lambda _{{\rm H}_{2}}$, over the range 10 K ⩽ T ⩽ 10⁴ K is given by (Galli & Palla 1998)

$$\begin{eqnarray} \log _{10}\left(\frac{\Lambda _{{\rm H}_{2}}(T)}{n_{\rm H}n_{{\rm H}_{2}}\ {\rm erg}\ {\rm cm}^{3}\ {\rm s}^{-1}}\right)\! &=&\! {-}103.0+97.59T_{\rm log}-48.05T_{\rm log}^{2} \nonumber \\ && {+}10.80T_{\rm log}^{3}\,{-}\,0.9032T_{\rm log}^{4}, \end{eqnarray} \tag{ 27 }$$

where T_log ≡ log₁₀(T/K). Glover & Abel (2008) have recently given the H₂ cooling rates, which include H − H₂ collision and H₂ − H₂ collision pathways, while the Galli & Palla (1998) rates include only H − H₂ collisions. We assume that the lower limit of gas temperature is the CMB temperature.

At temperature T ≳ 10⁴ K, collisional excitation, Λ_{H, ce}, and (less importantly) ionization of atomic hydrogen, Λ_{H, ci}, are more dominant cooling processes than molecular hydrogen cooling and are given by (Haiman et al. 1996)

$$\begin{equation} \frac{\Lambda _{\rm H,ce}(T)}{n_{{\rm e}^{-}}n_{\rm H}\ {\rm erg}\ {\rm cm}^{3}\ {\rm s}^{-1}}=7.50\times 10^{-19}\frac{1}{1+T_{5}^{\frac{1}{2}}}\ \exp ^{-\frac{1.183}{T_{5}}} \end{equation} \tag{ 28 }$$

and

$$\begin{equation} \frac{\Lambda _{\rm H,ci}(T)}{n_{{\rm e}^{-}}n_{\rm H}\ {\rm erg}\ {\rm cm}^{3}\ {\rm s}^{-1}}=4.02\times 10^{-19}\frac{T_{5}^{\frac{1}{2}}}{1+T_{5}^{\frac{1}{2}}}\ \exp ^{-\frac{1.578}{T_{5}}}, \end{equation} \tag{ 29 }$$

respectively, where T₅ is the gas temperature in units of 10⁵ K.

We consider fine-structure cooling by C_I, C_II, and O_I, which dominates the thermal evolution for the number density of gas of interest in this paper (Omukai et al. 2005). The related parameters of transitions are given in Hollenbach & McKee (1989).

The increasing rate of molecular fraction via dust surface reaction is estimated as

$$\begin{eqnarray} \left[\frac{{d}f_{{\rm H}_{2}}}{{d}t}\right]_{\rm dust}&=&2R_{\rm dust}{\cal D}n_{\rm H}f_{0} =\sum _{j}\int _{0}^{\infty } f_{0}f_{j}(a)\pi a^{2}\bar{v}S {d}a, \nonumber\\ \end{eqnarray} \tag{ 30 }$$

where $\bar{v}$ is the mean thermal speed of hydrogen and S is the sticking efficiency of hydrogen atoms. We assume that the gas follows a Maxwellian distribution so that thermal speed is given by (Krügel 2008)

$$\begin{equation} \bar{v}=\left(\frac{8}{\pi }\frac{k_{\rm B}T}{m_{\rm H}}\right)^{\frac{1}{2}}. \end{equation} \tag{ 31 }$$

Here, we define the dust-to-gas mass ratio, ${\cal D}$, as

$$\begin{equation} {\cal D}\equiv \sum _{j}\int _{0}^{\infty }\frac{4\pi a^{3}\rho _{j}f_{j}(a)}{3n_{\rm H}m_{\rm H}}{d}a. \end{equation} \tag{ 32 }$$

The reaction rate of the H₂ formation on grains, R_dust, can be estimated as

$$\begin{equation} R_{\rm dust}(a){\cal D}=\sum _{j}\int _{0}^{\infty }\left(\frac{3m_{\rm H}\bar{v}S}{8a\rho _{j}}\right)\left(\frac{4\pi a^{3}\rho _{j}f_{j}(a)}{3n_{\rm H}m_{\rm H}}\right){d}a. \end{equation} \tag{ 33 }$$

We adopt S = 0.2 for T < 300 K and S = 0 for T > 300 K (Hirashita & Ferrara 2002).

We follow photo-processes in chemical reactions and heating processes in thermal evolution by using the models developed by Kitayama & Ikeuchi (2000) and improved in Hirashita & Ferrara (2002). The intrinsic luminosity is assumed to be equal to the total luminosity of OB stars whose mass is larger than 3 M_☉ (Cox 2000):

$$\begin{equation} L_{\rm UV,0}(t)=\int _{3\, M_{\odot }}^{\infty }{d}m\int _{0}^{\tau _{m}}{d}t^{\prime }\ L(m)\phi (m)\Psi (t-t^{\prime }), \end{equation} \tag{ 34 }$$

where L(m) is the stellar luminosity as a function of stellar mass m. For L(m), we adopt the model of zero-metallicity stars without mass loss in Schaerer (2002). We assume the spectrum of the incident UV radiation from stars is a power law with an index α:

$$\begin{equation} I_{\rm UV}(\nu)=I_{0}(\nu _{\rm H\,\mathsc{i}})\left(\frac{\nu }{\nu _{\rm H\,\mathsc{i}}}\right)^{-\alpha }, \end{equation} \tag{ 35 }$$

where ν is the frequency of photons and $I_{0}(\nu _{\rm H\,\mathsc{i}})$ is the intensity at the ionization threshold frequency of neutral hydrogen $\nu _{\rm H\,\mathsc{i}}=3.3\times 10^{15}\ {\rm Hz}$. In this paper, we simply set α = 5 according to Hirashita & Ferrara (2002). The normalization of the intensity is determined by

$$\begin{equation} \frac{L_{\rm UV,0}\exp (-\tau _{\rm disk})}{4\pi r_{\rm disk}^{2}}=\int _{\nu _{\rm min}}^{\infty }I_{\rm UV}(\nu){d}\nu, \end{equation} \tag{ 36 }$$

where ν_min is the minimum frequency where OB stars dominate the radiative energy of star-forming galaxies, and τ_disk is the typical dust optical depth in the disk. We assume that ν_min = 10¹⁵ Hz. This typical optical depth can be simply estimated as

$$\begin{equation} \tau _{\rm disk}=r_{\rm disk}\sum _{j}\int _{0}^{\infty }\pi a^{2}f_{j}(a){d}a \end{equation} \tag{ 37 }$$

by assuming that extinction efficiency of dust is unity in UV. We calculate Γ₁₂, Γ₁₄, and Γ₁₅ from Equation (A20) of Kitayama & Ikeuchi (2000) and the heating rate from Equation (A21) of Kitayama & Ikeuchi (2000). We summarize the cross section for the photo-process in Table 2. The H₂ photodissociation cross section is given by Abel et al. (1997). However, if the H₂ column density becomes larger than 10¹⁴ cm⁻², self-shielding effects become important (Draine & Bertoldi 1996). Therefore, the H₂ dissociation rate, Γ₁₃, is given by (Hirashita & Ferrara 2002)

$$\begin{eqnarray} \Gamma _{13}&=&(4\pi)1.1\times 10^{8}I_{\rm UV}(3.1\times 10^{15}\ {\rm Hz}) \nonumber \\ &&\times \left(\frac{n_{\rm H}f_{{\rm H}_{2}}r_{\rm disk}}{10^{14}\ {\rm cm}^{-2}}\right)^{-0.75}\ {\rm s}^{-1}, \end{eqnarray} \tag{ 38 }$$

where I_UV(3.1 × 10¹⁵ Hz) is in the Lyman–Werner band. We should note that adoption of r_disk in Equation (38) gives an extreme upper bound to the self-shielding, so we may overestimate self-shielding to internal Lyman–Werner photons by H₂. In this paper, we focus on the effects of dust on the protogalaxy, so for simplicity, we set an extreme upper bound to the self-shielding.

Table 2. Cross Sections for the Photoionization and Photodissociation Process, where the Frequency, ν, is in Units of Hz

No.	Reaction	Cross Section	ν Range	Ref.
		(cm2)	(Hz)
12	${{\rm H}+\gamma \longrightarrow {\rm H}^++e^-}$	6.30 × 10−18(ν/3.3 × 1015)−3.0	ν > 3.3 × 1015	1
13	${\rm H_2+\gamma \longrightarrow H_2^*\longrightarrow 2H}$	see Equation (38)		2
14	${{\rm H}^-+\gamma \longrightarrow {\rm H}+e^-}$	3.486 × 10−16(x − 1)3/2/x3.11	ν > 1.8 × 1014	3
		(x ≡ ν/1.8 × 1014)
15	${\rm H_2^++\gamma \longrightarrow H+H^+}$	7.401 × 10−18	ν > 6.4 × 1014	3
		dex(− x2 − 0.0302x3 − 0.0158x4)
		(x ≡ 2.762ln (ν/2.7 × 1015)

References. (1) Kitayama & Ikeuchi 2000; (2) Abel et al. 1997; (3) Tegmark et al. 1997.

Download table as:  ASCIITypeset image

We do not consider the Lyman–Werner background, since we concentrate on the evolution of atomic line cooling halos (M_vir > 10⁸ M_☉ in z < 10) in which destruction of molecular hydrogen by the Lyman–Werner background is less efficient (O'Shea & Norman 2008; Susa 2008; Wise & Abel 2008; Wise & Cen 2009). However, 10⁸–10⁹ M_☉ halos are not immune to the Lyman–Werner background, just self-shielded at their very centers. Not all the baryons will be protected from external photodissociating flux, and this will affect H₂ production on dust outside the center of halo. In the lower mass halos, the Lyman–Werner background may be effective to dissociate H₂ molecules (Machacek et al. 2001, 2003; Yoshida et al. 2003; Susa 2007; Wise & Abel 2007; O'Shea & Norman 2008).

We first show the evolution of a galaxy with M_vir = 10⁹ M_☉ and z_vir = 10 for various dust destruction models. To clarify the dust destruction effects on galaxy evolution, we first show the result of the model without reverse and forward shocks (model C1m9), then compare the results of the models with only forward shocks (model B1m9) and with both forward and reverse shocks (model A1m9) to the model without both shocks (model C1m9).

In Figure 5, we show the evolution of various quantities without dust destruction (model C1m9). The figure shows the time evolution of the molecular fraction, $f_{{\rm H}_{2}}$; the SFR in units of M_☉ yr⁻¹, Ψ; the stellar mass fraction, M_star/(M_gas + M_star); the metallicity, Z, in units of Z_☉; total dust cross section per unit volume, σ_{d, −20}, in units of 10⁻²⁰ cm⁻¹; and the dust-to-gas mass ratio, D_{d, −2}, in units of 10⁻². The definition of D_{d, −2} is convenient for comparison with the MW value of the dust-to-gas mass ratio. In the MW, the dust-to-gas mass ratio is 0.5 × 10⁻² in the diffuse ISM (Draine 2009) and 0.9 × 10⁻² in molecular clouds (Pollack et al. 1994).

Zoom In Zoom Out Reset image size

The molecular fraction, $f_{{\rm H}_2}$, is very important, since it determines SFR and controls galaxy evolution. The molecular fraction reaches $f_{{\rm H}_{2}}\sim 1\times 10^{-3}$ around t ∼ 10⁷ yr. This result is robust for all models, since in this stage H₂ formation in the gas phase is dominant over that on dust grains (Tegmark et al. 1997; Hirashita & Ferrara 2002). The gas temperature rapidly drops below 200 K before 10⁷ yr. Then, the molecular fraction rapidly increases from t ∼ 10⁸ yr and reaches ∼0.83 at the galaxy age of ∼0.8 Gyr (z = 5). This is due to the enhancement of H₂ formation on dust grains by an increase of σ_{d, −20}. For t ≳ 3 × 10⁷ yr, σ_{d, −20} ≳ 1.1 × 10⁻⁴ and $[{d}f_{{\rm H}_{2}}/{d}t]_{\rm dust}$ exceeds $[{d}f_{{\rm H}_{2}}/{d}t]_{\rm star}$. For σ_{d, −20} ≳ 0.001, the increase of the molecular fraction enhances the star formation. The cycle of the H₂ formation on dust, the star formation, and the dust formation by SNe significantly accelerates galaxy evolution by, for example, rapidly increasing the stellar mass fraction, M_star/(M_gas + M_star). At the galaxy age ∼0.8 × 10⁹ yr, the stellar mass fraction goes up to M_star/(M_gas + M_star) ∼ 0.60. The SFR, Ψ(t), decreases from the time when M_star/(M_gas + M_star) ∼ 0.45, since gas mass decreases significantly. The active star formation causes the formation of dust grains and metals. At t ∼ 0.8 Gyr, the total dust cross section, σ_{d, −20}, goes up to 2.3 and the metallicity, Z, goes up to 3.3 × 10⁻¹ Z_☉.

Figure 6 shows the results of the galaxy model with the dust destruction by only the forward shocks (model B1m9) to illustrate the effects of dust destruction by forward shocks on the galaxy evolution. In this model n_SN = 1.0 cm⁻³. The dust destruction by forward shocks slightly affects the dust-to-gas mass ratio, D_{d, −2}, after the galaxy age of ∼5 × 10⁸ yr. This is because the destruction by forward shocks is roughly proportional to the dust-to-gas mass ratio (see Equation (9)), and dust grains are destroyed significantly for D_{d, −2} ≳ 0.1 in this case. The SFR decreases from the time when M_star/(M_gas + M_star) ∼ 0.4. At ∼0.8 Gyr, the molecular fraction reaches $f_{{\rm H}_{2}}\sim 0.51$ and the stellar mass fraction reaches M_star/(M_gas + M_star) ∼ 0.47.

Zoom In Zoom Out Reset image size

In Figure 7, we show the results of our fiducial model A1m9 in which we include the dust destruction by both reverse and forward shocks in the case of n_SN = 1 cm⁻³. The molecular fraction reaches $f_{{\rm H}_{2}}\sim 1\times 10^{-3}$ around t ∼ 10⁷ yr. This is similar to C1m9 and B1m9. After t ≳ 10⁷ yr, the molecular fraction evolution is quite different from models C1m9 and B1m9. The molecular fraction declines slowly until t ∼ 2 × 10⁸ yr. After ∼2 × 10⁸ yr, the molecular fraction increases slowly with increase of the dust mass. This is due to the H₂ formation on the dust grains. The molecular fraction reaches only ∼2.0 × 10⁻³ at t ∼ 0.8 Gyr. This is because dust destruction by reverse shocks is very effective and hence results in suppression of H₂ formation on dust grains. On the other hand, forward shocks hardly affect the evolution of dust size and dust mass, since the destruction of forward shocks can change dust mass only for a large dust-to-gas mass ratio, D_{d, −2} ≳ 10⁻¹. At t ∼ 0.8 Gyr, the stellar mass fraction reaches only M_star/(M_gas + M_star) ∼ 4.5 × 10⁻³, which is much less than the model without reverse shocks shown in Figure 5 (model C1m9) and 6 (model B1m9).

Zoom In Zoom Out Reset image size

We illustrate the difference in the H₂ formation rate among models C1m9, B1m9, and A1m9 as follows: The H₂ formation rate depends not only on the total dust mass but also on the dust size distribution. In models C1m9 and B1m9, the dust mass produced by an SN II without reverse shock is $\sum _{j}\int \overline{{\cal M}^{0}_{d,j}}(a){d}a=0.48\, M_{\odot }$ and in A1m9, the dust mass injection into ISM through a reverse shock with n_SN = 1.0 cm⁻³ is $\sum _{j}\int \overline{{\cal M}^{1.0}_{d,j}}(a){d}a=0.15\, M_{\odot }$. The ratio of the mean dust area to the mean dust volume is 〈a²〉/〈a³〉 = 1.5 × 10⁵ cm⁻¹ before the reverse shock destruction. After the reverse shock destruction, 〈a²〉/〈a³〉 = 4.2 × 10⁴ cm⁻¹. This is a measure of dust area per dust volume and is also a measure of the H₂ formation rate of the dust surface. Small 〈a²〉/〈a³〉 leads to a low H₂ formation rate. This is the reason why model A1m9 shows a smaller H₂ fraction than B1m9 and C1m9. The dust destruction by reverse shocks changes not only the dust mass but also the grain size distribution, and as a result drastically suppresses star formation in the galaxy.

In Figure 8, we show the evolution of the H₂ formation rate in gas phase, $[{d}f_{{\rm H}_2}/{d}t]_{\rm gas}$; the H₂ formation rate on dust grain, $[{d}f_{{\rm H}_2}/{d}t]_{\rm dust}$; the H₂ destruction rate by UV photons, $[{d}f_{{\rm H}_2}/{d}t]_{\rm UV}$; and the H₂ decreasing rate by the star formation, $[{d}f_{{\rm H}_2}/{d}t]_{\rm star}$, normalized to the total formation rate, $[{d}f_{{\rm H}_2}/{d}t]_{\rm gas}+[{d}f_{{\rm H}_2}/{d}t]_{\rm dust}$, in model A1m9. At t ∼ 4 × 10⁷ yr, the H₂ formation rate on dust grains exceeds the rate in the gas phase. However, the H₂ formation rate on dust grains is less than the H₂ decreasing rate by star formation at this epoch. At t ∼ 1.6 × 10⁸ yr, the H₂ formation rate on dust grains exceeds the H₂ decreasing rate by star formation. The formation on dust grains becomes the dominant process among all of H₂ formation and destruction processes. From this time when σ_{d, −20} ≳ 5.2 × 10⁻⁵, molecular fraction, $f_{{\rm H}_{2}}$, starts to increase. The molecular destruction rate by UV photons does not exceed the molecular decreasing rate by the star formation after t ∼ 4 × 10⁷ yr. After t ∼ 4 × 10⁷ yr, I_UV(3.1 × 10¹⁵ Hz) ∼ 2.0–4.8 × 10⁻²⁰ erg s⁻¹ cm⁻² Hz⁻¹ ster⁻¹, which corresponds to J₂₁ = 20–40 where J₂₁ is in units of I_UV(3.1 × 10¹⁵ Hz) = J₂₁ × 10⁻²¹ erg s⁻¹ cm⁻² Hz⁻¹ ster⁻¹. Note that J₂₁ in this model is higher than the values, J₂₁ ∼ 1, for the Lyman–Werner background in the redshift of 5 < z < 10 suggested in recent papers (e.g., Greif & Bromm 2006). The destruction process in the gas phase does not affect the evolution of molecular fraction significantly after t ≳ 1 × 10⁷ yr.

Zoom In Zoom Out Reset image size

The dependence of the time evolution of SFR on the ISM density around an SN is presented in Figure 9. We can see that after t ∼ 5 × 10⁷ yr higher density around SN progenitors results in lower molecular fraction and hence lower star formation efficiency. The SFR, Ψ(t), is independent of n_SN before t ∼ 5 × 10⁷ yr, because H₂ forms predominantly in the gas phase. In the model without reverse shock destruction, the SFR increases from ∼10⁸ yr and saturates around 5 × 10⁸ yr. This is because the gas is consumed by the star formation. In model A10m9 (n_SN = 10 cm⁻³), SFR is suppressed until t ∼ 0.8 × 10⁹ yr.

Zoom In Zoom Out Reset image size

We should note that it is probably true that n_SN < 1 cm⁻³ for Population III stars in the mass range of 20–40 M_☉, since Population III stars are massive and can photoevaporate the clouds in which they form, and ionized flows evacuate the dense gas around the stars to well below n_SN = 1 cm⁻³ (Whalen et al. 2004; Kitayama et al. 2004). In this case, consideration of circumstellar densities of 10 cm⁻³ and greater is not relevant to dust evolution in the SNR. In this paper, we consider n_SN > 5 cm⁻³ for completeness. If the stars are forming at lower redshift and are enriched, they will have stellar winds that also sweep away circumstellar gas to low densities.

We note that in a usual star formation recipe in both numerical simulations and analytic models, the SFR is assumed to increase with the local gas density, and our results show that the SFR is strongly affected by n_SN. We will discuss the effects of n_SN on the SFR in more detail in Section 5.

In Figure 10, we show the change of the molecular fraction, $f_{{\rm H}_{2}}$, with metallicity, Z, for various n_SN. In usual chemical evolution models, Z is a key indicator of galaxy evolution. However, as shown in this figure, $f_{{\rm H}_{2}}$ does not solely depend on the metallicity. For Z ≳ 5 × 10⁻⁴ Z_☉, $f_{{\rm H}_{2}}$ is large in models with small n_SN. This is because D_area/M_metal is large (small) in models with small (large) n_SN for the same Z as shown in Figure 3.

Zoom In Zoom Out Reset image size

The molecular fraction is well described by the total dust cross section per unit volume for σ_{d, −20} ≳ 0.001. In Figure 11, we show the evolution of the molecular fraction in terms of total dust cross section per unit volume. For σ_{d, −20} ≳ 1 × 10⁻³, H₂ formation on dust grains dominates $f_{{\rm H}_{2}}$ evolution as shown by the convergence of all models.

Zoom In Zoom Out Reset image size

Finally, we show the stellar mass fraction, M_star/(M_gas + M_star) for M_vir = 10⁸, 10⁹, 10¹⁰, and 10¹¹ M_☉ at z_vir = 10 in Figure 12. We note that dark halos of virial masses, M_vir = 10⁸, 10⁹, 10¹⁰, and 10¹¹ M_☉ correspond to the density fluctuation of 2.0σ, 2.5σ, 3.0σ, and 4.1σ, respectively. M_star/(M_gas + M_star) is large for large M_vir. This is explained as follows: the gas cools to the CMB temperature ∼30 K in all M_vir after t ∼ a few × 10⁷ yr, although T_vir increases with M_vir, so that the final gas density becomes higher because of smaller H/r_disk in larger M_vir (see Equation (15)). This results in more rapid molecular formation in a larger M_vir halo (see Equation (23)). The rapid molecular formation enhances the star formation and as a result causes the large stellar mass fraction. We should note that in higher z_vir, M_star/(M_gas + M_star) is larger for the galaxies with T_vir ≳ 10⁴ K, since the rotation timescale becomes smaller in our model.

Zoom In Zoom Out Reset image size