The Astrophysical Journal, 776:34 (9pp), 2013 October 10  C 2013.

doi:10.1088/0004-637X/776/1/34

The American Astronomical Society. All rights reserved. Printed in the U.S.A.

DWARF GALAXY FORMATION WITH H2 -REGULATED STAR FORMATION. II. GAS-RICH DARK GALAXIES AT REDSHIFT 2.5 1

Michael Kuhlen1,2 , Piero Madau2 , and Mark R. Krumholz2 Theoretical Astrophysics Center, University of California, Berkeley, CA 94720, USA; [email protected] 2 Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA Received 2013 May 23; accepted 2013 August 11; published 2013 September 24

ABSTRACT We present a cosmological hydrodynamic simulation of the formation of dwarf galaxies at redshifts z  2.5 using a physically motivated model for H2 -regulated star formation. Our simulation, performed using the Enzo code and reaching a peak resolution of 109 proper parsecs at z = 2.5, extends the results of Kuhlen et al. to significantly lower redshifts. We show that a star formation prescription regulated by the local H2 abundance leads to the suppression of star formation in dwarf galaxy halos with Mh  1010 M and to a large population of gas-rich “dark galaxies” at z = 2.5 with low star formation efficiencies and gas depletion timescales >20 Gyr. The fraction of dark galaxies is 60% at Mh  1010 M and increases rapidly with decreasing halo mass. Dark galaxies form late and their gaseous disks never reach the surface densities, 5700 M pc−2 (Z/10−3 Z )−0.88 , that are required to build a substantial molecular fraction. Despite this large population of dark galaxies, we show that our H2 -regulated simulation is consistent with both the observed luminosity function of galaxies and the cosmological mass density of neutral gas at z  2.5. Moreover, our results provide a theoretical explanation for the recent detection in fluorescent Lyα emission of gaseous systems at high redshift with little or no associated star formation. We further propose that H2 -regulation may offer a fresh solution to a number of outstanding “dwarf galaxy problems” in ΛCDM. In particular, H2 -regulation leads galaxy formation to become effectively stochastic on mass scales of Mh ∼ 1010 M , and thus these massive dwarfs are not “too big to fail.” Key words: cosmology: theory – galaxies: dwarf – galaxies: formation – galaxies: halos – methods: numerical Online-only material: color figures

(Ben´ıtez-Llambay et al. 2013), while able to quench star formation in satellite and isolated dwarfs at late times, must operate in concert with other suppression mechanisms that prevent the early conversion of gas into stars in such systems. Mass loss driven by supernova (SN) feedback can reduce the baryonic content of bright dwarfs (Dekel & Silk 1986; Mori et al. 2002; Governato et al. 2007) and flatten their central dark matter cusps (Read & Gilmore 2005; Mashchenko et al. 2008; Governato et al. 2010; Pontzen & Governato 2012; Teyssier et al. 2013), in all but the most star deficient dwarf spheroidals (Garrison-Kimmel et al. 2013; Zolotov et al. 2012; Boylan-Kolchin et al. 2012; Penarrubia et al. 2012). In Kuhlen et al. (2012, Paper I), we followed a different avenue to quench star formation in dwarfs, one based on a new understanding of the chemistry and thermodynamics of the interstellar gas that is actually forming stars. Spatially resolved observations of local galaxies have revealed that star formation correlates much more tightly with the density of molecular gas than with the total gas density (Leroy et al. 2008; Bigiel et al. 2008), even in regions where molecular gas constitutes only a trace component of the interstellar medium (ISM; Schruba et al. 2011). Even though the primary cooling agents are lines of CO or C ii (depending on the chemical state of the carbon), molecular hydrogen (H2 ) is expected to be good tracer of star formation (Krumholz et al. 2011; Glover & Clark 2012a, 2012b) because star formation occurs only where the gas is able to reach very low temperatures, which is possible only in regions where the extinction is high enough to block out the background interstellar radiation field that is responsible both for heating the gas and for dissociating H2 molecules. In Paper I, we showed that H2 -regulated star formation leads to the suppression of star formation in dwarf galaxy

1. INTRODUCTION For almost two decades, observations of dwarf galaxies have challenged the ΛCDM paradigm of cosmological structure formation and our understanding of the mapping from dark matter halos to their baryonic components. Dwarfs—either in the field or in the halos of larger systems like our own Milky Way—are much less abundant than the Mh  1010 M dark matter halos that should, in principle, be able to host them. Darkmatter-only simulations predict steep (“cuspy”) inner density profiles, while the observed rotation curves of dwarf galaxies instead suggest that they have near-constant density cores. Similarly, the most massive subhalos found in dark-matteronly simulations of Milky-Way-sized systems appear to be too dense to be consistent with existing constraints on the classical Galactic dwarf satellites. While the discrepancies above may be truly reflective of a fundamental failure of the otherwise remarkably successful ΛCDM model, numerous astrophysical solutions have been proposed to explain this “dwarf galaxy problem” and are being actively investigated. To date, they all appear to provide only a partial, often environment-dependent solution to the puzzle. Photo-heating from the cosmic ultraviolet background may suppress gas infall onto dwarf galaxy halos in the low-redshift universe (Efstathiou 1992; Thoul & Weinberg 1996), but the importance of photoionization feedback is greatly reduced at high redshift because dwarf galaxy-sized objects can either self-shield against the ionizing background (Dijkstra et al. 2004) or accrete substantial amount of gas prior to reionization (Bullock et al. 2000; Madau et al. 2008). Gas removal by ram pressure in massive galaxy halos (Mayer et al. 2006) or through hydrodynamical interactions with the cosmic web 1

The Astrophysical Journal, 776:34 (9pp), 2013 October 10

Kuhlen, Madau, & Krumholz

halos at z  4. We discussed how such a novel quenching mechanism—one that modifies the efficiency of star formation of cold gas directly, rather than indirectly reducing the cold gas content with SN feedback—may contribute to alleviate some of the issues faced by theoretical galaxy formation models. Similar models have been investigated numerically by Gnedin et al. (2009), Gnedin & Kravtsov (2010, 2011), and Christensen et al. (2012), and semi-analytically by Fu et al. (2010) and Krumholz & Dekel (2012). Our work in Paper I and here differs from the previous numerical simulations in that we investigate a cosmologically representative volume rather than single galaxies; we are therefore able to discuss the statistical properties of the galaxy population produced by H2 -regulated star formation. At the same time, we are still performing full three-dimensional simulations, rather than relying on semianalytic model prescriptions. Here, we extend our state-of-theart simulations to lower redshifts and confirm the existence at these epochs of a large population of very inefficient star formers below Mh  1010 M . We show that such “dark galaxies” form late and that their gaseous disks never reach the surface densities required to build a substantial molecular fraction.

We enforce a minimum stellar mass of mmin = 104 M , since even Δt0 can occasionally become very small. Below this mass, we implement a stochastic star formation criterion as follows: if m∗ < mmin , we form a particle of mass equal to mmin if a randomly generated number is smaller than (m∗ /mmin ). To obtain the molecular hydrogen mass fraction fH2 in a given grid cell, we use the two-phase equilibrium model of Krumholz et al. (2008, 2009) and McKee & Krumholz (2010), hereafter referred to as the KMT model. This model is based on a radiative transfer calculation of an idealized spherical giant atomic–molecular complex, subject to a uniform and isotropic Lyman–Werner radiation field. The H2 abundance is calculated assuming formation–dissociation balance and a twophase equilibrium between a cold neutral medium (CNM) and a warm neutral medium (WNM; Wolfire et al. 2003): fH2  1 −

ln(1 + 0.6χ + 0.01χ 2 ) , 0.6 τc   σd,−21 1 + 3.1 (Z/ZSN )0.365 . χ = 2.3 R−16.5 φCNM s=

2. SIMULATIONS AND MOLECULAR CHEMISTRY We summarize here the main features of our cosmological adaptive mesh refinement hydrodynamics simulation (see Kuhlen et al. 2012 for a more extensive discussion) with the Enzo v2.2 code3 . The computational domain covers a (12.5 Mpc)3 box with a root grid of 2563 grid cells. The dark matter density field is resolved with 2563 particles of mass 3.1 × 106 M . Adaptive mesh refinement is allowed to occur throughout the entire domain for a maximum of seven levels of refinement, resulting in a maximum spatial resolution of Δx7 = 109 × 3.5/(1 + z) proper parsecs. Mesh refinement is triggered by a grid cell reaching either a dark matter mass equal to four times the mean root grid cell dark matter mass or a baryonic mass equal to 8 × 2−0.4l times the mean root grid cell baryonic mass, where l is the grid level. The simulation is initialized at z = 99 with cosmological parameters consistent with the WMAP7 year results (Komatsu et al. 2011). It includes radiative cooling from both primordial and metal enriched gas, as well as photo-heating from an optically thin, uniform metagalactic UV background. Artificial pressure support is applied to cells that have reached the maximum refinement level in order to stabilize these cells against artificial fragmentation. We implement an H2 -regulated star formation prescription, as follows. The local star formation rate (hereafter SFR) in a grid cell is proportional to the molecular hydrogen density by the free-fall time, tff , determined from the total gas density, ρgas ρ˙SF =  fH2 , (1) tff

(3)

(4) (5)

τc  0.067 (Z/ Z ) (ΣHI / M pc−2 ) is the dust optical depth of the cloud, σd,−21 is the dust cross section per H nucleus to 1000 Å radiation, normalized to a value of 10−21 cm−2 , R−16.5 is the rate coefficient for H2 formation on dust grains, normalized to the Milky Way value of 10−16.5 cm3 s−1 (Wolfire et al. 2008), ZSN is the gas phase metallicity in the solar neighborhood, and φCNM is the ratio of the typical CNM density to the minimum density at which a two-phase CNM–WNM equilibrium can be established. We set ZSN = Z (Rodr´ıguez & Delgado-Inglada 2011), Z = 0.0204, φCNM = 3, and (σd,−21 /R−16.5 ) = 1. Note that in this model the H2 fraction is independent of the intensity of the ambient UV radiation field (which is responsible for both heating and H2 dissociation), since under the assumption of a two-phase CNM–WNM equilibrium the CNM density scales approximately linearly with the UV intensity and the equilibrium H2 abundance depends on the ratio of the two. Krumholz & Gnedin (2011) demonstrated that this simple analytical model very accurately reproduces the H2 abundance determined in full non-equilibrium radiative transfer calculations when Z/ Z  10−2 . At lower metallicities, the equilibrium model tends to overpredict the H2 abundance, but at such low metallicities star formation should nevertheless scale with the equilibrium H2 abundance, not with the nonequilibrium value (Krumholz 2012). While in the KMT model fH2 drops to zero at a metallicitydependent critical column density of

where fH2 = ρH2 /ρgas is the H2 fraction, computed as described below, and the star formation efficiency is  = 0.01, the value favored by observations (Krumholz & Tan 2007; Krumholz et al. 2012). No density threshold for star formation is applied. Stars form only once per root grid time step and only in cells at the highest refinement level (here, l = 7), but the star particle mass is proportional to the root grid time step Δt0 (see Kravtsov 2003), i.e., Δt0 . (2) m∗ =  fH2 ρgas (Δx7 )3 tff 3

s 3 , 4 1 + 0.25s

log(1 + 0.6χ + 0.01χ 2 ) M pc−2 0.0804 Z/Z   Z/ Z −0.88 ≈ 5700 M pc−2 , 10−3

Σcrit (Z) =

(6)

observations indicate that such a hard cutoff is too extreme. In particular, Bigiel et al. (2010) have demonstrated that the scaling between total gas surface density (Σgas ) and SFR surface density (ΣSFR ) continues below the turnover commonly attributed to the transition of the atomic to molecular gas phase ( few M pc−2 at solar metallicity), albeit at a ∼50–100 times lower amplitude.

http://enzo-project.org/

2

The Astrophysical Journal, 776:34 (9pp), 2013 October 10

Kuhlen, Madau, & Krumholz

Figure 1. Stellar mass fraction f = M /Mh (left panel), SFR (middle panel), and gas depletion time tdep = MHI+H2 /SFR (right panel) vs. total halo mass, all at z = 2.5. In the left panel, the dotted horizontal line indicates an f equal to the cosmic baryon fraction Ωb /ΩM , i.e., a 100% gas-to-star conversion efficiency. The dashed lines in the left and middle panels indicate the lowest f and SFR that our simulation can resolve, given our minimum stellar particle mass of (1 − fej ) msp = 7.5 × 103 M and star formation averaging timescale of 10 Myr. In the right panel, the horizontal dotted line marks our definition of a “dark galaxy,” i.e., one with tdep > 20 Gyr. Completely dark halos (M = 0, SFR= 0, or tdep = ∞) are indicated by arrows at the location corresponding to the lowest resolvable M or SFR (for clarity, we only plotted 1/10 of all such halos). (A color version of this figure is available in the online journal.)

Assuming the molecular gas–SFR scaling relation measured by Schruba et al. (2011) continues to correspondingly low values of ΣSFR (10−4 M yr−1 kpc−2 ), the gas at these extremely low total gas columns has an effective molecular gas fraction of fH2 ≈ 0.01. In order to capture this empirical low surface density behavior, we impose a floor of fH2 = 0.01 in cells with T < 104 K even if they have Σgas < Σcrit (Z). The resolution of our simulation is not sufficient to resolve the formation sites of the first generation of stars, the so-called Population III stars. In order to capture the metal enrichment resulting from the SN explosions of this primordial stellar population, we instantaneously introduce a metallicity floor of [Zfloor ] ≡ log10 (Zfloor / Z ) = −3.0 at z = 10, as motivated by recent high-resolution numerical simulations of the transition from Population III to Population II star formation (Wise et al. 2012). This method ensures the presence of a minimum amount of metals, which seed subsequent star formation and further metal enrichment. The simulation was run on the NASA supercomputer Pleiades and required ∼200,000 core hours to run to z = 2.5.

by the difficulty that low metallicity gas has in transitioning to a cold, molecular phase that is capable of forming stars. In our simulation, the stellar mass fraction drops rapidly below 1010 M and becomes highly stochastic, with f spanning about two orders of magnitude. The mean star formation efficiency drops from 0.046 to 1.6×10−3 when averaged over all halos with mass 1010 –1011 M and 109 –1010 M , respectively. A small number of halos (8 out of 1010) have not been able to form any stars and we have marked these halos as upper limits at the lowest f that our simulation can resolve given our minimum stellar particle mass of msp = 104 M , of which fej = 0.25 is ejected back into the ISM as part of the SN feedback prescription. The middle panel of Figure 1 shows the total SFR versus halo mass. Again, the effect of our metallicity-dependent, H2 -regulated star formation prescription is clearly visible as a drop in SFR around Mh = 1010 M . At higher masses, the SFR falls on a narrow “SFR main sequence” corresponding to a roughly constant halo mass-specific SFR of ≈4 × 10−11 yr−1 . At lower halo masses, the SFR drops by almost two orders of magnitude, but remains non-zero owing to our 1% H2 floor. Slightly less than half of our halos (396 out of 1010) have a zero SFR, and one-tenth of these are plotted as upper limit arrows at the lowest SFR our simulation can resolve, 7.5 × 10−4 M yr−1 (i.e., forming 1 star particle over the star formation averaging timescale of 10 Myr). Finally, the right panel of Figure 1 shows the neutral hydrogen gas depletion time, defined as tdep = (MH i + MH2 )/SFR. On the “SFR main sequence” (Mh  1010 M ), we find gas depletion times ranging from 0.4 to a few Gyr, consistent with the observational constraints from Leroy et al. (2008), Bigiel et al. (2008), and Genzel et al. (2010). The star formation-suppressed galaxies make up a second population with depletion times in excess of 20 Gyr. Owing to their low stellar content and suppressed SFR, such galaxies would be extremely difficult to detect with continuum IR, optical, or UV observations, and hence we refer to them as “dark galaxies.” In this work, we use tdep > 20 Gyr as a working definition of a dark galaxy. With this definition, 789 out of our 1010 halos are dark galaxies. To first order, dark galaxies lie in halos that rely on the 1% H2 floor to form any stars, while halos that become metal enriched enough to enter the KMT regime (Σgas > Σcrit (Z)) quickly form a significant stellar population and are luminous.

3. DARK GALAXIES Our previous work in Paper I was limited to z  4, but we have since extended the simulations to z = 2.5. The left panel of Figure 1 shows a scatter plot of the stellar mass fraction (f = M /Mh ) against the total halo mass Mh . Only halos in which the refinement reached the maximum refinement level of max = 7 are analyzed and there are 1010 such halos in our simulation at z = 2.5. At high halo masses (>1011 M ), about half of the available gas has been converted to stars. Such a large star formation efficiency appears to be in conflict with observational constraints derived from abundance matching (e.g., Moster et al. 2010; Behroozi et al. 2012), which imply a peak efficiency of a few percent for 1012 M halos. This discrepancy is a well-known problem with hydrodynamic galaxy formation simulations that include only a weak form of thermal SN feedback (see, e.g., Katz 1992) and much contemporary work is concerned with solving this problem (Stinson et al. 2013; Hopkins et al. 2012; Agertz et al. 2012; Guedes et al. 2011), primarily by improving the treatment of feedback from young massive stars, SNe, and black holes. Here, we defer addressing this problem and instead wish to draw attention to the fact that a suppression of the gas-to-star conversion efficiency in low-mass halos can be effected simply 3

The Astrophysical Journal, 776:34 (9pp), 2013 October 10

Kuhlen, Madau, & Krumholz

Figure 3. Maximum column density, Σmax = ρmax × Δx7 = ρmax × 109 pc (3.5/(1 + z)), reached in the simulated halos as a function of redshift. The size and color of the plotted symbols correspond to the mass of the halo. Luminous galaxies are shown in color; dark galaxies are shown in white. The horizontal dotted lines show the critical column density required for H2 formation, at metallicities (from top to bottom) Z/ Z = 10−3 , 10−2 , 10−1 , and 1. The curved solid lines show estimates of the maximum central surface density (averaged over Δx7 ) calculated using the analytical model by Mo et al. (1998) for halos of mass 1011 , 1010 , and 109 M (see the text for details). The solid lines are for a spin parameter λ = 0.035 and the gray shaded band delimits the region between λ = 0.02 and 0.05 for the 1010 M case. The dashed line gives the analytical estimate for a fixed central averaging scale of 50 pc. This plot shows that (1) at a given redshift, more massive halos are more easily able to convert gas to the molecular phase and hence form stars and (2) at a given halo mass, it is easier for gas to become molecular at higher redshift, which results in an increase, toward lower redshifts, of the critical halo mass below which star formation is not possible in very low metallicity systems. (A color version of this figure is available in the online journal.) Figure 2. Top panel: the fraction of dark galaxies (tdep > 20 Gyr) vs. total halo mass at different redshifts. The circles mark the halo mass at which fdark = 0.5. Bottom panel: the average mass fraction of total baryons, neutral, ionized, and molecular hydrogen, and stars vs. total halo mass, at z = 2.5. These lines were calculated using a sliding window of width 0.5 dex. The cosmic baryon fraction (Ωb /ΩM ) is indicated by the dotted horizontal line. (A color version of this figure is available in the online journal.)

systems. At lower halo masses, the effects of the metagalactic UV background are apparent based on a lowering of the neutral relative to ionized hydrogen abundance and a reduction in the total baryon fraction. In our model, a halo’s ability to become luminous depends on the maximum H i surface density its gas can reach. To demonstrate this fact, we show in Figure 3 a plot of the maximum H i density in each halo, for a range of redshifts. To obtain a simple estimate of the column density, we multiply the gas density by the width of the most refined grid cell, Σmax = ρmax × Δx7 = ρmax × 109 pc (3.5/(1 + z)). This quantity is plotted as circles whose sizes indicate the halo’s mass; the circles are filled blue (white) if the halo is luminous (a dark galaxy). Σmax should be compared with the KMT critical column density that allows the atomic-to-molecular conversion to occur, which is metallicity dependent and plotted with dotted horizontal lines for metallicities Z/ Z = 10−3 , 10−2 , 10−1 , and 1: Σcrit = 5700, 690, 92.8, and 13.6 M pc−2 . The lowest metallicity (10−3 Z ) corresponds to the metallicity floor that we impose at z = 10, implying that halos must at some point exceed this threshold in order for a substantial number of stars to form. Subsequent star formation is much easier, since the metal enrichment from these stars rapidly lowers the required density. Halos that never exceed the 10−3 Z column density threshold are dark galaxies. The figure illustrates two points: (1) that more massive halos are more easily able to convert atomic into molecular gas and hence support star formation and (2) that the

The top panel of Figure 2 shows fdark , the fraction of dark galaxies (tdep > 20 Gyr), at a range of redshifts. At z = 2.5, fdark rises from zero at 5 × 1010 M to near unity at Mh  5 × 109 M . The transition from mostly luminous to mostly dark halos shifts to lower halos masses at higher redshifts. A remarkable prediction of our work then is that the universe should be filled with a large number of dark galaxies of mass 1010 M that were unable to convert their gas into stars simply because their metal content never rose to the level required for a conversion from the atomic to the molecular phase. Unlike in the more conventional feedback-based star formation suppression models, in our picture many of these dark galaxies would retain nearly the cosmic baryon budget in the form of neutral hydrogen, making them potentially observable in 21 cm or Lyα fluorescence. Indeed, Cantalupo et al. (2012) may have already discovered examples of such dark galaxies in Lyα fluorescence illuminated by a z ≈ 2.4 quasar. The bottom panel of Figure 2 shows the mean baryonic content of our halos as a function of total halo mass. As discussed above, a drop in f occurs around 1010 M , and it closely follows the transition between predominantly atomic gas in lower mass systems and mostly molecular gas in higher mass 4

The Astrophysical Journal, 776:34 (9pp), 2013 October 10

Kuhlen, Madau, & Krumholz

threshold halo mass for star formation increases toward lower redshifts (as also demonstrated in the top panel of Figure 2).

We compare this metallicity with Zcrit (Σmax ), the KMT critical metallicity for atomic-to-molecular gas conversion given the maximum density reached in each of our dark galaxies. If Zwind > Zcrit , then we judge this halo to have been externally enriched to the point where stars would have formed, and hence it would no longer be a dark galaxy. Using this procedure, we find that for fwind = 3 only 7 out 789 dark galaxies would be externally enriched to above their critical metallicity. For fwind = 1 (5), this number becomes 44 (1). The fraction of externally enriched halos becomes smaller for larger fwind , because even though the enriched region extend farther and intercepts more dark halos, the wind SOI’s metallicity becomes correspondingly smaller. It thus appears that the majority of our dark halos are isolated enough to remain unenriched, and hence dark, even in the presence of strong galactic outflows. Note that external enrichment will result in a strong spatial clustering of the small fraction of low mass halos that are able to become luminous in this way.

4. DISCUSSION 4.1. Analytical Model To gain an understanding of the trends seen in Figure 3, we have calculated the expected central surface density using the analytical model of Mo et al. (1998). In this model, the baryonic disk is completely specified by its host halo’s properties and is modeled as an exponential disk embedded in the potential of a Navarro–Frenk–White dark matter halo.4 The disk mass is set to md = 0.05 times the total halo mass and the disk scale radius is determined in an iterative manner from the assumption that the baryons contain a fraction jd = 0.05 of the total angular momentum J of the dark matter halo, which itself is related to the halo spin parameter λ = J |E|1/2 G−1 M −5/2 . Numerical simulations have shown that λ is distributed approximately lognormally, with a mean of λ ≈ 0.035 and a dispersion of ∼0.5 (Knebe & Power 2008). Given this analytic disk model, we can estimate the central gas surface density on a length scale comparable to the most refined grid cells in our simulation. R The solid lines in Figure 3 show ΣR = [ 0 2r Σ(r) dr]/R 2 for R = Δx7 in halos of mass 1011 , 1010 , and 109 M with a concentration c = 5. The shaded band denotes the region spanned by λ = 0.02 and 0.05. The analytical model agrees with our simulation results quite well and explains why a halo of a fixed mass and metallicity has a harder time converting its gas to the molecular phase at lower redshifts—the disk surface density is set by the mean density of the halo, which decreases toward lower redshift. It follows that a halo with an initial metallicity equal to our metallicity floor of Z/ Z = 10−3 can only become luminous if either it grows to a large enough mass that its disk surface density exceeds Σcrit = 5700 M pc−2 or, alternatively, if it becomes externally enriched from a nearby star forming system. The star formation afforded by the low molecular gas content resulting from our 1% H2 floor by itself is generally too low to allow the halo to become luminous.

4.3. z = 2.5 UV Luminosity Function Next, we confront our simulated galaxy population with the observationally determined rest-frame UV luminosity function (LF) of Lyman break galaxies at z = 1.9–2.7 from Reddy & Steidel (2009) and the very recent photometric measurement at z ≈ 2 by Alavi et al. (2013). For this purpose, we calculate UV luminosities from the simulated SFR assuming a Z = 0.2 Z Chabrier initial mass function (IMF) from 0.08 to 120 M and a constant SFR for 100 Myr (Madau et al. 1998): LUV = 1.25×1028 (SFR/M yr−1 ) erg s−1 Hz−1 . A Salpeter IMF would result in ∼1.7 times lower luminosities. We obtain SFRs from our simulated galaxies by summing the mass of all young star particles with ages less than τ = 10 Myr and dividing by this star formation timescale, SFR = age < τ m /τ . Star forming galaxies at z ≈ 2.5 are known to be significantly reddened and extincted by internal dust obscuration. Reddy et al. (2012) use a combination of Herschel, Very Large Array, and Spitzer data of 146 UV-selected galaxies with spectroscopic redshifts in the GOODS-North field to derive a median dust correction factor of 5.2 ± 0.6 between the bolometric luminosity and the unobscured (i.e., observed) UV luminosity. Consistent with this measurement, we dim our SFR-derived UV luminosities by a factor of five (ΔMdust = +1.75) and then calculate a rest-frame LF by binning our galaxies in bins of width 0.5 MUV (Figure 4). The agreement between the simulated and observed LFs is remarkably good down to the spectroscopic observational limit of −18 mag from Reddy & Steidel (2009). Owing to our fairly small box size, we do not have any simulated galaxies with dust-corrected UV luminosities brighter than −21 mag. The suppression of star formation due to the inability of low metallicity gas to become molecular in low-mass halos results in a flattening and eventual cutoff at the faint end of the UV LF, which occurs 1–2 mag below the current observational limit from spectroscopic studies. The agreement with the observational LF at z ≈ 2.5 may be puzzling, given the unrealistically high f in our most massive halos. The only way for these two facts to be consistent is for the SFR to have been too high in the past. Indeed, Figure 18 of Paper I indicates that the high redshift (z  4) LFs in our simulations exceed the observational constraints. Very recent photometric work by Alavi et al. (2013) has extended the z ∼ 2 LF to systems that are several magnitudes fainter, by taking advantage of the large magnifications provided by the strongly lensing galaxy cluster A1689. Their work

4.2. External Metal Enrichment? The possibility of external enrichment deserves further consideration. As discussed above, our simulation has fairly weak stellar feedback and, as a result, does not exhibit strong metal enriched galactic outflows. In order to estimate how many of our dark galaxies would have been externally enriched by more realistic galactic outflows from neighboring galaxies, we perform the following simple calculation: for every luminous halo, we calculate its escape velocity (vesc = (2GMh /Rvir )1/2 ), as well as the mean age of all its star particles, t . Assuming a galactic wind outflow with a velocity vwind = fwind vesc , we determine the radius of its wind sphere of influence (SOI), Rwind = vwind t . We estimate the mean metallicity of the wind by uniformly spreading the total metal mass contained in the luminous halo over the wind SOI. Next, we determine for each dark galaxy how many SOIs it intersects and we sum the metallicities of each of these to produce an estimate of the metallicity Zwind that this dark galaxy would have been externally enriched to. 4

We neglect the treatment of adiabatic contraction, since observationally there does not seem to be much evidence for it in late type galaxies (e.g., Dutton et al. 2011).

5

The Astrophysical Journal, 776:34 (9pp), 2013 October 10

Kuhlen, Madau, & Krumholz

UV background, which causes a declining neutral gas fraction toward lower mass halos (see Figure 2), thus compensating for the higher abundance of such halos. Dark galaxies make up the vast majority of all MH i < 109 M systems. Comparing with the observational measurement of the H i mass function in the local universe from the ALFALFA survey (Haynes et al. 2011), we see an excess of about one order of magnitude at MH i between 107 and a few ×109 M in our simulation. Furthermore, since our simulation overproduces the stellar content of the brightest galaxies (see Section 3), those same galaxies have values of M /MH i that are too high compared with the relation observed in present-day galaxies (Huang et al. 2012). Note that we are comparing here z ≈ 0 observational measurements with simulation results from an epoch ∼11 Gyr earlier, so this comparison is not entirely fair. Nevertheless, the excess of low MH i systems is somewhat troubling, since most physical processes that are commonly proposed to lower the H i content in low-mass halos (e.g., radiative feedback from young stars or SN driven galactic outflows; Altay et al. 2011; Dav´e et al. 2013) rely on a substantial SFR, which is not present in our dark galaxies. It is possible that over the following 11 Gyr the low residual star formation in our dark galaxies (stemming from the H2 floor) or external enrichment from nearby luminous halos will raise the metal content of these systems to the point where the molecular transition and, hence, more efficient star formation, becomes possible. Additionally, interactions with the cosmic web could ram pressure strip gas from these dark galaxies (Ben´ıtez-Llambay et al. 2013) and reduce their H i content. In the middle panel of Figure 5, we show the differential contribution to ΩH i per decade of MH i (i.e., MH i φ(MH i )/ρcrit ) from all halos and from dark halos only. In our simulation, the dominant contribution to ΩH i comes from systems with MH i between 3 × 108 and 109 M and luminous halos contribute about 1.5 times as much to the total H i budget than dark galaxies. The shape of our distribution closely matches that of the observational distribution from the ALFALFA survey (see Figure 9 of Martin et al. 2010), which, however, is peaked at about half a dex higher MH i . Again, we attribute this shift to the comparatively small box size of our simulation. In the right panel, we show a comparison of our simulation results with observational measurements of ΩH i from quasar absorption line studies from Rao et al. (2006) at z < 2 and Prochaska & Wolfe (2009) and Noterdaeme et al. (2012) at z > 2. The total H i content of our simulations is in reasonable agreement with the data. At z = 2.5, the content is about 0.2–0.4 dex (a factor of 1.5–2.5) too high; this fact is likely a consequence of the weak SN feedback in our simulation (see also Altay et al. 2011; Erkal et al. 2012; Dav´e et al. 2013). We re-emphasize that the dark galaxies in our simulation only contribute sub-dominantly to the total H i content of the simulation (see the open circles in Figure 5) and the absence of radiative feedback or SN-driven outflows in these systems is not a problem for the H i budget.

Figure 4. Simulated UV luminosity function at z = 2.5 compared with observational determinations from Reddy & Steidel (2009) and Alavi et al. (2013). We calculate LUV from the SFR according to LUV = 1.25 × 1028 (SFR/ M yr−1 ) erg s−1 Hz−1 (MUV = 51.63–2.5 log10 (LUV / erg s−1 Hz−1 )), corresponding to a Z = 0.2 Z Chabrier IMF from 0.08 to 120 M . The simulated galaxies’ luminosities have been dimmed by a factor of five (ΔM = +1.75) to account for the expected attenuation from dust, consistent with the observational constraints from Reddy et al. (2012). Note that the faintest two data points from Alavi et al. (2013) consist of only four systems and are likely subject to considerable systematic uncertainties arising from the lens mass modeling. (A color version of this figure is available in the online journal.)

indicates a steeply rising faint end of the LF down to MUV ≈ −13, in apparent conflict with the predictions of our H2 -regulated star formation simulation. We caution, however, that the faintest two data points (MUV > −15) of the Alavi et al. (2013) LF (the only ones disagreeing with our predictions) are comprised of only four systems and are likely subject to considerable systematic uncertainties, due to the difficulty of determining accurate lensing magnifications. These uncertainties affect both the inferred intrinsic magnitudes and the effective volume probed by the survey. Furthermore, the disagreement with the faintest two data points would be greatly reduced if we applied a luminosity-dependent dust attenuation law (see, e.g., Reddy & Steidel 2009) and dimmed MUV > −15 galaxies by less than a factor of five. For these reasons, we do not consider the current observational LF constraints to rule out the existence of a large number of dark galaxies at z ≈ 2.5. Our simulations, however, predict that the faint end of the UV LF should not continue to rise much beyond MUV ≈ −16. The Alavi et al. (2013) study highlights the potential of future deeper multicolor imaging surveys to constrain the H2 -regulated star formation scenario. 4.4. H i Abundance and Mass Function In Figure 5, we consider the abundance and mass distribution of neutral hydrogen in our simulation. The left panel shows the H i mass function at z = 2.5. We find halos with H i masses as low as 50 M and as high as 4.6 × 109 M . From MH i = 107 to 109 M , the H i mass function is remarkably flat (in contrast with the dark matter mass function, which rises steeply toward low masses) and exhibits a sharp cutoff at higher H i masses. This cutoff is a result of the limited box size of our simulation, which prohibits us from resolving halos more massive than ≈1012 M at z = 2.5. The flat shape of the mass function at lower masses is due to the ionizing and heating effects of the metagalactic

4.5. H i Depletion Time versus MH i Finally, in Figure 6, we show a comparison of the relation between H i depletion time and H i mass in our simulation at z = 2.5 with observational data both from high redshifts (z ≈ 2.5–3, top panel) and in the local universe (bottom panel). In the top panel, we include the “dark galaxy” systems reported by Cantalupo et al. (2012), which were detected as Lyα fluorescence in a z ≈ 2.4 quasar field. Assuming 6

The Astrophysical Journal, 776:34 (9pp), 2013 October 10

Kuhlen, Madau, & Krumholz

Figure 5. Left: the simulated H i mass function at z = 2.5 compared with the z ≈ 0 determination from the ALFALFA survey (Haynes et al. 2011). The solid line is for all simulated halos and the dashed is for all dark galaxies. Middle: the differential contribution to ΩH i per decade of MH i from all halos (solid) and dark galaxies (dashed). Right: ΩH i vs. redshift. Comparison of observational constraints from Rao et al. (2006; orange triangles) at z < 2 and Prochaska & Wolfe (2009; cyan circles) and Noterdaeme et al. (2012; red symbols; the filled triangles are systematics-corrected) at z > 2 with ΩH i determined in our simulation, summing over all halos (solid circles) and dark halos only (open circles). (A color version of this figure is available in the online journal.)

a negligible contribution to the Lyα luminosity from internal star formation or cooling, these authors derive a total hydrogen mass of ∼109 M (assumed to be mostly ionized). Without a nearby quasar to ionize it, most of this hydrogen would be neutral and so we equate it with MH i . The range of MH i derived from the Lyα luminosities is indicated by the horizontal error bar. Since no continuum emission was measured in these systems, Cantalupo et al. (2012) report only an upper limit of SFR < 0.01 M yr−1 derived from a stacking analysis and a corresponding lower limit on tdep,H i . This lower limit is consistent with the long depletion times of our simulated dark galaxies, but systems with MH i  109 M in our simulation tend to lie in >1010 M halos with more active star formation and shorter depletion times. Note, however, that the H i masses derived from the unresolved Lyα sources are quite uncertain since they depend on modeling the unknown spatial distribution of the fluorescing gas. We also plot in this panel a band corresponding to H i depletion times determined in z ≈ 3 damped Lyα (DLA) systems without central bulges of star formation from Rafelski et al. (2011). Since no length scales are known for these DLAs, we cannot obtain MH i estimates and instead show a horizontal band. The vertical extent of the band corresponds to the range of smoothing parameters employed by Rafelski et al. (2011), which result in different lower limits on tdep,H i . The band lies above the cloud of simulation points corresponding to actively star forming halos with short tdep,H i and is consistent with the much longer tdep,H i of our dark galaxies. In the lower panel of this plot, we compare with measurements from the local universe. Again, we caution that this comparison is not really fair, since a lot of evolution in both SFR and MH i may occur in the ∼11 Gyr between z = 2.5 and today. Nevertheless, it is interesting to see how our z = 2.5 dark galaxies compare to observational constraints in z = 0 dwarf galaxies. We plot data points from low surface brightness galaxies (Wyder et al. 2009), ALFALFA dwarf galaxies (Huang et al. 2012), the THINGS survey (Walter et al. 2008), and the newly discovered Leo P dwarf galaxy (Rhode et al. 2013). We see that a small fraction of the local dwarfs overlaps with the cloud of high tdep,H i simulation points, indicating that some counterparts of dark galaxies may have already been observed in the local universe. Most of the z ≈ 0 dwarfs, however, fall in between the bimodal distribution of simulation points, with intermediate tdep,H i of ∼10 Gyr, where there are few systems in our simulation. The sharpness of the division between low and

high tdep,H i systems is likely to be a consequence of the unique value of our H2 floor and our adoption of a single Zfloor ; the latter issue has been emphasized by Tassis et al. (2012) and Krumholz & Dekel (2012). In reality, there is probably some scatter in both the minimum H2 content in low column density gas and the metal content provided by Population III star formation; incorporating this scatter in our simulation would fill in some of the intermediate tdep,H i region. Furthermore, any reduction in the H i content of our dark galaxy halos, which we already know to be necessary from our comparison with the z = 0 ALFALFA H i mass function (see Section 4.4), would move the high tdep,H i points down and to the left, bringing them into better agreement with the local dwarf data. A pure downward shift of these points could be effected by a slight increase in the SFR at a fixed H i content, which may occur through gradual internal or external metal enrichment. 5. SUMMARY AND CONCLUSIONS We have extended the work of Kuhlen et al. (2012) by evolving simulations with a metallicity-dependent, H2 -regulated star formation prescription for an additional one billion years to z = 2.5. We have improved the simulations by imposing a 1% floor in the H2 fraction of cold metal poor gas, where the KMT prescription predicts zero molecular gas content. While the previous simulations exhibited a sharp cutoff in the stellar content below 1010 M , our new simulations result in a more gradual turnover, with star formation becoming effectively stochastic at halo masses below 1010 M . We have shown that this stochasticity reflects the halo formation time, with more massive (at a fixed time) and earlier collapsing (at a fixed mass) halos more easily exceeding the metallicity-dependent surface density threshold (Σgas  5700 M pc−2 (Z/10−3 Z )−0.88 ) for molecular gas formation. If the resulting ∼2 dex scatter in the stellar mass content of 1010 M halos persists to z = 0, it may explain the puzzling dearth of bright dwarf satellite galaxies in the Local Group. In this picture, massive halos that previously have been deemed “too big to fail” (Boylan-Kolchin et al. 2012), would instead just be halos lying in the tails of the f distribution. An interesting implication of our results is that the z ≈ 2.5 universe may be filled with a large population of gas-rich halos with very low stellar content and close to zero ongoing star formation, which we call “dark galaxies.” We define dark galaxies as systems with a neutral gas depletion time of 7

The Astrophysical Journal, 776:34 (9pp), 2013 October 10

Kuhlen, Madau, & Krumholz

mass density of neutral gas at z  2.5 to within a factor of approximately two. We conclude by highlighting some areas where our simulations are challenged by observational data and which demand additional work. The biggest concern is the stellar mass content in halos sufficiently massive (Mh > 1011 M ) to avoid being constrained by the molecular gas transition. The lack of efficient feedback processes in our simulation allows almost half of all gas to be converted into stars, resulting in a stellar mass fraction that is 5–10 times higher than observational constraints based on abundance matching (Moster et al. 2010; Behroozi et al. 2012). A second, possibly related, concern is that, compared with the H i mass function measured by the ALFALFA survey in the local universe (Martin et al. 2010; Haynes et al. 2011), we overproduce by almost one order of magnitude the abundance of MH i < 109 M systems. The majority of this excess is contributed by dark galaxies, for which more effective stellar or SN feedback would not likely help much, given their low SFRs. It remains to be seen whether internal or external processes in the subsequent 11 Gyr between z = 2.5 to z = 0 can reduce the neutral gas content in low mass halos and, if so, whether these processes maintain the stochastic nature of the stellar mass content in the present universe. We thank Naveen Reddy for assistance with the dust attenuation of high-redshift star forming galaxies and Sebastiano Cantalupo, Phil Hopkins, and Eliot Quataert for stimulating discussions. Support for this work was provided by the NSF through grants OIA-1124453 and AST-0955300, by NASA through ATP grants NNX12AF87G, NNX13AB84G, and a Chandra telescope grant, and by the Alfred P. Sloan Foundation. We acknowledge computational support from the NASA Advanced Supercomputing Division, on whose Pleiades supercomputer this simulation was carried out. Figure 6. Comparison of the relation between H i depletion time tdep,H i and MH i with observational data. Top panel: high-redshift data. Simulation points are from the z = 2.5 output. As in Figure 1, we mark halos with SFR = 0 as lower limits at the depletion time corresponding to the minimum resolvable SFR in our simulation. The cyan star is the lower limit from the stack of dark galaxies in Cantalupo et al. (2012, C12) assuming a SFR < 10−2 M yr−1 and their MH i estimates (width of the error bar). The green shaded band indicates the range of lower limits on the depletion times reported by Rafelski et al. (2011, R11) in z ∼ 3 DLAs without a central star formation component. Bottom panel: Low-redshift observational data. Simulation data are the same as in the top panel. The data points (star symbols) are for low surface brightness galaxies (red; Wyder et al. 2009, W09), ALFALFA dwarf galaxies (green; Huang et al. 2012, H12), the THINGS survey (magenta; Walter et al. 2008, W08), and for the newly discovered Leo P dwarf galaxy (cyan; Rhode et al. 2013, R13; Giovanelli et al. 2013, G13). All data, both local and high redshift, have been standardized to the star formation rate–luminosity calibrations recommended by Kennicutt & Evans (2012). (A color version of this figure is available in the online journal.)

REFERENCES Agertz, O., Kravtsov, A. V., Leitner, S. N., & Gnedin, N. Y. 2012, ApJ, 770, 25 Alavi, A., Siana, B., Richard, J., et al. 2013, arXiv:1305.2413 Altay, G., Theuns, T., Schaye, J., Crighton, N. H. M., & Vecchia, C. D. 2011, ApJL, 737, L37 Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2012, ApJ, 770, 57 Ben´ıtez-Llambay, A., Navarro, J. F., Abadi, M. G., et al. 2013, ApJL, 763, L41 Bigiel, F., Leroy, A., Walter, F., et al. 2008, AJ, 136, 2846 Bigiel, F., Leroy, A., Walter, F., et al. 2010, AJ, 140, 1194 Boylan-Kolchin, M., Bullock, J. S., & Kaplinghat, M. 2012, MNRAS, 422, 1203 Bullock, J. S., Kravtsov, A. V., & Weinberg, D. H. 2000, ApJ, 539, 517 Cantalupo, S., Lilly, S. J., & Haehnelt, M. G. 2012, MNRAS, 425, 1992 Christensen, C., Quinn, T., Governato, F., et al. 2012, MNRAS, 425, 3058 Dav´e, R., Katz, N., Oppenheimer, B. D., Kollmeier, J. A., & Weinberg, D. H. 2013, MNRAS, 434, 2645 Dekel, A., & Silk, J. 1986, ApJ, 303, 39 Dijkstra, M., Haiman, Z., Rees, M. J., & Weinberg, D. H. 2004, ApJ, 601, 666 Dutton, A. A., Conroy, C., van den Bosch, F. C., et al. 2011, MNRAS, 416, 322 Efstathiou, G. 1992, MNRAS, 256, 43P Erkal, D., Gnedin, N. Y., & Kravtsov, A. V. 2012, ApJ, 761, 54 Fu, J., Guo, Q., Kauffmann, G., & Krumholz, M. R. 2010, MNRAS, 409, 515 Garrison-Kimmel, S., Rocha, M., Boylan-Kolchin, M., Bullock, J., & Lally, J. 2013, MNRAS, 433, 3539 Genzel, R., Tacconi, L. J., Gracia-Carpio, J., et al. 2010, MNRAS, 407, 2091 Giovanelli, R., Haynes, M. P., Adams, E. A. K., et al. 2013, AJ, 146, 15 Glover, S. C. O., & Clark, P. C. 2012a, MNRAS, 421, 9 Glover, S. C. O., & Clark, P. C. 2012b, MNRAS, 426, 377 Gnedin, N. Y., & Kravtsov, A. V. 2010, ApJ, 714, 287 Gnedin, N. Y., & Kravtsov, A. V. 2011, ApJ, 728, 88 Gnedin, N. Y., Tassis, K., & Kravtsov, A. V. 2009, ApJ, 697, 55 Governato, F., Brook, C., Mayer, L., et al. 2010, Natur, 463, 203 Governato, F., Willman, B., Mayer, L., et al. 2007, MNRAS, 374, 1479

tdep = MH i+H2 /SFR greater than 20 Gyr. In our simulation, we find that the fraction of dark galaxies rises from zero at Mh > 5 × 1010 M to near unity at Mh < 1 × 1010 M . In total, 78% of resolved halos (789 out of 1010) are classified as dark galaxies. Remarkably, the existence of a substantial population of dark galaxies at z = 2.5 appears to be fully compatible with existing observational constraints. In particular our simulated galaxies match the LF determined by Reddy & Steidel (2009), once their luminosities as dimmed by a factor of five (consistent with measurements by Reddy et al. 2012) to account for dust dimming. Our simulation also reproduces the cosmological 8

The Astrophysical Journal, 776:34 (9pp), 2013 October 10

Kuhlen, Madau, & Krumholz

Guedes, J., Callegari, S., Madau, P., & Mayer, L. 2011, ApJ, 742, 76 Haynes, M. P., Giovanelli, R., Martin, A. M., et al. 2011, AJ, 142, 170 Hopkins, P. F., Quataert, E., & Murray, N. 2012, MNRAS, 421, 3522 Huang, S., Haynes, M. P., Giovanelli, R., et al. 2012, AJ, 143, 133 Katz, N. 1992, ApJ, 391, 502 Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531 Knebe, A., & Power, C. 2008, ApJ, 678, 621 Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18 Kravtsov, A. V. 2003, ApJL, 590, L1 Krumholz, M. R. 2012, ApJ, 759, 9 Krumholz, M. R., & Dekel, A. 2012, ApJ, 753, 16 Krumholz, M. R., Dekel, A., & McKee, C. F. 2012, ApJ, 745, 69 Krumholz, M. R., & Gnedin, N. Y. 2011, ApJ, 729, 36 Krumholz, M. R., Leroy, A. K., & McKee, C. F. 2011, ApJ, 731, 25 Krumholz, M. R., McKee, C. F., & Tumlinson, J. 2008, ApJ, 689, 865 Krumholz, M. R., McKee, C. F., & Tumlinson, J. 2009, ApJ, 693, 216 Krumholz, M. R., & Tan, J. C. 2007, ApJ, 654, 304 Kuhlen, M., Krumholz, M. R., Madau, P., Smith, B. D., & Wise, J. 2012, ApJ, 749, 36 Leroy, A. K., Walter, F., Brinks, E., et al. 2008, AJ, 136, 2782 Madau, P., Kuhlen, M., Diemand, J., et al. 2008, ApJL, 689, L41 Madau, P., Pozzetti, L., & Dickinson, M. 1998, ApJ, 498, 106 Martin, A. M., Papastergis, E., Giovanelli, R., et al. 2010, ApJ, 723, 1359 Mashchenko, S., Wadsley, J., & Couchman, H. M. P. 2008, Sci, 319, 174 Mayer, L., Mastropietro, C., Wadsley, J., Stadel, J., & Moore, B. 2006, MNRAS, 369, 1021 McKee, C. F., & Krumholz, M. R. 2010, ApJ, 709, 308 Mo, H. J., Mao, S., & White, S. D. M. 1998, MNRAS, 295, 319

Mori, M., Ferrara, A., & Madau, P. 2002, ApJ, 571, 40 Moster, B. P., Somerville, R. S., Maulbetsch, C., et al. 2010, ApJ, 710, 903 Noterdaeme, P., Petitjean, P., Carithers, W. C., et al. 2012, A&A, 547, L1 Penarrubia, J., Pontzen, A., Walker, M. G., & Koposov, S. E. 2012, ApJL, 759, L42 Pontzen, A., & Governato, F. 2012, MNRAS, 421, 3464 Prochaska, J. X., & Wolfe, A. M. 2009, ApJ, 696, 1543 Rafelski, M., Wolfe, A. M., & Chen, H. 2011, ApJ, 736, 48 Rao, S. M., Turnshek, D. A., & Nestor, D. B. 2006, ApJ, 636, 610 Read, J. I., & Gilmore, G. 2005, MNRAS, 356, 107 Reddy, N., Dickinson, M., Elbaz, D., et al. 2012, ApJ, 744, 154 Reddy, N. A., & Steidel, C. C. 2009, ApJ, 692, 778 Rhode, K. L., Salzer, J. J., Haurberg, N. C., et al. 2013, AJ, 145, 149 Rodr´ıguez, M., & Delgado-Inglada, G. 2011, ApJL, 733, L50 Schruba, A., Leroy, A. K., Walter, F., et al. 2011, AJ, 142, 37 Stinson, G. S., Brook, C., Macci`o, A. V., et al. 2013, MNRAS, 428, 129 Tassis, K., Gnedin, N. Y., & Kravtsov, A. V. 2012, ApJ, 745, 68 Teyssier, R., Pontzen, A., Dubois, Y., & Read, J. I. 2013, MNRAS, 429, 3068 Thoul, A. A., & Weinberg, D. H. 1996, ApJ, 465, 608 Walter, F., Brinks, E., de Blok, W. J. G., et al. 2008, AJ, 136, 2563 Wise, J. H., Turk, M. J., Norman, M. L., & Abel, T. 2012, ApJ, 745, 50 Wolfire, M. G., McKee, C. F., Hollenbach, D., & Tielens, A. G. G. M. 2003, ApJ, 587, 278 Wolfire, M. G., Tielens, A. G. G. M., Hollenbach, D., & Kaufman, M. J. 2008, ApJ, 680, 384 Wyder, T. K., Martin, D. C., Barlow, T. A., et al. 2009, ApJ, 696, 1834 Zolotov, A., Brooks, A. M., Willman, B., et al. 2012, ApJ, 761, 71

9