METALLICITY-DEPENDENT QUENCHING OF STAR FORMATION AT HIGH REDSHIFT IN SMALL GALAXIES

Following the "bath-tub" model of Bouché et al. (2010), we assume that the gas and stellar mass in a halo are determined by three processes: inflow, star formation, and ejection of material by star formation feedback. The continuity equations that govern this system are

$$\begin{equation} \dot{M}_{g} = \dot{M}_{g,\rm in} - (1-R) \dot{M}_{*,\rm form} - \epsilon _{\rm out} \dot{M}_{*,\rm form} \end{equation} \tag{ 1 }$$

$$\begin{equation} \dot{M}_* = \dot{M}_{*,\rm in} + (1-R) \dot{M}_{*,\rm form}. \end{equation} \tag{ 2 }$$

Here $\dot{M}_{g,\rm in}$ and $\dot{M}_{*,\rm in}$ are the rates of cosmological inflow of gas and stars (Section 2.1.1), $\dot{M}_{*,\rm form}$ is the SFR (Section 2.1.2), R is the stellar return fraction, and _out is the amount of mass ejected into the intergalactic medium (IGM) per unit mass of stars formed (Appendix A). We summarize the parameters that appear in these equations, as well as others that appear below, in Table 1. Below, we discuss our estimates for the quantities appearing in these equations.

Table 1. Model Parameters

Parameter	Meaning	Fiducial	$f_{\rm H_2} = 1$a	fg, in = 1b	Mret = 0.03c	ζlo = 0.8d	ZIGM = 2 × 10−4e	ΣSF = 10f
$f_{\rm H_2}$	H2 mass fraction	Equation (18)	1	Equation (18)	Equation (18)	Equation (18)	Equation (18)	Equation (24)
fg, in	Inflow gas fraction	Equation (B4)	Equation (B4)	1	Equation (B4)	Equation (B4)	Equation (B4)	Equation (B4)
Mret	Halo mass for metal retention	0.3	0.3	0.3	0.03	0.3	0.3	0.3
ζlo	Small halo metal ejection fraction	0.9	0.9	0.9	0.9	0.8	0.9	0.9
ζ	Overall Z ejection fraction	Equation (22)	Equation (22)	Equation (22)	Equation (22)	Equation (22)	Equation (22)	Equation (22)
in	Inflow fraction reaching disk	Equation (9)	Equation (9)	Equation (9)	Equation (9)	Equation (9)	Equation (9)	Equation (9)
fb	Universal baryon fraction	0.17	0.17	0.17	0.17	0.17	0.17	0.17
out	Star formation gas expulsion multiplier	1.0	1.0	1.0	1.0	1.0	1.0	1.0
R	Stellar return fraction	0.46	0.46	0.46	0.46	0.46	0.46	0.46
λ	Galaxy spin parameter	0.07	0.07	0.07	0.07	0.07	0.07	0.07
y	Stellar metal yield	0.069	0.069	0.069	0.069	0.069	0.069	0.069
ZIGM	IGM metallicity	2 × 10−5	2 × 10−5	2 × 10−5	2 × 10−5	2 × 10−5	2 × 10−4	2 × 10−5

Notes. ^aMetallicity-independent star formation model. ^bModel where cosmological inflow is purely gas, no stars. ^cModel where halos begin to retain most of their metals at a mass of 3 × 10¹⁰ M_☉, instead of the fiducial 3 × 10¹¹ M_☉. ^dModel where small halos eject 80% of the metals they produce, rather than the fiducial 90%. ^eModel where the IGM metallicity is 10 times that in the fiducial model. ^fModel where we allow star formation everywhere above a metallicity-independent threshold surface density Σ_SF = 10 M_☉ pc⁻².

Download table as:  ASCIITypeset image

2.1.1. Accretion Rates

The average total (baryonic plus dark matter) IR into a halo of virial mass M_h at redshift z due to cosmological instreaming from the cosmic web is approximated by

$$\begin{equation} \dot{M}_{h,12} = -\alpha M_{h,12}^{1+\beta } \dot{\omega }, \end{equation} \tag{ 3 }$$

where M_{h, 12} = M_h/10¹² M_☉, α = 0.628, β = 0.14. Here, ω = 1.68/D(t) is the self-similar time variable of the extended Press-Schechter (EPS) formalism, in which D(t) is the linear fluctuation growing mode, and it is approximated by

$$\begin{equation} \dot{\omega } = -0.0476[1+z+0.093(1+z)^{-1.22}]^{2.5} \mbox{ Gyr}^{-1}. \end{equation} \tag{ 4 }$$

This approximation has been derived by Neistein et al. (2006) and Neistein & Dekel (2008) based on EPS theory and was fine-tuned using the Millennium cosmological N-body simulation. We adopt the cosmological parameters from the Wilkinson Microwave Anisotropy Probe: Ω_m = 0.27, $\Omega _\Lambda = 0.73$, h = 0.7, and σ₈ = 0.81.⁵ This expression is accurate to better than 5% for z = 0.2–5, and to ∼10% for z = 0–10. The mass dependence is accurate to 5% for halos from mass M_{h, 12} = 0.1–10². The distribution of accretion rates about this mean can be fit by a lognormal with a standard deviation of ∼0.2, plus a tail representing mergers that extends to a factor of 10 above the average (Dekel et al. 2009; T. Goerdt et al. 2012, in preparation), though we do not implement this scatter in our model. Although we use Equation (3) for all our numerical computations, a simpler expression, which we will use to make scaling arguments later, is

$$\begin{equation} \dot{M}_h \approx 34\,M_\odot \mbox{ yr}^{-1} M_{h,12}^{1.14} (1+z)^{2.4}. \end{equation} \tag{ 5 }$$

This matches Equation (3) to ∼10%. We note that if the power 1.14 is approximated by unity, and if the 2.4 is replaced by 5/2 (which is the exact predicted value in the Einstein–de Sitter cosmology phase that is approximately valid at z > 1), this accretion rate corresponds to a mass growth rate M_h∝exp (− γz), with γ ≃ 0.9.

The corresponding gas and stellar IRs into a halo are

$$\begin{equation} \dot{M}_{g,\rm in} = 0.17 \epsilon _{\rm in} f_{b,0.17} f_{g,\rm in} \dot{M}_{h} \end{equation} \tag{ 6 }$$

$$\begin{equation} \dot{M}_{*,\rm in} = 0.17f_{b,0.17} (1-f_{g,\rm in}) \dot{M}_{h}, \end{equation} \tag{ 7 }$$

where f_b is the cosmic baryon fraction, f_{b, 0.17} = f_{b, in}/0.17, f_{g, in} is the fraction of the inflow that reaches the halo as gas rather than stars, and _in is the fraction of accreted gas that reaches the galactic disk rather than being shock-heated and going into the galactic halo. Note that we have not attempted to include the effects of a reduction in accretion rates onto small halos with virial velocities due to suppression of cooling by the UV background after the universe is reionized. This is expected to be a large effect for halos with virial velocities below 30 km s⁻¹, and a factor of a few effect for halos with virial velocities of 30–50 km s⁻¹ (Thoul & Weinberg 1996). Our omission of the effect will make some difference for the very smallest halos that we model, but the majority of the models we present below are in the regime where suppression of accretion is unimportant or is only a factor of a few effect.

In our fiducial model, we compute f_{g, in} using a self-consistent approximation described in Appendix B, and we explore how this affects our results by also considering a variant model in which we take f_{g, in} = 1 for all halos at all redshifts.

We have estimated the fraction of gas inflow from outside the virial radius that actually reaches the galaxy at the halo center, _in, using analytic arguments and hydro cosmological simulations (Birnboim & Dekel 2003; Kereš et al. 2005, 2009b; Dekel & Birnboim 2006; Ocvirk et al. 2008; Dekel et al. 2009). Well below a critical halo mass M_{h, 12} ∼ 1 there is no virial shock, and the gas flows in cold, so _in ≃ 1. At z > 2, in halos above the critical mass, cold streams bring in most of the gas along the dark matter filaments of the cosmic web, and they penetrate efficiently through the hot medium, yielding _in ∼ 1 also there (Dekel et al. 2009). A statistical analysis of 400 adaptive mesh refinement simulated halos at z ∼ 2.5 indeed reveals _in ∼ 0.9, with the exact value depending on how the averaging is performed (Danovich et al. 2012; A. Dekel et al. 2012, in preparation). Another estimate based on smoothed particle hydrodynamics simulations (Faucher-Giguere et al. 2011) gives similar values at high-z but somewhat smaller values, _in ∼ 0.5, at z ∼ 2.5. At z < 2, in halos more massive than the shock-heating scale, the cold flows are broader and their penetration power is weaker, so the accretion rate into the disk is suppressed (Dekel & Birnboim 2006; Ocvirk et al. 2008). This turns out to explain the evolution of the red sequence of galaxies (Dekel & Birnboim 2006; Cattaneo et al. 2006). We crudely model this using a simple analytic form following Dekel et al. (2009):

$$\begin{equation} M_{\rm 12, max} = \max \left(2,\frac{1}{3 M_{*,\rm PS}(z)}\right) \end{equation} \tag{ 8 }$$

$$\begin{equation} \epsilon _{\rm in} = \left\lbrace \begin{array}{@{}ll}1, \quad & M_{h,12} < M_{\rm 12, max} \\ \ms 0, & M_{h,12} > M_{\rm 12,max} \end{array} \right.. \end{equation} \tag{ 9 }$$

Here, M_{*, PS}(z) is the Press–Schechter mass at a given redshift, which we compute from the Sheth & Tormen (1999) ellipsoidal collapse model following the procedure outlined in Mo & White (2002). We have also run models in which we adopt the analytic fitting formulae proposed by Faucher-Giguere et al. (2011), and found that the results are qualitatively unchanged. Of course the use of a sharp mass cutoff for accretion is an oversimplification of reality. Clearly some galaxies with masses below M_{12, max} are early types without accretion or star formation, while some with larger masses are late types that have ongoing star formation. Indeed, based on cosmological simulations, the transition from a population of galaxies dominated by full cosmological accretion rate to a population dominated by shutdown is stretched over more than a decade in halo mass (e.g., Figure 1 of Birnboim et al. 2007; Ocvirk et al. 2008; Kereš et al. 2009a). However, given the simplicity of the remainder of our model, there would be little point in attempting to include this scatter. Moreover, given the difficulty in carrying high-resolution simulations of cosmological inflow to redshifts z ≲ 1, theoretical estimates of the accretion rates there must be considered highly uncertain in any event. In this paper, we will be more concerned with redshifts z ≳ 2, where Equation (9) may be regarded as reasonably accurate.

2.1.2. Metallicity-dependent Star Formation

The model presented in the previous section is similar to that of Bouché et al. (2010). We now diverge from that model by taking into account the atomic–molecular phase transition in computing the rate at which gas turns into stars. On scales ranging from individual clouds to entire galaxies, the SFR is well described by $\dot{M}_* \sim \epsilon _{\rm ff}M_{\rm H_2} / t_{\rm ff}$, where $M_{\rm H_2}$ is the total molecular gas mass, t_ff is the free-fall time of the gas, and _ff ∼ 0.01 (Krumholz & Tan 2007; Krumholz & Thompson 2007). This result is an integrated version of the classical Kennicutt (1998) relation based on observations.

However, galaxies at low surface densities and metallicities are observed to form stars at a rate considerably below that predicted by the Kennicutt relation (e.g., Wolfe & Chen 2006; Wild et al. 2007; Leroy et al. 2008; Bigiel et al. 2008; Wyder et al. 2009), a result that has been successfully explained by noting that stars form only in molecular gas, and that systems with low surface density and metallicity tend to have little of their gas in molecular form (Robertson & Kravtsov 2008; Krumholz et al. 2009a, 2009c; Gnedin et al. 2009; Gnedin & Kravtsov 2010). We therefore adopt an SFR in a given halo:

$$\begin{equation} \dot{M}_{*,\rm form} = 2\pi \int _0^{\infty } f_{\rm H_2} \frac{\epsilon _{\rm ff}}{t_{\rm ff}}\Sigma _g r \, dr, \end{equation} \tag{ 10 }$$

where Σ_g is the gas surface density and $f_{\rm H_2}$ is the fraction of that gas in molecular form. All quantities inside the integral are functions of the distance from the disk center.

For the surface density, we assume that the gas in disks follows an exponential profile $\Sigma =\Sigma _c e^{-r/R_d}$. The scale length R_d is assumed to be proportional to the virial radius of the halo in which it resides,

$$\begin{equation} R_d = 0.05 \lambda _{0.1} R_v, \end{equation} \tag{ 11 }$$

where $\lambda = (1/\sqrt{2})(J/M_h)/({\it VR})$ is the spin parameter for a halo of angular momentum J, mass M_h, virial radius R_v, and circular velocity V (Bullock et al. 2001), and λ_0.1 = λ/0.1. Note that we choose to make R_d half of λR_v because for an exponential disk the half-mass radius is 1.7 scale radii, and given the uncertainties our ansatz places the half-mass radius at roughly λR_v. Tidal torque theory and N-body simulations give on average λ ≃ 0.04, but the observed radii of z ∼ 2 galaxies (Genzel et al. 2006, 2008) suggest that the gas in them has a somewhat higher angular momentum. Following Dutton et al. (2011), we adopt λ ≃ 0.07 as our typical value. The virial radius is related to the virial mass and the expansion factor a = 1/(1 + z) by (Dekel & Birnboim 2006):

$$\begin{equation} R_{v,100} \simeq 1.03 M_{h,12}^{1/3} A_{1/3}, \end{equation} \tag{ 12 }$$

where R_{v, 100} = R_v/100 kpc, A_1/3 = A/(1/3),

$$\begin{equation} A = \left(\Delta _{200} \Omega _{m,0.3} h_{0.7}^2\right)^{-1/3} a, \end{equation} \tag{ 13 }$$

and Ω_{m, 0.3} = Ω_m/0.3. Here,

$$\begin{equation} \Delta (a) \simeq (18\pi ^2 - 82 \Omega _{\Lambda }(a) - 39 \Omega _{\Lambda }(a)^2)/\Omega _m(a) \end{equation} \tag{ 14 }$$

is the overdensity within the virial radius at a given epoch, Δ₂₀₀ = Δ/200, $\Omega _m(a) = \Omega _m a^{-3}/(\Omega _\Lambda +\Omega _m a^{-3})$, and $\Omega _\Lambda (a) = 1-\Omega _m(a)$. At redshifts ≳ 1, in the Einstein–de Sitter regime, Δ₂₀₀ ≈ 1 and so A ≈ a for our standard cosmology. In this cosmology, Δ₂₀₀ grows to 1.7 at redshift 0. Given this surface density profile, the central surface density is related to the disk gas mass by

$$\begin{equation} \Sigma _{c,0} = \frac{M_g}{2\pi R_d^2} = 0.125 M_{g,11} \lambda _{0.1}^{-2} M_{h,12}^{-2/3} A_{1/3}^{-2}, \end{equation} \tag{ 15 }$$

where Σ_{c, 0} = Σ_c/(1 g cm⁻²). Note that, if M_g∝M_h, the surface density scales with mass and redshift like ∝M^1/3(1 + z)².

Also note that our treatment of disk sizes is extremely simple. We have ignored effects like the variation of halo concentration with mass and redshift, and the possible evolution of galactic angular momenta with redshift or other galaxy properties (e.g., Burkert et al. 2010). The disk size affects our models mainly by changing the surface density and thus the radius at which gas transitions from atomic to molecular, as discussed in the next section. This transition is also affected by metallicity, and the uncertainties associated with metallicity evolution are probably at the order of magnitude level. Given the size of these uncertainties, a more detailed treatment of disk radii seems unnecessary.

The local free-fall time depends on the gas surface density Σ of the galactic disk. Krumholz et al. (2009c) have analyzed the structure of molecular clouds and developed a theoretical model for star formation that gives

$$\begin{equation} \frac{\epsilon _{\rm ff}}{t_{\rm ff}} \approx \left(2.6\mbox{ Gyr}\right)^{-1} \left\lbrace \begin{array}{@{}ll}(\Sigma _0/0.18)^{-0.33}, \quad & \Sigma _0 < 0.18 \\ (\Sigma _0/0.18)^{0.33}, & \Sigma _0 \ge 0.18 \end{array} \right. \!\!, \end{equation} \tag{ 16 }$$

where Σ₀ = Σ/(1 g cm⁻²). The low Σ regime corresponds to galaxies like the present-day Milky Way where molecular clouds are confined by self-gravity, while the high Σ regime describes galaxies like low-redshift ultraluminous infrared galaxies where they are confined by external pressure. This formula agrees well with the SFR observed in nearby galaxies (cf. Figures 1 and 2 of Krumholz et al. 2009c), and we adopt it here.

Most importantly, the fraction of the ISM in the cold, molecular phase, $f_{\rm H_2}$, is determined by the balance between grain photoelectric heating and UV photodissociation on one hand, and collisionally excited metal line cooling and H₂ formation on dust grains on the other hand. (Near zero metallicity other processes become important, but, as we discuss below, we will not consider this regime.) Krumholz et al. (2008, 2009b) and McKee & Krumholz (2010) show that these processes produce a molecular fraction that depends primarily on the surface density and metallicity of the galactic disk, and relatively little on any other parameters. A crude approximation to their result is that the molecular fraction is

$$\begin{equation} f_{\rm H_2} \sim \Sigma / \big(\Sigma + 10 \,Z_0^{-1}\, M_\odot \,\mbox{pc}^{-2} \big), \end{equation} \tag{ 17 }$$

where Σ is the total gas column density and Z₀ = (M_Z/M_g)/Z_☉ is the metallicity normalized to the solar neighborhood value Z_☉ = 0.02. Thus, regions with Σ ≪ 10 Z⁻¹₀ M_☉ pc⁻² are primarily atomic, and those with Σ ≫ 10 Z⁻¹₀ M_☉ pc⁻² are primarily molecular. A more accurate expression, which we adopt here, is

$$\begin{equation} f_{\rm H_2} = \left\lbrace \begin{array}{@{}ll}1 - \frac{3}{4}\left(\frac{s}{1+0.25s}\right), \quad & s < 2 \\ \ms 0, & s\ge 2\\ \end{array} \right. \end{equation} \tag{ 18 }$$

$$\begin{equation} s = \frac{\ln (1+0.6\chi +0.01\chi ^2)}{0.6 \tau _c} \end{equation} \tag{ 19 }$$

$$\begin{equation} \chi = 3.1 \frac{1+Z_0^{0.365}}{4.1} \end{equation} \tag{ 20 }$$

$$\begin{equation} \tau _c = 320 c Z_0 \Sigma _0, \end{equation} \tag{ 21 }$$

where c is a clumping factor that accounts for smoothing of the surface density on scales larger than that of a single atomic–molecular complex. If Σ is measured on 100 pc scales then there is no averaging and c ≈ 1, while on ∼1 kpc scales it is ∼5 (Krumholz et al. 2009c), and we adopt this as a fiducial parameter in our models. This analytic model agrees very well with numerical simulations that follow the full time-dependent chemistry of H₂ formation and dissociation (Krumholz & Gnedin 2011).

It is also possible to form stars in galaxies that are essentially free of molecules and dust, a topic that has been studied extensively since the pioneering work of Bromm et al. (2002) and Abel et al. (2002). In these cases the cooling required for star formation is driven by the tiny fraction of H₂ that is able to form by gas-phase processes, rather than via grain catalysis, the dominant process over most of cosmic time. However, this Population III process is extremely slow and inefficient compared to the normal mode of star formation. While it is important for providing the seed metals that enable the normal mode to begin, we will neglect the contribution of Population III stars to the total SFR of the universe.

The final ingredient to our model is the metallicity of the gas, which affects its ability to form a cold molecular phase. The first seed metals that begin the process of H₂ formation and the transition to normal star formation are produced by these Population III stars. A single pair-instability supernova from one of these stars pollutes its host halo up to a metallicity of ∼10⁻³ solar (Wise et al. 2012). We therefore adopt Z_IGM = 2 × 10⁻⁵ as the fiducial "IGM" metallicity at which all halos start, and which characterizes newly accreted gas. We test our sensitivity to this choice below by also considering a model with Z_IGM 10 times higher.

Once the molecule fraction becomes non-zero and normal star formation begins, we set the metallicity based on the expected yield of the resulting stars. Following the standard practice in chemical evolution models (e.g., Maeder 1992), we use the instantaneous recycling approximation and write the metal production rate in terms of the yield y of the stellar population. The net metal production rate corresponding to an SFR $\dot{M}_{\rm *,form}$ is $y(1-R) \dot{M}_{\rm *,form}$; we adopt y = 0.069, as discussed in Appendix A. However, only a fraction of those metals will be retained in the galaxy rather than lost in supernova explosions. This fraction is a very sharply increasing function of halo mass, and based on a rough fit to the simulations of Mac Low & Ferrara (1999), we take the fraction of metals ejected by supernovae to be

$$\begin{equation} \zeta = \zeta _{\rm lo} e^{-M_{h,12}/M_{\rm ret}}, \end{equation} \tag{ 22 }$$

where M_ret is the halo mass (in units of 10¹² M_☉) at which supernovae become unable to eject most of their metals, and ζ_lo is the fraction of metals retained in halos much smaller than M_ret. As a fiducial value, based on the simulations of Mac Low & Ferrara (1999), we adopt M_ret = 0.3. Choosing a fiducial value of ζ_lo is more difficult, because even if almost all metals are initially ejected from the galactic disk in low-mass halos (as Mac Low & Ferrara find), many of these will not escape the halo entirely, and will later be re-accreted—as must be the case, since halos below M_{h, 12} = 0.3 are not completely devoid of metals. In the absence of firm theoretical guidance we adopt as a fiducial value ζ_lo = 0.9, i.e., we assume that 90% of metals will be ejected completely from small halos, while 10% will be retained or re-accreted. Below, we test the sensitivity of our results to both M_ret and ζ_lo by varying these fiducial values.

The two other factors capable of changing the metal content are accretion of pristine material, which adds metals at a rate $Z_{\rm IGM} \dot{M}_{g,\rm in}$, and expulsion of existing ISM by feedback, which removes metals at a rate $(M_Z/M_g) \epsilon _{\rm out} \dot{M}_{\rm *,form}$. This expression assumes that the expelled gas has a metallicity equal to the mean metallicity of the galaxy. Combining metal production, accretion of pristine material, and ejection of existing interstellar material, the evolution equation for the mass of metals is

$$\begin{eqnarray} \dot{M}_Z & = & y (1-R) (1-\zeta) \dot{M}_{\rm *,form} \nonumber \\ & & {} + Z_{\rm IGM} \dot{M}_{g,\rm in} - \epsilon _{\rm out} \frac{M_Z}{M_g} \dot{M}_{\rm *,form}. \end{eqnarray} \tag{ 23 }$$

Note that we distinguish between the ejection of hot metals that are already part of a 10,000 km s⁻¹ supernova shock and the driving of outflow in the cold dense ISM, the former being easier to eject and more strongly dependent on the halo mass via ζ. We have verified that our results are rather insensitive to the value of _out and to a possible halo-mass dependence in it; for more discussion of this issue, see Appendix A.

Finally, we caution that our treatment of metals neglects the existence of metallicity gradients within galactic disks. In the local universe, these are observed to be relatively modest. For spiral galaxies, the center-to-edge metallicity difference is typically at most half a dex (e.g., Oey & Kennicutt 1993), with the metallicity gradient disappearing completely outside R₂₅ (e.g., Werk et al. 2011). Dwarf galaxies also have essentially no metallicity gradient (e.g., Croxall et al. 2009). Thus, metallicity gradients seem unlikely to very significantly alter our results. Any potential effects are certainly likely to be smaller than those induced by changing Z_IGM by an order of magnitude, as we do below.

We consider four example halos, with present-day masses of M_{h, 12} = 0.008, 0.016, 1.0, and 10.0; these masses are chosen to illustrate several interesting behaviors in the model. We note that each of these halos has a virial velocity above 25 km s⁻¹ at all redshifts below z ∼ 10, and above 30 km s⁻¹ from z = 2–3, when the bulk of the gas is accreted and most star formation occurs. Thus, these halos will not be significantly affected by ionization-driven winds (Shaviv & Dekel 2003), though for the smallest ones the accretion rates may be reduced by a factor of a few due to reduction in cooling by the UV background (Efstathiou 1992; Thoul & Weinberg 1996; Hoeft et al. 2006). Note, however, that below we will plot results for smaller halos for which photoionization effects are likely be important. As noted above, we have not attempted to include these simply to avoid further complicating the model. Figures 1–4 summarize the evolutionary history of our four example halos using the fiducial model. Figure 1 shows the mass of each component, Figures 2 and 3 show the total and specific rates of star formation (including only quiescent star formation, not stars formed in accretion-driven bursts), as well as gas and stellar inflow, and Figure 4 shows the evolution of the gas-phase metallicity and column density.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1.1. Star Formation Cannot Catch up at High Redshift

We see that, at very high redshift (z ≫ 5) all halos are dominated by gas rather than stellar mass, and they are forming stars more slowly than they are accreting both gas and stars from the IGM. However, there is star formation during this phase, which is critical for building up metals in the ISM. Although the halos' starting metallicities are quite low, this is compensated for by their high column densities, and as a result much of their mass is able to form molecules and participate in star formation. They remain gas dominated simply because their SFR is unable to keep up with the IR. To see this, note that the characteristic timescale for a halo to double its mass via inflow is (using the approximate Equation (5) for the mass IR)

$$\begin{equation} t_{\rm acc} = \frac{M_h}{\dot{M}_h} = 2.1\mbox{ Gyr } M_{h,12}^{0.14} (1+z)_3^{-2.4}, \end{equation} \tag{ 25 }$$

where (1 + z)₃ = (1 + z)/3. The timescale for gas inflow to change the gas mass is similar during the era when star formation has not yet significantly reduced the gas mass. In comparison, for a halo with M_g = 0.17f_{b, 0.17}f_gM_h, where f_g is the fraction of the baryonic mass in gas, we can estimate the time required to turn an order of unity of the molecular gas in a galaxy into stars by plugging M_g into Equation (15) and then Equation (16). This gives a timescale for star formation at the disk center:

$$\begin{equation} t_{\rm SF} \equiv \frac{t_{\rm ff}}{\epsilon _{\rm ff}} = 2.5\mbox{ Gyr } \lambda _{0.1}^{0.67} f_{b,0.17}^{-0.33} f_g^{-0.33} M_{h,12}^{-0.1} (1+z)_3^{-0.67}, \end{equation} \tag{ 26 }$$

where we have taken A ≈ a and we have assumed that we are in the Σ_{c, 0} > 0.18 regime; both of these approximations hold well at high-z. Thus, we see that, for z ≫ 2, t_acc ≪ t_SF, so we expect that the SFR will not be able to keep up with the IR, more so at higher redshifts. This effect is even further enhanced by the fact that only a portion of the gas is cold and molecular. We note that even if t_SF scaled with the disk and halo crossing times, which scale with the Hubble time and thus ∝(1 + z)^−1.5, the SFR was still unable to catch up with the IR at high enough redshift. We should emphasize that the inability of the SFR to catch up with the IR is a result of having the disk radius given by Equation (11). During a major merger the gas can be driven to a much smaller radius, leading to a shorter free-fall time and more rapid star formation. However, there are not enough major mergers to make a substantial difference in gas consumption on cosmological scales (Neistein & Dekel 2008; Dekel et al. 2009).

3.1.2. Metal Buildup and the Peak in SFR

3.1.3. Drop Beyond the Peak

3.2.1. The Effect of Metallicity on SFR History

In Figure 5, we compare the SFR in our fiducial model to the SFRs in our variant models. We see that the variants and the fiducial model are all very similar for the two higher mass halos. These are able to form molecules with relatively little difficulty due to their consistently high column densities, so their star formation histories are not significantly affected by metallicity-dependent star formation. The story is different for the lower mass halos. For them, making the accretion purely gaseous (f_{g, in} = 1) or reducing the metal retention mass (M_ret = 0.03) still has relatively modest effects—the star formation lasts slightly longer in the M_{h, 12} = 0.008 halo, and peaks slightly sooner in the M_{h, 12} = 0.016 one—but qualitatively there is little difference. In contrast, assuming that gas can form stars regardless of whether it is in a cold, molecular phase ($f_{\rm H_2} = 1$), and that this phase transition occurs at a metallicity-independent threshold of Σ_SF = 10 M_☉ pc⁻², decreasing the maximum fraction of metals from supernovae that are ejected (ζ_lo = 0.8), or increasing the IGM metallicity (Z_IGM = 2 × 10⁻⁴) has dramatic effects in these small halos. All these changes allow much more rapid star formation at high-z, reduce the star formation peak behavior at z ∼ 1–2 in the M_{h, 12} = 0.016 halo, and allow star formation to continue up to the present day in the M_{h, 12} = 0.008 halo. Thus, we see that including a realistic treatment of ISM chemistry and thermodynamics in models is likely to have a dramatic effect on the star formation history of the universe, and to make star formation sensitive to exactly how metals are distributed in halos, while changing the stellar inflow fraction or the rate of star formation in accretion shocks has relatively little effect.

Zoom In Zoom Out Reset image size

The insensitivity of our results to the choice of M_ret might at first seem somewhat surprising, given how dramatically metallicity-dependent star formation changes the halos' star formation history. However, note that the halos where metallicity makes a difference have masses much smaller than M_ret even if we reduce M_ret by a factor of 10 compared to Mac Low & Ferrara's value. This suggests that what matters in determining the mass below which star formation will be suppressed is not the precise mass at which halos start retaining most of their metals, but instead the ability of halos below this retention mass to build up a non-negligible level of metallicity, which in turn depends on the ratio of metal production rate to accretion rate in these halos. This is a function of the SFR, the yield, the metallicity of the accreting gas, and the fraction of metals that are retained in small halos.⁷

3.2.2. The Origin of an sSFR Plateau

In Figure 6, we show the sSFR for our fiducial model and for all the variants. Here, we see that setting f_{g, in} = 1 has a major effect on the results. It is easiest to understand this by focusing on the two higher mass halos, since for these both the fiducial model and the variants have about the same total SFR. We find that the fiducial model and all the variant models except f_{g, in} = 1 have sSFRs in the two higher mass halos that are relatively flat at high-z. This is fairly easy to understand. At redshifts above z ∼ 4–5, in all of these models the stellar accretion rate exceeds the in situ SFR (cf. Figures 2 and 3). By examining where in Figure 1 the stellar mass curve begins to deviate from simply tracking the total baryonic mass, one can see that stars formed in situ do not dominate the total stellar mass until z ∼ 3. Thus, at z ≳ 3, the stellar mass in a halo of mass M_h is approximately M_* = 0.17f_{b, 0.17}(1 − f_{g, in})M_h. Assuming the gas mass is M_g = 0.17f_{b, 0.17}f_{g, in}M_h (i.e., that gas has not been significantly depleted by star formation yet), we therefore have an sSFR

$$\begin{equation} \mbox{sSFR} = \frac{\dot{M}_{*,\rm form}}{M_*} = \frac{f_{g,\rm in}}{1-f_{g,\rm in}} \left\langle f_{\rm H_2} \frac{\epsilon _{\rm ff}}{t_{\rm ff}}\right\rangle, \end{equation} \tag{ 27 }$$

where the angle brackets indicate a surface-density-weighted average over the disk. Since the latter is generally ∼(2–3 Gyr)⁻¹, we obtain a nearly constant sSFR of 1–3 Gyr⁻¹ at high redshift. The emergence of an sSFR plateau at high redshifts is insensitive to the exact value of f_{g, in}, as long as the fraction of ex situ stars is non-negligible, but it affects the amplitude of the sSFR plateau, and its extent toward lower redshifts. Metallicity-dependent star formation flattens the sSFR as a function of z even further, because it makes high-z halos even more dominated by accreted/accretion shock-formed stars than in the model with $f_{\rm H_2} = 1$. In contrast, in the models with f_{g, in} = 1 halos are dominated by in situ stars at all epochs and masses. As a result, M_* is very small at high-z, giving rise to a larger sSFR. Since t_SF gets close to t_acc as we approach the present day, M_* rises faster than M_h, and the sSFR declines toward the present, rather than remaining flat.

Zoom In Zoom Out Reset image size

At intermediate redshifts, z ∼ 2–4, where the in situ SFR is already a major contributor to the stellar mass, the fact that the SFR is growing roughly exponentially with time leads to a similar exponential growth of its integral, the stellar mass, and thus helps extending the plateau toward z ∼ 2.

We now turn our attention to the ensemble properties of our entire model grid. Figure 7 shows the SFRs and gas accretion rates in the halos as they evolve in mass with redshift. To help interpret this plot, in Figure 8 we show two ratios:

$$\begin{equation} s_{\rm H_2} \equiv \frac{\dot{M}_{*,\rm fiducial}}{\dot{M}_{*,f_{\rm H_2}=1}} \end{equation} \tag{ 28 }$$

$$\begin{equation} s_{\rm IR} \equiv \frac{\dot{M}_{*,\rm fiducial}}{\dot{M}_{g,\rm in,fiducial}}. \end{equation} \tag{ 29 }$$

The first of these describes the factor by which the SFR is reduced by our inclusion of metallicity quenching relative to a model that does not include. The second describes the factor by which the SFR is suppressed relative to the rate of gas inflow in our fiducial model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

These plots show that a star formation law in which stars form only in cold molecular gas dramatically suppresses star formation in low-mass halos. Notice that contours of constant SFR and large $s_{\rm H_2}$ are close to horizontal at z ≳ 2, indicating that metallicity quenching acts much like a halo-mass threshold for star formation at z ≳ 2. The value of this threshold mass is sensitive to how well small halos are able to retain their metals and to the metallicity of accreting IGM gas, changing from ∼10¹⁰ M_☉ in our fiducial model to ∼10⁹M_☉ in our variant where small halos retain twice as much of their metal production, and to slightly smaller masses in the model where the IGM metallicity is raised to 1% of solar. As we will see below, this induces a ∼0.5 dex shift in the total SFR budget of the universe at high-z. In contrast, both the gas IR and the star formation in the model where stars form equally well in warm atomic or cold molecular gas ($f_{\rm H_2}=1$) or where the atomic–molecular transition is taken to be metallicity-independent (Σ_SF = 10 M_☉ pc⁻²) continue all the way down to the smallest halos in our model set. For halos above the threshold, the story is quite different. At z ≳ 4, metallicity quenching reduces the SFR in these halos too, but only by factors of $s_{\rm H_2}\sim 2$. In contrast, the SFR lags the gas IR by factors of s_IR ∼ 3–10, simply because star formation cannot keep up with the IR. Conversely, at z ≲ 4 we find that $s_{\rm H_2}<1$ in massive halos, indicating that SFRs are actually higher in the fiducial model than in the $f_{\rm H_2} = 1$ or Σ_SF = 10 M_☉ pc⁻² models. The reason for this behavior is that metallicity quenching reduces the SFR at high-z when the metallicity is low, but it does not remove the gas. At lower-z when the metallicity rises, this gas is able to form stars. Thus, metallicity quenching in large halos simply delays star formation, helping to produce the peak of cosmic SFR at z ∼ 2. In contrast, the choice of f_{g, in} or M_ret = 0.03 clearly makes little difference. We therefore conclude that the star formation history of the universe is sensitive to the need to form a cold phase, and that the ability of different halos to do so depends on their ability to retain metals. The star formation history of the universe does not depend strongly on stars being formed externally or in accretion shocks or on the exact mass at which halos become efficient at retaining metals.

Figure 9 shows the corresponding specific star formation and specific IRs. We see that the fiducial model exhibits a very wide range of halos masses and redshifts over which the sSFR is ∼2 Gyr⁻¹. At lower-z the sSFR starts to fall, as star formation both increases the stellar mass and decreases the gas mass, but this process takes some time to complete, and in most halos the sSFR does not drop below 1 Gyr⁻¹ until after z = 2, and does not fall below 0.5 Gyr⁻¹ until z ∼ 1. The model with $f_{\rm H_2}=1$ exhibits a slightly larger range of sSFRs, but except for very low mass halos, it is qualitatively similar to the fiducial model. Thus, we conclude that molecules do not dramatically change the sSFR in massive halos.

Zoom In Zoom Out Reset image size

In contrast, the model with f_{g, in} = 1 reaches sSFRs of many tens per Gyr, and the specific IR reaches similarly large values. This is the same effect we identified in the previous section in studying the histories of individual halos. When f_{g, in} ≠ 1, halos at very high-z tend to be dominated by stars that are directly accreted or formed in accretion shocks, so their stellar mass simply scales with the halo mass. Since the SFR scales with the gas mass (and thus also with the halo mass), the sSFR is constant. For f_{g, in} = 1, on the other hand, the stellar mass is no longer proportional to the halo mass at high-z, and the sSFR becomes large.

It is interesting to notice the contrasting roles played by metallicity quenching (i.e., the fact that $f_{\rm H_2} \ne 1$) and stellar inflow (i.e., the fact that f_{g, in} ≠ 1). In order to have an SFR that increases with time but an sSFR that is flat, as the observations discussed in Section 1 appear to demand, one requires that the SFR deviate from the IR, but that the total stellar mass does follow the time-integrated IR. These seemingly contradictory requirements are met by a combination of suppression of star formation by metallicity quenching and increase in stellar mass due to stellar inflow.

We can also use our models to predict the evolution of SFR or other galaxy properties with observational samples selected in a variety of ways. One observational approach is to select galaxies based on their comoving number density. In this strategy, at each redshift z one selects galaxies with luminosities L near a cutoff luminosity L(z) chosen based on the condition that the comoving number density of galaxies with luminosities L > L(z) satisfy the condition that n(> L(z)) = n₀, where n₀ is some fixed number density. This approach is useful because the comoving number density of main halo progenitors should, to good approximation, remain constant with redshift. Thus, this technique should produce a selection that is close to following one of our theoretical main-progenitor tracks in Figure 7. It is therefore a straightforward prediction of our model.

To compare our models to samples of this sort, we first compute the halo abundance n(M_h, z), the number density of halos at redshift z with masses in the range M_h to M_h + dM_h, using the Sheth & Tormen (1999) approximation, following the computation procedure outlined by Mo & White (2002). At each z, we then numerically solve for the halo mass M_h(z) for which $\int _{M_h(z)}^\infty n(M_h,z) \, dM_h = n_0$, where n₀ is the target number density. We then interpolate between the two model halos in our grid that most closely bracket M_h(z).⁸

In Figures 10 and 11, we compare the SFR and stellar mass as a function of look-back time for our model halos to three number density-selected samples, taken from Marchesini et al. (2009), Papovich et al. (2011), and B. Lundgren et al. (2012, in preparation). The observational samples use a threshold density n_{0, obs} = 2 × 10⁻⁴ Mpc⁻³ for Papovich et al. and Marchesini et al., and 1.8 × 10⁻⁴ Mpc⁻³ for Lundgren et al. For the latter sample, only stellar masses are available, not SFRs. The choice of n₀ for the theoretical comparison requires some care, because depending on the selection method the sample is likely to be at least somewhat incomplete. Halos at a given mass will show some scatter in their luminosity, and downward scatter (e.g., a low UV luminosity associated with a temporarily low SFR) will remove galaxies from the sample. Upward scatter can also add lower mass galaxies to the sample, but the effects are not symmetric because, if star formation histories are bursty, then galaxies will tend to spend more time with SFRs below their long-term average than with SFRs above their average. Van Dokkum et al. (2008) estimate that ∼50% of galaxies at z ∼ 2 will have suppressed SFRs, so that the appropriate value of n₀ for a UV-selected sample such as that of Papovich et al. is a factor of ∼2 above the nominal value. The samples of Marchesini et al. are stellar mass selected and those of Lundgren et al. are H-band selected, and should therefore be substantially more complete; for these samples the appropriate value of n₀ to choose is probably close to the nominal one for the survey. Thus, we plot curves for our model galaxies selected using n_{0, model} = 2 × 10⁻⁴, 4 × 10⁻⁴, and 8 × 10⁻⁴ Mpc⁻³, corresponding to perfect completeness, factor of two incompleteness, and factor of four incompleteness in the observations, respectively. As the figures show, for reasonable estimates of observational completeness, the models do a good job reproducing the observed trends in SFR and stellar mass. In particular, we see that in both the models and the data, SFRs and stellar masses rise roughly exponentially with −z at high-z before leveling off near z ≈ 2–3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In fact, the evolution of SFR and stellar mass in the model halo with M_{12, 0} = 10 shown in Figures 1 and 2, which has a halo mass of 10¹¹ M_☉ at z ∼ 6, is very similar to the observed evolution seen in Figures 10 and 11.

Another common observational method is to examine the properties of galaxies at fixed stellar mass across a range in redshift (e.g., Papovich et al. 2006; Stark et al. 2009; González et al. 2010). To compare to observations of this sort, we pick a target stellar mass in the range ∼10⁹–10¹¹ M_☉ for which high-redshift observations are available, and at each z in our model grid we select the halo with stellar mass closest to the target value. We then plot the corresponding sSFR in Figure 12. For comparison, we also plot a region illustrating the range of observationally determined sSFR values based on a data compilation from Weinmann et al. (2011). We see the same general result as in Figure 9: In our fiducial model the sSFR has a nearly constant value of ∼2 Gyr⁻¹ from z ∼ 8 to z ∼ 2. This is in good agreement with the observations. In contrast, the specific gas infall rate $\mbox{sIR}=\dot{M}_{g,\rm in}/M_*$ in these halos, and the sSFR, if f_{g, in} = 1, is more than an order of magnitude higher at z ∼ 8 than at z ∼ 2, in disagreement with the observations. The underlying reason for this trend is the same as for the individual galaxies: at high-z, stellar mass is roughly proportional to halo mass because most stars are accreted, and SFR is also roughly proportional to halo mass because t_SF has only a weak redshift dependence. This provides an explanation for the otherwise puzzling observed sSFR plateau, as discussed in Sections 1 and 5.

Zoom In Zoom Out Reset image size

Yet another important observable is the correlation between stellar mass and SFR for star-forming galaxies, the "SFR sequence," at different redshifts. This is related to the sSFR versus redshift discussed in Section 4.2, but measured for a range of stellar masses at fixed redshifts rather than for a fixed stellar mass at a range of redshifts.

To extract this quantity from our model grid, we pick a redshift slice to correspond to a chosen observational study, then plot stellar mass versus sSFR for all halos at that redshift. Figure 13 shows the result, compared to several observational surveys: a sample of Lyman break galaxies at z ≈ 3.7 from Lee et al. (2011), a sample of BzK galaxies at z ≈ 2 from Daddi et al. (2007), a sample of z ≈ 1 galaxies from the GOODS survey studied by Elbaz et al. (2007), and z ≲ 0.1 galaxies from the Sloan Digital Sky Survey (SDSS) studied by Brinchmann et al. (2004). For the Daddi et al. (2007) and Lee et al. (2011) samples, we use the observational constraint regions given in Lee et al. (2011). For the Brinchmann et al. (2004) and Elbaz et al. (2007) samples, we use the fitting formulae given in Elbaz et al. (2007); since they do not give explicit confidence regions or stellar mass ranges, we set our mass range based on eyeball estimates of the region occupied by the data. Similarly, we adopt a scatter of 0.4 dex for the Elbaz et al. (2007) sample and 0.5 dex for the Brinchmann et al. (2004) sample, based on an eyeball estimate of the scatter in the data. We correct all observations to a Chabrier (2005) initial mass function (IMF), which is a factor of 1.59 lower in SFR than a Kennicutt (1998) IMF.

Zoom In Zoom Out Reset image size

The plot shows that our fiducial model does a reasonably good job of reproducing the observed sSFR sequence, particularly at high redshift. At lower redshift, our models fall slightly below the observed range, consistent with the similar undershoot of sSFR we saw below z = 2 in Section 4.2. Nonetheless, the level of agreement is gratifying given the extreme simplicity of our model. In particular, we recover the basic results that the sSFR is roughly independent of stellar mass at all redshifts, and that the value of the sSFR varies fairly slowly with redshift above z ∼ 2. The model also makes clear predictions for the slope and normalization of the SFR sequence at higher redshift than has yet been observed. We do not show the corresponding predictions for our variant model with f_{g, in} = 1, but, not surprisingly, based on the discussion in Section 4.2, it provides a much worse fit to the higher redshift data (and a somewhat better fit to the low-redshift data). The failure at high redshift is because of the lack of stellar accretion, which produces systematically smaller stellar masses at the same SFR.

Another observational constraint to which we can compare our model is the mass–metallicity relation. This is not tremendously sensitive to the inclusion or exclusion of metallicity-dependent star formation, because most of the galaxies for which this relation has been observed are in the regime where there is little metallicity quenching. However, it is important as a consistency check. Since the model depends on metallicity affecting the SFR, it is critical that the metallicity evolution the model generates be at least roughly consistent with the observed metallicities of galaxies as a function of mass and redshift.

Such an observational comparison is unfortunately very difficult because of uncertainties in the absolute zero points of the strong-line metallicity indicators that are generally used to infer galaxy metallicities. While the various metallicity measurements in use can be calibrated against each other to a scatter of ∼0.03 dex, the zero point is uncertain by as much as 0.7 dex (Kewley & Ellison 2008). Since our models predict an absolute metallicity, i.e., a mass of metals divided by a mass of gas, this means that there will be a large systematic uncertainty in any observational comparison. We handle this problem by fixing the models to match observations at a particular stellar mass and redshift, thereby fixing the zero point. We can then compare the mass and redshift dependence of the models to the observations. It is important to note that we are not looking for more than a very rough consistency here, because the slopes of metallicity versus mass and metallicity versus redshift depend on which metallicity calibration one adopts. Different choices of calibration can make a significant difference to how well the data fit the model.

Figure 14 shows a comparison between the mass–metallicity relation in our models and a compilation of observations from Tremonti et al. (2004) at z = 0.07, Erb et al. (2006) at z = 2.27, and Maiolino et al. (2008) at z = 3.7. As with other comparisons, we generate the model predictions by finding the redshift slice in our model grid nearest to the target value and plotting mass versus metallicity at that slice. We correct all the observed masses to a Chabrier (2005) IMF, and all the observed metallicities to the KD02 calibration scheme using the conversions of Kewley & Ellison (2008). We fix the absolute calibration of the models by forcing agreement between the models and the Tremonti et al. (2004) data at z = 0.07 and log M_{*, 11} = −0.5. This corresponds to adopting a conversion

$$\begin{equation} 12 + \log ({\rm O}/{\rm H}) = 10.75 + \log (M_Z/M_g), \end{equation} \tag{ 30 }$$

where M_Z and M_g are the masses of metals and gas, respectively.

Zoom In Zoom Out Reset image size

As the figure shows, the agreement between the observed and model metallicity evolution is generally good at high-z. In the present-day universe, our model produces a somewhat steeper mass–metallicity relationship than the observed SDSS sample. While this could indicate a real defect in the models, it could equally well be a result of the particular metallicity calibration scheme we have used, since the slope of the SDSS sample (and the other samples) can vary by factors of ∼2 depending on the calibration scheme used. If our models really are too flat compared to reality, this may be related to our models' tendency to somewhat underpredict the SFR in small galaxies at low redshift, as we will see in Sections 4.2 and 4.3. However, since the agreement between our models and the observations is good at z ≳ 2, our results should be robust there.

Another feature of our model is its ability to predict the H₂ content of galaxies as a function of redshift and other galactic properties. Figure 15 shows the H₂ mass fraction versus stellar mass predicted by our models. For comparison, at z = 0 we show the results of COLD GASS (Saintonge et al. 2011a, 2011b), a volume-limited survey of CO line emission (used as a proxy to infer H₂) in several hundred galaxies at distances of 100–200 Mpc, selected to have stellar masses M_* > 10¹⁰ M_☉. As the plot shows, we find generally good agreement between our model and the observed amount of H₂ in galaxies at z = 0, albeit with quite a large spread in the data. We recover the trend that the H₂ fraction is (very weakly) declining with stellar mass. We also note that the data contain a significant number of non-detections. Some of these represent star-forming galaxies with relatively modest molecular content whose true value of $M_{\rm H_2}/M_*$ is likely only a bit below the reported upper limit. However, some also represent quenched galaxies that are almost entirely devoid of star formation and molecular gas. Our model does not correctly produce galaxies of this type due to our assumption that accretion cuts off sharply at a single, relatively high, halo mass. The observed coexistence of star-forming and passive galaxies at the same stellar (and presumably halo) mass indicates that this is clearly an oversimplification.

Zoom In Zoom Out Reset image size

At higher redshift, we predict that the H₂ fraction at fixed stellar mass should rise, reaching ∼10% by mass at z = 1 and roughly 50% by mass at z = 2. This is qualitatively consistent with the observed trend toward high molecular gas fractions at z ≈ 2 (e.g., Tacconi et al. 2010), but given the selection biases inherent in the current high-z data, it is hard to draw any general conclusions yet. However, these predictions will be testable with ALMA in the next few years (e.g., see reviews by Combes 2010, 2011).

Another useful quantity to plot from our models are the mass functions of stars, H i, and H₂ as a function of epoch. Such a comparison can provide some physical insight into the behavior of the model. However, we strongly caution that we should not expect good agreement between the model and observations here. Our models are very simple. Unlike in a full SAM, we have not made any effort to tune parameters to reproduce observations, and we omit much physics that is known to be necessary to reproduce observed mass functions. For example, we assume that there is one galaxy per halo at all halo masses, clearly not a reasonable assumption for high-mass halos that should host clusters. We make the mass loading factor of galactic winds independent of halo mass or other properties, which is probably not reasonable for low masses. Finally, we assume that all halos at a given mass and redshift are the same, ignoring scatter, which is not reasonable given the steepness of the halo-mass function.

It is important to note that these omissions are unlikely to dramatically affect the other observational comparisons we present. In particular, we are not particularly concerned with galaxies larger than L_*, and these do not contribute greatly to the cosmic SFR budget, so our incorrect treatment of them is not that important. Similarly, our neglect of scatter can have dramatic effects on the mass function, but it is unlikely to affect correlations between, for example, mass and metallicity. The one factor that could have a significant effect is our assumption of weak galactic winds at low masses, since these can have much the same effect as metallicity quenching in suppressing star formation in small galaxies. In reality both metallicity quenching and strong winds are likely to be important, but we have deliberately kept the winds weak and mass-independent in order to perform a cleaner experiment on the effects of metallicity alone.

To construct our mass functions, we compute the halo abundance n(M_h, z) from the Sheth & Tormen (1999) approximation as in Section 4.1, and for each halo we assign the H i, H₂, and stellar mass of the corresponding halo in our grid. Since the stellar mass is a monotonic function of the halo mass in our models, at any redshift z the number density of galaxies n(M_*, z) in the mass range from M_* to M_* + dM_* is simply given by

$$\begin{equation} n(M_*,z) = n(M_h,z) \frac{dM_*}{dM_h}, \end{equation} \tag{ 31 }$$

where on the right-hand side M_h is the halo mass for which the stellar mass is M_*, and the derivative is evaluated as this mass.

The atomic and molecular gas masses are not strictly monotonic in the halo mass in our model, so there may be multiple halo masses M_{h, 1}, M_{h, 2}, ... for which the corresponding atomic or molecular mass is $M_{{\rm H\,\mathsc{i}}}$ or $M_{\rm H_2}$. The equivalent expression for the number densities of galaxies with a given H i mass is therefore

$$\begin{equation} n(M_{{\rm H\,\mathsc{i}}}, z) = \sum _i n(M_{h,i},z) \left.\frac{dM_{{\rm H\,\mathsc{i}}}}{dM_h}\right|_{M_{\rm h,i}}, \end{equation} \tag{ 32 }$$

where the sum runs over all halo masses M_{h, i} for which the H i mass is $M_{{\rm H\,\mathsc{i}}}$. The expression for H₂ is analogous.

Figure 16 shows the model mass functions at three redshifts, compared to data where it is available. The plots reveal several interesting features of the model. First, examining the stellar mass function at z = 0, we see the sharp drop at low stellar masses corresponding to halos where star formation begins to be suppressed by metallicity effects. The mass at which this cutoff occurs moves to smaller values at higher redshifts. Well above this cutoff, our model provides very little suppression of star formation, which is not surprising given that we have not included strong galactic winds or similar mechanisms, and that we do not have any metallicity scatter in our models. Thus, we suppress star formation too much below the cutoff mass and too little above it. This shows up in reverse in the H i mass function, where we overpredict the number of low H i mass galaxies, and underpredict the number of large H i mass galaxies. Clearly a sharp cutoff in star formation such as the one produced by our model is at best a crude approximation to reality, and the real effects of metallicity should be spread out much more in halo mass.

Zoom In Zoom Out Reset image size

We also fail to produce galaxies with very large H i masses, indicating that we are probably converting H i to H₂ and thence into stars too efficiently in some cases. The origin of this failure is not entirely clear. One possibility is that it comes from our lack of scatter in spin parameter. Galaxies with large H i masses tend to have very large, extended H i disks that are non-star forming. Our model does not produce these.

By combining the SFR as a function of halo mass and redshift from our model grid with a calculation of the halo abundance as a function of mass and redshift, we are able to produce the cosmic star formation history for our toy model. We caution that, since the agreement between our simple model and observed sSFRs at redshifts <2 is only approximate, we should expect no better agreement here. Our goal is simply to understand the qualitative effect of metallicity quenching on cosmic star formation history with an emphasis on the high redshifts where its effects are greatest.

We compute the halo number density n(M_h, z) as in Section 4.1. The cosmic SFR density at redshift z is then given by

$$\begin{equation} \dot{\rho }_* = \int \dot{M}_{*,\rm form}(M_h,z) n(M_h,z) M_h \, d\ln M_h. \end{equation} \tag{ 33 }$$

The cosmic density of gas accretion into galaxies is given by the analogous expression with $\dot{M}_{*,\rm form}$ replaced with $\dot{M}_{g,\rm in}$. Note that this expression omits stars formed in the accretion shock when f_{g, in} ≠ 1, but that is a modest effect.

We plot the integrand of Equation (33), and the analogous accretion rate for the cosmic gas accretion density, in Figure 17. The figure illustrates several points. First, our fiducial model range of masses does a good job of sampling the halo masses that dominate the SFR, and thus we can make an accurate estimate of ρ_* from it for the fiducial model. Second, it is obvious upon inspection that, for the fiducial model, the SFR density of the universe will peak at z ∼ 2–3 and fall off sharply on either side of it. On the low-z side, this fall-off is driven by the decline in SFR (which is in turn driven by the fall off in gas accretion rate). On the high-z side, the decline in SFR is driven by the suppression of star formation due to the difficulty in forming a cold, molecular phase in small halos. In comparison, if we ignore the need to form such a phase and set $f_{\rm H_2} = 1$, or if we were to assume that halos could form stars as quickly as they could accrete gas, the cosmic star formation history would be much less peaked toward z ∼ 2–3. At high redshift, it would be dominated by small halos undergoing rapid accretion and star formation.

Zoom In Zoom Out Reset image size

Figure 18 shows the cosmic SFR density as a function of redshift integrated over all halos, both as actually observed and in our variant models, and the cosmic accretion density into halos above a given mass. These latter lines show what we would expect if halos simply turned gas into stars as quickly as they accreted it from the IGM; the model with a minimum halo mass M_{h, 12} = 0.1 is a rough approximation of the model presented by Bouché et al. (2010), although it is simplified in that it assumes a star formation time of zero as opposed to some number of galactic orbital periods as Bouché et al. assume.

Zoom In Zoom Out Reset image size

As the plot shows, all the models recover the observed SFR from z = 0–2 reasonably well. Given the long timescales from z = 2 to the present, the SFR in galaxies at this epoch becomes limited only by their gas supply, and it is the decline in gas supply from the z = 2 to the present that drives this trend. Note that our model SFR is slightly above the best-fit observations at z = 2, but the exact shape and value of the peak in our models depend on how we choose to truncate cold gas accretion in massive halos near z = 2.

In contrast, the only models that are able to reproduce the observed SFR at z = 2–8 are those that include metallicity quenching (fiducial, f_{g, in} = 1 and M_ret = 0.03) and the ad hoc model with a threshold for accretion at M_{h, 12} = 0.1 (as in Bouché et al. 2010). The run with $f_{\rm H_2} = 1$ provides a more gradual rise of the SFR density at high redshift, similar to the ad hoc model with a threshold for accretion at M_{h, 12} = 0.01, indicating that the suppression because t_acc is shorter than t_SF is effective below a characteristic halo mass of ∼10¹⁰ M_☉. This is indeed comparable to the Press–Schechter typical halo mass at the end of this period, z ∼ 2, where the SFR finally catches up with the IR.⁹ In our fiducial model and its variants that also include metallicity quenching, we recover the steeper rise in time because star formation at high redshift is also suppressed by the difficulty in forming a cold molecular phase in small galaxies that have not yet enriched themselves with metals. The ζ_lo = 0.8 model is in between the fiducial model and the $f_{\rm H_2} = 1$ model, but is still marginally consistent with the observational constraints. The Z_IGM = 2 × 10⁻⁴ model (not shown in the plot) is very similar. Thus, we see that the amount by which star formation is suppressed when we adopt a realistic treatment of ISM thermodynamics is somewhat but not tremendously sensitive to how metals are retained in halos as a function of mass and to the metal enrichment history of the IGM.

In the simplified model where the SFR is set equal to the accretion rate into halos with masses below M_{h, 12} = 0.1, the mass limit has much the same effect as including metallicity in the star formation law on the integrated SFR of the universe. Both suppress star formation at high-z. This demonstrates that quenching below a threshold mass captures some of the effects of metallicity-dependent star formation. The model presented here provides a physical mechanism that induces a strong but more gradual mass dependence, which boils down to a similar effect on the cosmic SFR history. Our physical model, while it does reduce the SFR in small halos, does not make them entirely devoid of stars, and it still allows galaxies like the Small Magellanic Cloud which lives in a halo of mass M_{h, 12} ∼ 0.01 (Bekki & Stanimirović 2009). In reality, the suppression of star formation is probably even less sharp than in our models, since there are almost certainly large stochastic variations in the metal enrichment history of small halos that our simple deterministic model does not capture. The difference between our fiducial model and the ζ_lo = 0.8 variant suggests what effect this stochastic variation is likely to have. We should therefore emphasize that the value of the effective threshold mass is predicted by our model only as an order of magnitude estimate.