We present a series of highresolution cosmological simulationsSUP1SUP of galaxy formation to z  0 spanning halo masses 10SUP8SUP10SUP13SUP MSUBSUB and stellar masses 10SUP4SUP10SUP11SUP MSUBSUB. Our simulations include fully explicit treatment of the multiphase interstellar medium and stellar feedback. The stellar feedback inputs energy momentum mass and metal fluxes are taken directly from stellar population models. These sources of feedback with zero adjusted parameters reproduce the observed relation between stellar and halo mass up to MSUBhaloSUB  10SUP12SUP MSUBSUB. We predict weak redshift evolution in the MSUBSUBMSUBhaloSUB relation consistent with current constraints to z gt 6. We find that the MSUBSUBMSUBhaloSUB relation is insensitive to numerical details but is sensitive to feedback physics. Simulations with only supernova feedback fail to reproduce observed stellar masses particularly in dwarf and highredshift galaxies radiative feedback photoheating and radiation pressure is necessary to destroy giant molecular clouds and enable efficient coupling of later supernovae to the gas. Star formation rates SFRs agree well with the observed Kennicutt relation at all redshifts. The galaxyaveraged Kennicutt relation is very different from the numerically imposed law for converting gas into stars and is determined by selfregulation via stellar feedback. Feedback reduces SFRs and produces reservoirs of gas that lead to rising latetime star formation histories significantly different from halo accretion histories. Feedback also produces large shorttimescale variability in galactic SFRs especially in dwarfs. These properties are not captured by common subgrid wind models.
['stars: formation', 'galaxies: active', 'galaxies: evolution', 'galaxies: formation', 'cosmology: theory', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'Astrophysics - Astrophysics of Galaxies']
Galaxies on FIRE Feedback In Realistic Environments stellar feedback explains cosmologically inefficient star formation
{'Galaxies on FIRE (Feedback In Realistic Environments): Stellar Feedback Explains Cosmologically Inefficient Star Formation': 'Philip F. Hopkins ∗ 1 , 2 , Dušan Kereš 3 , José Oñorbe 4 , Claude-André Faucher-Giguère 2 , 5 , Eliot Quataert 2 , Norman Murray 6 , 7 , & James S. Bullock 4  \n1 TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA  \n- 2 Department of Astronomy and Theoretical Astrophysics Center, University of California Berkeley, Berkeley, CA 94720\n- 3 Department of Physics, Center for Astrophysics and Space Science, University of California at San Diego, 9500 Gilman Drive, La Jolla, CA 92093\n- 4 Department of Physics and Astronomy, University of California at Irvine, Irvine, CA 92697, USA\n- 5 Department of Physics and Astronomy and CIERA, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA\n- 6 Canadian Institute for Theoretical Astrophysics, 60 St. George Street, University of Toronto, ON M5S 3H8, Canada\n- 7 Canada Research Chair in Astrophysics  \nSubmitted to MNRAS, November, 2013', 'ABSTRACT': "We present a series of high-resolution cosmological simulations 1 of galaxy formation to z = 0, spanning halo masses ∼ 10 8 -10 13 M glyph[circledot] , and stellar masses ∼ 10 4 -10 11 M glyph[circledot] . Our simulations include fully explicit treatment of the multi-phase ISM & stellar feedback. The stellar feedback inputs (energy, momentum, mass, and metal fluxes) are taken directly from stellar population models. These sources of feedback, with zero adjusted parameters, reproduce the observed relation between stellar and halo mass up to M halo ∼ 10 12 M glyph[circledot] . We predict weak redshift evolution in the M ∗ -M halo relation, consistent with current constraints to z > 6. We find that the M ∗ -M halo relation is insensitive to numerical details, but is sensitive to feedback physics. Simulations with only supernova feedback fail to reproduce observed stellar masses, particularly in dwarf and high-redshift galaxies: radiative feedback (photo-heating and radiation pressure) is necessary to destroy GMCs and enable efficient coupling of later supernovae to the gas. Star formation rates agree well with the observed Kennicutt relation at all redshifts. The galaxy-averaged Kennicutt relation is very different from the numerically imposed law for converting gas into stars, and is determined by self-regulation via stellar feedback. Feedback reduces star formation rates and produces reservoirs of gas that lead to rising late-time star formation histories, significantly different from halo accretion histories. Feedback also produces large short-timescale variability in galactic SFRs, especially in dwarfs. These properties are not captured by common 'sub-grid' wind models.  \nKey words: galaxies: formation - galaxies: evolution - galaxies: active - stars: formation - cosmology: theory", '1 INTRODUCTION': "It is well-known that feedback from stars is a critical, yet poorlyunderstood, component of galaxy formation. Within galaxies, star formation is observed to be inefficient in both an instantaneous and an integral sense.  \nInstantaneously, the Kennicutt-Schmidt (KS) relation implies gas consumption timescales of ∼ 50 dynamical times (Kennicutt 1998), while the total fraction of GMC mass converted into stars is only a few percent (Zuckerman & Evans 1974; Williams & McKee 1997; Evans 1999; Evans et al. 2009). Without strong stellar feedback, however, gas inside galaxies cools efficiently and collapses on a dynamical time, predicting order-unity star formation efficiencies on all scales (Hopkins et al. 2011; Tasker 2011; Bournaud et al. 2010; Dobbs et al. 2011; Krumholz et al. 2011; Harper-Clark & Murray 2011).  \nIn an integral sense, without strong stellar feedback, gas in cosmological models cools rapidly and inevitably turns into stars, predicting galaxies with far larger masses than are observed (e.g. Katz et al. 1996; Somerville & Primack 1999; Cole et al. 2000; Springel & Hernquist 2003b; Kereš et al. 2009, and references therein). Decreasing the instantaneous star formation efficiency does not eliminate this integral problem: the amount of baryons  \n- ∗ E-mail:phopkins@caltech.edu  \nin real galactic disks is much lower than that predicted in models absent strong feedback (essentially, the Universal baryon budget; see White & Frenk 1991; Kereš et al. 2009). Constraints from intergalactic medium (IGM) enrichment require that many of those baryons must have entered galaxy halos and disks at some point to be enriched, before being expelled (Aguirre et al. 2001; Pettini et al. 2003; Songaila 2005; Martin et al. 2010). Galactic superwinds with mass-loading ˙ M wind of many times the star formation rate (SFR) are therefore generally required to reproduce observed galaxy properties (e.g. Oppenheimer & Davé 2006). Such winds have been observed ubiquitously in local and high-redshift starforming galaxies (Martin 1999, 2006; Heckman et al. 2000; Newman et al. 2012; Sato et al. 2009; Chen et al. 2010; Steidel et al. 2010; Coil et al. 2011).  \nHowever, until recently, numerical simulations have been unable to produce winds with large-mass loading factors from an a priori model (let alone the correct scalings of wind mass-loading with galaxy mass or other properties), nor to simultaneously predict the instantaneous inefficiency of star formation within galaxies. This is particularly true of models which invoke only energetic feedback via supernovae (SNe), which is efficiently radiated in the dense gas where star formation actually occurs (see e.g. Guo et al. 2010; Powell et al. 2011; Brook et al. 2011; Nagamine 2010; Bournaud et al. 2011, and references therein). More recent simulations, using higher resolution and invoking stronger feedback prescrip-  \n<!-- image -->  \nFigure 1. Gas in a representative simulation of a Milky Way-mass halo ( m12i in Table 1). Image shows the projected gas density, log-weighted ( ∼ 4dex stretch). Magenta shows cold molecular/atomic gas ( T < 1000K). Green shows warm ionized gas (10 4 glyph[lessorsimilar] T glyph[lessorsimilar] 10 5 K). Red shows hot gas ( T glyph[greaterorsimilar] 10 6 K). 2 Each image shows a box centered on the main galaxy. Left: Box 200kpc (physical) on a side at high redshift. The galaxy has undergone a violent starburst, leading to strong outflows of hot and warm gas that have blown away much of the surrounding IGM (even outside the galaxy). Note that the 'filamentary' structure of cool gas in the IGM is clearly affected by the outflows. Right: Near present-day, with a ∼ 50kpc box. A more relaxed, well-ordered disk has formed, with molecular gas tracing spiral structure, and a halo enriched by diffuse hot outflows.  \n<!-- image -->  \nFigure 2. Stars in the m12i simulation at z ∼ 0, in a box 50kpc on a side near present-time. Image is a mock u / g / r composite. The disk is approximately face-on, and the spiral structure is visible. (The image uses STARBURST99 to determine the SED of each star particle given its known age and metallicity, then ray-traces the line-of-sight flux following Hopkins et al. (2005), attenuating with a MW-like reddening curve with constant dust-to-metals ratio for the abundances at each point.)  \n<!-- image -->  \ntions, have seen strong winds, but have generally found it necessary to include simplified prescriptions for 'turning off cooling' in the SNe-heated gas and/or include some adjustable parameters representing 'pre-SNe' feedback (see Governato et al. 2010; Macciò et al. 2012; Teyssier et al. 2013; Stinson et al. 2013; Agertz et al. 2013). This is physically motivated since feedback processes other than SNe - protostellar jets, HII photoionization, stellar winds, and radiation pressure - both occur and are critical in suppressing star formation in dense gas, as well as 'pre-processing' gas prior to SNe explosions so that SNe occur at densities where thermal heating can have much larger effects (Evans et al. 2009; Hopkins et al. 2011; Tasker 2011; Lopez et al. 2011; Stinson et al. 2013; Kannan et al. 2013).  \nAnd in fact, there have been many studies with enormously higher resolution (enough to evolve each star explicitly) and a full treatment of the radiation-magnetohydrodynamics and time dependence of these multiple feedback mechanisms. Because of computational limitations, however, these have necessarily been restricted to very small systems, either single molecular clouds/star clusters (e.g. Krumholz et al. 2007, 2011; Offner et al. 2009, 2011; HarperClark & Murray 2011; Bate 2012), or the 'first stars' (e.g. Wise et al. 2012; Pawlik et al. 2013; Muratov et al. 2013). But these studies, without exception , have found that the non-linear interaction of the feedback mechanisms above - especially the dual roles of HII photoionization and radiation pressure in concert with SNe - is absolutely critical to explain the generation of large local outflows, the self-regulation of star formation, and the shape of the stellar initial mass function.  \nDespite these breakthroughs, given limited resolution and the complexity of the baryonic physics, many cosmological models have treated galactic wind generation and the inefficiency of star  \n<!-- image -->  \nFigure 3. Gas, as Fig. 1, for a dwarf galaxy ( m10 in Table 1). Top: 40kpc (physical) box, at high redshift. Bottom: 20kpc box at intermediate redshift. Strong outflows are still present, though they are more spherical, because the galaxy halo is itself small and embedded within a much larger filament.  \n<!-- image -->  \nformation in a tuneable, 'sub-grid' fashion. This is not to say that the models have not tremendously improved our understanding of galaxy formation! They have demonstrated that stellar feedback can plausibly lead to (globally) inefficient star formation, constrained the parameter space of allowed feedback models, made predictions for the critical role of outflows and recycling in enriching the IGM, provided possible baryonic solutions to apparent dark matter 'problems' (e.g. Pontzen & Governato 2012), demonstrated the need for 'early' feedback from radiative mechanisms beyond SNe alone, and generally created the framework for our interpretation of observations. However, with wind models often relying on adjustable parameters, the integrated efficiency of star formation in galaxies is to some extent tuned 'by hand' and the predictive power is inherently limited. This is particularly true for studies of gas in the circum-galactic medium (CGM), a current area of much observational progress - measurements of the CGM are sensitive  \nto the phase structure of the gas, which is not faithfully represented in models which simply 'turn off' hydrodynamics or cooling, or mimic strong feedback via pure thermal energy injection or 'particle kicks' (see e.g. Hummels et al. 2013, for an explicit demonstration of this).  \nAccurate treatment of star formation and galactic winds ultimately requires realistic treatment of the stellar feedback processes that maintain the multi-phase ISM. Motivated by this philosophy (and building on the studies with single-star resolution), in Hopkins et al. (2011) (Paper I) and Hopkins et al. (2012d) (Paper II), we developed a new set of numerical models to follow stellar feedback on scales from sub-GMC star-forming regions through galaxies. These simulations include the energy, momentum, mass, and metal fluxes from stellar radiation pressure, HII photo-ionization and photo-electric heating, SNe Types I & II, and stellar winds (O-star and AGB). Critically, the feedback is directly tied to the young stars, with the energetics and time-dependence taken from stellar evolution models. In our previous work, we showed, in isolated galaxy simulations, that these mechanisms produce a quasisteady ISM in which GMCs form and disperse rapidly, with phase structure, turbulence, and disk and GMC properties in good agreement with observations (for various comparisons, see Narayanan & Hopkins 2012; Hopkins et al. 2012b,a, 2013c). In Paper I, Hopkins et al. (2013d), and Hopkins et al. (2013a) we showed that this leads naturally to 'instantaneously' inefficient SF (predicting the KS-law), regulated self-consistently by feedback and independent of the numerical prescription for star formation in very dense gas. In Hopkins et al. (2012c) (Paper III) and Hopkins et al. (2013b) we showed that the same feedback models reproduce the galactic winds invoked in previous semi-analytic and cosmological simulations, and that the combination of multiple feedback mechanisms is critical to produce massive, multi-phase galactic winds.  \nHowever, our simulations have thus far been limited to idealized studies of isolated galaxies and galaxy mergers. These previous calculations thus cannot follow accretion from or interaction of outflows with the IGM, realistic galaxy merger histories, and many other important processes. In this paper, the first of a series, we present the FIRE (Feedback In Realistic Environments) simulations: 1 a suite of fully cosmological 'zoom-in' simulations developed to study the role of feedback in galaxy formation. To test the models and understand feedback in a wide range of environments, we study a wide range in galaxy halo and stellar mass (as opposed to focusing just on MW-like systems), and follow evolution fully to z = 0. Our suite of calculations includes several of the highestresolution galaxy formation simulations that have been run to z = 0. Our simulations utilize a significantly improved numerical implementation of SPH (which has resolved historical discrepancies with grid codes), as well as the full physical models for feedback and ISM physics introduced and tested in Paper I-Paper III. Here, we explore the consequences of stellar feedback for the inefficiency of star formation, perhaps the most basic consequence of stellar feedback for galaxy formation. In companion papers, we will investigate the properties of outflows and their interactions with the IGM, the effect of those outflows on dark matter structure, the differences between numerical methods in treating feedback, the role  \nof feedback in determining galaxy structure, and many other open questions.  \nIn § 2-4, we describe our methodology. § 2 describes the initial conditions for the simulations; § 3 outlines the implementation of the key baryonic physics of cooling, star formation, and feedback (a much more detailed description is given in Appendix A); § 4 briefly describes the improvements in the numerical method compared to past work (again, more details are in Appendix B). And in Appendix C we test and compare these algorithms with higherresolution simulations of isolated (non-cosmological) galaxies.  \nWe describe our results in § 5. We examine the predicted galaxy stellar masses (§ 5.1), and how this depends on both numerical algorithms (§ 5.3) and feedback physics (§ 5.4), as well as how it compares to previous theoretical work (§ 5.5). We show that the treatment of feedback physics overwhelmingly dominates these results, and discuss the distinct roles of multiple independent feedback mechanisms. We also explore the predictions for the Kennicutt-Schmidt relation (§ 5.6), the shape of galaxy star formation histories (§ 5.7), the star formation 'main sequence' (§ 5.8), and the 'burstiness' of star formation (§ 5.9). We summarize our important conclusions and discuss future work in § 6.", '2 INITIAL CONDITIONS & GALAXY PROPERTIES': "The simulations presented here are a series of fully cosmological 'zoom-in' simulations of galaxy formation; some images of the gas and stars in representative stages are shown in Figs. 1-3. 2 The technique is well-studied; briefly, a large cosmological box is simulated at low resolution to z = 0, and then the mass within and around halos of interest at that time is identified, traced back to the starting redshift, and the Lagrangian region containing this mass is re-initialized at much higher resolution (with gas added) for the ultimate simulation (Porter 1985; Katz & White 1993).  \nWe consider a series of systems with different masses. Table 1 describes the initial conditions. All simulations begin at redshifts ∼ 100 -125, with fluctuations evolved using perturbation theory up to that point. 3  \nThe specific halos we re-simulate are chosen to represent a broad mass range and be 'typical' in most properties (e.g. sizes, formation times, and merger histories) relative to other halos of the same z = 0 mass. The simulations m09 and m10 are constructed using the methods from Onorbe et al. (2013); they are isolated dwarfs. Simulations m11 , m12q , m12i , and m13 are chosen to match a subset of initial conditions from the AGORA project (Kim et al. 2013), which will enable future comparisons with a wide range of different codes. These are chosen to be somewhat quiescent merger histories, but lie well within the typical scatter in such histories at each mass (and each has several major mergers). Simulation m12v , for contrast, is chosen to have a relatively violent merger history (several major mergers since z ∼ 2), and is based on the initial conditions studied in Kereš & Hernquist (2009) and Faucher-Giguère &Kereš (2011).  \nFigure 4. Galaxy stellar mass-halo mass relation at z = 0. Top : M ∗ ( M halo ) . Bottom: M ∗ relative to the Universal baryon budget of the halo ( fb M halo). Each simulation (points) from Table 1 is shown; large point denotes the most massive halo in each box. We compare the relation if all baryons became stars ( M ∗ = fb M halo; dotted) and the observationally inferred relationship as determined in Moster et al. (2013, magenta) and Behroozi et al. (2012, cyan) (dashed lines denote extrapolation beyond the observed range). 9 The agreement with observations is excellent at M halo glyph[lessorsimilar] 10 13 M glyph[circledot] , including dwarf though MW-mass galaxies. We stress that there are zero adjusted parameters here: stellar feedback, with known mechanisms taken from stellar population models, is sufficient to explain galaxy stellar masses at/below ∼ L ∗ .  \n<!-- image -->  \n/circledot  \nIn each case, the resolution is scaled with the simulated mass, so as to achieve the optimal possible force and mass resolution. It is correspondingly possible to resolve much smaller structures in the low-mass galaxies. The critical point is that in all our simulations with mass < 10 13 M glyph[circledot] , we resolve the Jeans mass/length of gas in the galaxies, corresponding to the size/mass of massive molecular cloud complexes. This is necessary to resolve a genuine multi-phase ISM and for our ISM feedback physics to be meaningful. Fortunately, because most of the mass and star formation in GMCs in both observations (Evans 1999; Blitz & Rosolowsky 2005) and simulated systems (Paper II) is concentrated in the most massive GMCs, the resolution studies in Paper I-Paper II confirm that resolving small molecular clouds makes little difference. We refer interested readers to Paper II for a detailed discussion of the scales that must be resolved for feedback to operate appropriately, but note here that all our simulations are designed to be approxi-  \nTable 1. Simulation Initial Conditions  \nParameters describing the initial conditions for our simulations (units are physical):  \n(1)  \nName: Simulation designation.  \n- (2) M 0 halo : Approximate mass of the z = 0 'main' halo (most massive halo in the high-resolution region). (3) mb : Initial baryonic (gas and star) particle mass in the high-resolution region, in our highest-resolution simulations.  \n(4) glyph[epsilon1] b : Minimum baryonic gravity/force softening (minimum SPH smoothing lengths are comparable or smaller). Recall, force softenings are adaptive (mass resolution is fixed); for more details see Appendix B. (3) mdm : Dark matter particle mass in the high-resolution region, in our highest-resolution simulations. (4) glyph[epsilon1] dm : Minimum dark matter force softening (fixed in physical units at all redshifts).  \nFigure 5. M ∗ -M halo relation as Fig. 4 (points follow the legend therein), at different redshifts. Observational constraints are also shown at each redshift (each pair of lines shows the ± 1 σ fit to the observations at that redshift). With no tuned parameters, the simulations predict M ∗ -M halo and, by extension, the stellar MF and galaxy clustering, at all z . Redshift evolution in M ∗ -M halo is weak, with the sense that low-mass dwarfs become higherM ∗ , leading to a steeper faint-end galaxy MF, in agreement with constraints from reionization (see Kuhlen & Faucher-Giguère 2012, and references therein).  \n<!-- image -->  \n/circledot  \n/circledot  \n/circledot  \nmately comparable to the 'high-resolution' simulations of isolated galaxies and the ISM in Paper I-Paper II, within the range of resolution where the results in those studies (star formation rates, wind outflow rates, GMC lifetimes, etc.) were numerically converged (unfortunately, it is not possible to evolve cosmological simulations to z = 0 with the 'ultra-high' sub-pc resolution therein).  \nIn terms of the Jeans mass/length of the galaxies, our resolution is broadly comparable between different simulations. Our worst resolution in units of the Jeans length/mass occurs in the more massive galaxies at late times, when they are relatively gas poor, and so (despite the large total galaxy mass) the Jeans length  \ncan become relatively small. 4 Every galaxy identified in this paper contains at least glyph[greatermuch] 10 5 bound particles.  \nWe adopt a 'standard' flat Λ CDM cosmology with h ≈ 0 . 7, Ω M = 1 -Ω Λ ≈ 0 . 27, and Ω b ≈ 0 . 046 for all runs. 5", '3 BARYONIC PHYSICS': 'The simulations here use the physical models for star formation and stellar feedback developed and presented in a series of papers studying isolated galaxies (Hopkins et al. 2012c, 2013a, 2012b,a), adapted for fully cosmological simulations. We summarize their properties below, but refer to Appendix A for a more detailed explanation and list of improvements. Readers interested in further details (including resolution studies and a range of tests of the specific numerical methodology) should see Paper I & Paper II.', '3.1 Cooling': 'Gas follows an ionized+atomic+molecular cooling curve from 10 -10 10 K, including metallicity-dependent fine-structure and molecular cooling at low temperatures, and high-temperature ( glyph[greaterorsimilar] 10 4 K) metal-line cooling followed species-by-species for 11 separately tracked species. At all times, we tabulate the appropriate ionization states and cooling rates from a compilation of CLOUDY runs, including the effect of the photo-ionizing background, accounting for gas self-shielding. Photo-ionization and photo-electric heating from local sources are accounted for as described below.', '3.2 Star Formation': "Star formation is allowed only in dense, molecular, self-gravitating regions above n > n crit ( n crit = 100cm -3 for our primary runs, but we also tested from ∼ 10 -1000cm -3 ). This threshold is much higher than that adopted in most 'zoom-in' simulations of galaxy formation (the high value allows us to capture highly clustered star formation). We follow Krumholz & Gnedin (2011) to calculate the molecular fraction f H 2 in dense gas as a function of local column density and metallicity, and allow SF only from molecular gas. We also follow Hopkins et al. (2013d) and restrict star formation to gas which is locally self-gravitating, i.e. has α ≡ δ v 2 δ r / Gm gas ( < δ r ) < 1 on the smallest available scale ( δ r being our force softening or smoothing length). This forms stars at a rate ˙ ρ ∗ = ρ mol / t ff (i.e. 100% efficiency per free-fall time); so that the galaxy and even kpcscale star formation efficiency is not set by hand, but regulated by feedback (typically at much lower values). Because of this, in Paper I, Paper II, and Hopkins et al. (2013d) we show that the galaxy structure and SFR are basically independent of the small-scale SF law, density threshold (provided it is high), and treatment of molecular chemistry.", '3.3 Stellar Feedback': 'Once stars form, their feedback effects are included from several sources. Every star particle is treated as a single stellar population, with a known age, metallicity, and mass. Then all feedback quantities (the stellar luminosity, spectral shape, SNe rates, stellar wind mechanical luminosities, metal yields, etc.) are tabulated as a function of time directly from the stellar population models in STARBURST99, assuming a Kroupa (2002) IMF.  \n- (1) Radiation Pressure: Gas illuminated by stars feels a momentum flux ˙ P rad ≈ ( 1 -exp ( -τ UV / optical ))( 1 + τ IR ) L incident / c along the optical depth gradient, where 1 + τ IR = 1 +Σ gas κ IR accounts for the absorption of the initial UV/optical flux and multiple scatterings of the re-emitted IR flux if the region between star and  \nfor numerical comparisons, they feature very small cosmological parameter differences. These are percent-level, smaller than the observational uncertainties in the relevant quantities (Planck Collaboration et al. 2013) and produce negligible effects compared to differences between randomly chosen halos.  \ngas particle is optically thick in the IR (see Appendix A). We assume that the opacities scale linearly with gas metallicity. 6  \n- (2) Supernovae: Wetabulate the SNe Type-I and Type-II rates from Mannucci et al. (2006) and STARBURST99, respectively, as a function of age and metallicity for all star particles and stochastically determine at each timestep if an individual SNe occurs. If so, the appropriate mechanical luminosity and ejecta momentum is injected as thermal energy and radial momentum in the gas within a smoothing length of the star particle, along with the relevant mass and metal yield (for all followed species). When the Sedov-Taylor phase is not fully resolved, we account for the work done by hot gas inside the unresolved cooling radius (converting the appropriate fraction of the SNe energy into momentum). We discuss this in detail in Appendix A, but emphasize that it is particularly important that SNe momentum not be neglected in massive halos whose mass resolution ∼ 10 4 M glyph[circledot] is much larger than the ejecta mass of a single SNe.\n- (3) Stellar Winds: Similarly, stellar winds are assumed to shock locally and so we inject the appropriate tabulated mechanical power L ( t , Z ) , wind momentum, mass, and metal yields, as a continuous function of age and metallicity into the gas within a smoothing length of the star particles. The integrated mass fraction recycled in winds (including both fast winds from young stars and slow AGB winds) and SNe is ∼ 0 . 3.\n- (4) Photo-Ionization and Photo-Electric Heating: Knowing the ionizing photon flux from each star particle, we ionize each neighboring neutral gas particle (provided there are sufficient photons, given the gas density, metallicity, and prior ionization state), moving outwards until the photon budget is exhausted; this alters the heating and cooling rates appropriately. The UV fluxes are also used to determine photo-electric heating rates following Wolfire et al. (1995).  \nExtensive numerical tests of the feedback models are presented in Paper II.', '4 SIMULATION NUMERICAL DETAILS': "All simulations are run using a newly developed version of TreeSPH which we refer to as 'P-SPH' (Hopkins 2013), in the code GIZMO. 7 This adopts the Lagrangian 'pressure-entropy' formulation of the SPH equations developed in Hopkins (2013); this eliminates the major differences between SPH, moving mesh, and grid (adaptive mesh) codes, and resolves the well-known issues with fluid mixing instabilities in previously-used forms of SPH (e.g. Agertz et al. 2007; Sijacki et al. 2012). The gravity solver is a heavily modified version of the GADGET-3 code (Springel 2005); but GIZMO also includes substantial improvements in the artificial viscosity, entropy diffusion, adaptive timestepping, smoothing kernel, and gravitational softening algorithm, as compared to the 'previous generation.' These are all described in detail in Appendix B.  \n6 There has been some debate in the literature regarding whether or not the full τ IR 'boost' applies to the infrared radiation pressure when τ IR glyph[greatermuch] 1 (see e.g. Krumholz & Thompson 2012, but also Kuiper et al. 2012 and Davis et al. (2014), who find much stronger effects in the infrared). We have considered alternatives, discussed in Paper I. However, in the simulations here we never resolve the extremely high densities where τ IR glyph[greaterorsimilar] 1 (where this distinction is important), and so if anything are under -estimating the IR radiation pressure, even compared to the most conservative studies.  \nwww.tapir.caltech.edu/~phopkins/Site/GIZMO.html  \nWe emphasize that our version of SPH has been tested extensively and found to give good agreement with analytic solutions as well as well-tested grid codes on a broad suite of test problems. Many of these are presented in Hopkins (2013). This includes Sod shock tubes; Sedov blastwaves; wind tunnel tests (radiative and adiabatic, up to Mach ∼ 10 4 ); linear sound wave propagation; oscillating polytropes; hydrostatic equilibrium 'deformation'/surface tension tests (Saitoh & Makino 2013); Kelvin-Helmholtz and Rayleigh-Taylor instabilities; the 'blob test' (Agertz et al. 2007); super-sonic and sub-sonic turbulence tests (from Mach ∼ 0 . 1 -10 3 ); Keplerian gas ring, disk shear, and shearing shock tests (Cullen & Dehnen 2010); the Evrard test; the Gresho-Chan vortex; spherical collapse tests; and non-linear galaxy formation tests. For additional tests showing the improvements relative to previousgeneration SPH, see Hu et al. (2014). Since it is critical for the problems addressed here that a code be able to handle high dynamic range situations, the numerical method and parameters such as SPH 'neighbor number' were not modified for these tests individually, but are similar to what we use in our production runs in this paper.  \nIn Appendix B, we note that we have explicitly tested many of the purely numerical elements of the gravity and hydrodynamic solvers in the simulations shown here: for example, whether to use adaptive or fixed gravitational softenings, the choice of SPH smoothing kernel, and the timestepping algorithm. However we do not discuss these in the main text because they produce extremely small ( glyph[lessorsimilar] 10%-level) differences in the quantities plotted in this paper.", '5.1 Galaxy Masses as a Function of Redshift': "Fig. 4 plots the z = 0 stellar mass-halo mass relation for our main set of simulations from Table 1 (highest-resolution, with all physics enabled). Note that although each high-resolution region at z = 0 contains one 'primary' halo (the focus of that region), there are several smaller-mass, independent halos also in that region. We therefore identify and plot all such halos. 8 We exclude those halos that are outside the high-resolution region (more than 1% masscontaminated by low-resolution particles; although varying this between 0 . 5 -10%makes little difference to our comparisons here) or insufficiently resolved ( < 0 . 01 times the primary halo mass, or with < 10 5 dark matter particles). We also exclude subhalos/satellite galaxies.  \nThe known sources of stellar feedback we include, with no adjustment, automatically reproduce a relation between galaxy stellar  \n8 We use the HOP halo finder (Eisenstein & Hut 1998) to automatically identify halos (which combines an iterative overdensity identification with a saddle density threshold criterion to merge subhalos and overlapping halos). Halo masses are defined as the mass within a spherical aperture about the density maximum with mean density > 200 times the critical density at each redshift (this is chosen to be similar to the choice used in abundancematching models, which define the observations to which we compare). Stellar mass plotted is the total stellar mass within 20kpc of the center of the central galaxy in the halo (we do not include satellite galaxy masses). However, we have compared with the results of a basic friends-of-finds routine or simple by-eye identification, and find that for the results here (focused on simple, integral halo quantities, and ignoring subhalos), this makes no significant difference. Likewise defining the mass within ∼ 0 . 1 R vir instead of 20kpc makes no significant difference.  \nand halo mass consistent with the observations 9 from M halo ∼ 10 7 -10 13 M glyph[circledot] . Specifically, the distribution of points for all M halo ≤ 10 12 M glyph[circledot] is statistically consistent (in a χ 2 sense) with having been drawn from the Moster et al. (2013) curve; if we allow for the observed scatter ( ∼ 0 . 15dex at ∼ 10 12 M glyph[circledot] , the width between the plotted lines, increasing to ∼ 0 . 3dex at the lowest observed masses) then all our primary galaxies lie within the 2 σ scatter. 10  \nDespite the fact that this relation implies a non-uniform (and even non-monotonic) efficiency of star formation as a function of galaxy mass, we do not need to invoke different physics or distinct parameters at different masses. This is particularly impressive at low masses, where the integrated stellar mass must be suppressed by factors of ∼ 1000 relative to the Universal baryon fraction. Unfortunately, at high masses ( > 10 13 M glyph[circledot] ), the large Lagrangian regions (hence large number of required particles) limit the resolution we can achieve; we have experimented with some low-resolution test runs which appear to produce overly massive galaxies, but higher-resolution studies are required to determine if that owes to a need for additional physics or simply poor numerical resolution.  \nInterestingly, the scatter in M ∗ at fixed M halo may decrease weakly with mass, from ∼ 0 . 5 dex in dwarf galaxies ( M halo glyph[lessorsimilar] 10 10 M glyph[circledot] ) to ∼ 0 . 1 -0 . 2 dex in massive ( ∼ L ∗ ) galaxies. But given the limited number of halos we study here, further investigation allowing more diverse merger/growth histories is needed.  \nFig. 5 shows the M ∗ -M halo relation at various redshifts. At each z , we compare with observationally constrained estimates of the M ∗ -M halo relation. Implicitly, if they agree in M ∗ ( M halo ) , our models are consistent with the observed stellar MF (given, of course, the limited statistics by which we are 'sampling' the MF). At high redshifts, the halos we simulate are of course lower-mass, so eventually we have no high-mass galaxies; this limits the extent to which our results can be compared to observations above z ∼ 2.", '5.2 Other (Basic) Galaxy Properties': "We wish to focus here on galaxy masses and star formation histories. Companion papers (in preparation) will examine the galaxy morphologies and other observables in more detail. It is important, when studying those properties, to construct a meaningful comparison (e.g. using the same methods and wavelengths observed), and this is non-trivial. Moreover, it is by no means obvious that these properties are as robust to numerical details as the galaxy stellar masses (discussed further below), and it is completely outside the scope of this paper to fairly explore those dependencies.  \nThat said, we can briefly note the basic properties of the specific simulations in Table 1 at z = 0, with the caveat that these may not be robust to changes in either the initial conditions (the particular halo simulated) or our numerical methods. Morphologically, at z = 0, run m09 resembles an ultrafaint dwarf; m10 a thick, but rotating dwarf irregular; and m11 a more 'fluffy' dwarf spheroidal.  \n10 We can also fit the points here to the same power-law functional form used empirically: if we do so, the best-fit slope and normalization are both within 0 . 5 σ of the fit to observations in Moster et al. (2013) (the error bar is dominated by the small-number statistics in our halo sampling). The simulations with M halo glyph[lessmuch] 10 10 M glyph[circledot] are statistically inconsistent with the extrapolation of the flatter slope from Behroozi et al. (2012), but this is entirely below the region actually observed, where Behroozi et al. (2012) and Moster et al. (2013) agree well, and there the simulations do not significantly 'prefer' either fit.  \nRuns m12v , m12q , m12i produce bulge+disk systems, with m12v showing a prominent bulge at all times z glyph[lessorsimilar] 2; m12q is more diskdominated until a late major merger at z < 0 . 5 destroys the disk; and run m12i produces a stellar disk with little bulge. Run m13 is totally bulge-dominated. Each galaxy has an approximately flat rotation curve outside of the central couple kpc; those with M ∗ < 10 10 slowly rise with radius to V max, and the more massive systems are flat to within the central ∼ kpc, except for m12v (where the compact bulge leads to a central-kpc spike at ∼ 250kms -1 ). The galaxy sizes, measured as the half-stellar mass effective radii, are ( 0 . 3 , 0 . 52 , 3 . 5 , 2 . 8 , 3 . 2 , 4 . 2 , 4 . 8 ) kpc for ( m09, m10, m11, m12v, m12q, m12i, m13 ), consistent with the observed stellar size-mass relation (Shen et al. 2003; Wolf et al. 2010).", '5.3 (Lack of) Dependence on Numerical Methods': "In Fig. 6 we investigate how the M ∗ -M halo relation depends on numerical parameters and feedback. First we repeat Fig. 4 for simulations with different purely numerical parameters. These can and do, indeed, have significant quantitative effects - they can easily shift the predicted stellar masses by factors ∼ 2 -3. However, we stress that they do not qualitatively change our conclusions.  \nModest changes in resolution (our 'low-resolution' runs correspond to one power of two step in spatial resolution, and a corresponding factor of 2 3 = 8 change in mass resolution) lead to significant, but not order-of-magnitude, changes in M ∗ : generally we obtain larger M ∗ by factors of ∼ 1 . 5 at high masses ( M halo glyph[greaterorsimilar] 10 11 M glyph[circledot] ) and ∼ 2 -3 at the lowest masses ( M halo < 10 10 M glyph[circledot] ) at lower resolution, owing to a combination of (a) artificially enhanced mixing and thus cooling of diffuse gas, since ISM phases are less wellresolved, and (b) the fact that the coupling of feedback energy and momentum is necessarily spread over larger mass elements. If we downgrade our resolution more substantially - by a factor of ∼ 100 in mass, or > 10 in spatial scale (i.e. using the > 100pc spatial resolution which is typical of most previous cosmological simulations), the results diverge more substantially: galaxy masses at z ∼ 0 are a factor of ∼ 3 -5 higher at high masses and ∼ 10 higher at low masses. This makes sense, because at that resolution, we simply cannot meaningfully resolve even the most massive structures in the ISM. 11  \nSome of our numerical tests are not plotted here because their effects are not significant. We have, for example, re-run several simulations with twice and five times larger dark matter softening lengths (same baryonic softening); using or de-activating adaptive gravitational softenings (which ensure there are always ∼ 100 neighbor particles in the softening kernel); varying the number of SPH 'neighbors' in the hydrodynamic kernel and number of SPH particles to which energy and momentum are coupled; using a single timestep or Strang-split integration scheme in the code; varying the Courant factor of the hydrodynamic solver; changing the order of operator-splitting for the cooling and feedback steps; or forcing equal vs. allowing separate gravitational softenings for baryons and dark matter. These produce very small ( < 10%) differences. We also varied the sizes of the high-resolution 'zoom-in' Lagrangian regions of the halos; the results here are insensitive to the region size if we choose sizes glyph[greaterorsimilar] 2 R vir (at the redshift of interest), but the  \n11 We have run a couple tests with 30 times higher particle numbers than our production-quality runs (for m12i and m10 ), to z = 2, and found that the stellar masses at this time and earlier vary by ∼ 10 -50% from those quoted here. However, this appears to be primarily stochastic, rather than systematic, so we suspect the masses will not change much further at still higher resolution.  \nFigure 6. M ∗ -M halo relation at z = 0, as Fig. 4. Top: Simulations with different numerical parameters: we show the effects of varied resolution, artificial viscosity, and the algorithmic implementation of feedback. We also compare a completely different version of SPH (with a different set of hydrodynamic equations), which is known to differ significantly in certain idealized hydrodynamic test problems. These have little effect on our predictions. Bottom: Effect of physical variation in stellar feedback properties. We compare runs with no stellar feedback, with no supernovae (but stellar winds, radiation pressure, and photo-ionization heating included), or with no radiative feedback (radiation pressure and local HII-heating). 'No feedback' runs generally predict M ∗ ∼ fb M halo, in severe conflict with the observations. 12 Removing radiative or SNe feedback also produce order-ofmagnitude too-large stellar masses. The non-linear combination of feedback mechanisms (not any one in isolation) is critical to drive winds and regulate galaxy masses.  \n<!-- image -->  \ncooling of halo gas and star formation are artificially suppressed if the high-resolution region is much smaller.  \nChanging the small-scale star formation prescription in the simulations has very little effect on our predictions. This is expected based on all of our previous studies using isolated (noncosmological) simulations (for explicit examples where we vary the density threshold and instantaneous 'efficiency' of star formation in dense gas by factors of > 1000, as well as the density, temperature, chemical, and virial-state dependence of star formation, see Paper I, Paper II, and Hopkins et al. 2013d). Globally, star formation is feedback-regulated : a certain number of young stars are required to balance gravitational collapse/dissipation, independent  \nof how those stars form. So long as cooling can proceed, they will form (what will change, via this self-regulation, if we change e.g. the density threshold above which stars form, is the amount of gas which 'builds up' above that threshold; see Hopkins et al. 2012b). In Fig. 6 we show examples where we vary the density threshold for star formation by factors of ∼ 100 (producing < 0 . 2 dex random/non-systematic changes in stellar mass), or impose a much stricter local virial criterion for star formation (local virial parameter < 0 . 5 instead of 1; producing a ∼ 20% difference in stellar mass); we have also experimented with removing the virial parameter entirely or dividing the instantaneous efficiency of star formation in dense gas by a factor of 100 (both produce < 10% changes).  \nWe have also investigated different purely algorithmic methods for coupling the same feedback physics. Subtle differences in the algorithmic implementation of feedback have little systematic effect on the stellar mass, provided the same mechanisms are included; however they can only be compared statistically, since the stochastic nature of feedback means that even very subtle changes can produce significant differences in the exact time history of bursts, for example. Of what we have considered, the most important parameter is how we implement the momentum gained during the Sedov-Taylor phase of SNe remnant expansion when the cooling radius is unresolved (see Appendix A). For example, one of the 'mod. SNe coupling algorithm' examples changes the particle weights (using a standard SPH kernel weight - effectively massweighted in the smoothing kernel - instead of a volumetric weighting) used to determine the coupling of SNe energy and momentum in the kernel. This can have dramatic effects on test problems: for a SNe in an infinitely thin, adiabatic disk with a low-density exterior, a mass-weighting couples all the momentum in the disk plane, instead of the vertical direction (the correct solution). Nevertheless we see this has relatively weak ( ∼ 20%) effects on the stellar mass and star formation history (in part because, in the average over many SNe over large volumes, all that matters is the total feedback input); however, it can significantly effect the morphological structure of e.g. the dense gas in a thin disk. We have also experimented with different functional forms for the ratio of the SNe energy and momentum coupling (producing small effects). The 'mod. RP algorithm' choice discretizes our radiation pressure term (which is usually a continuous force) into intentionally very large ( > 500kms -1 ) 'kicks' (this keeps the same total momentum flux, but makes each such particle 'kicked' unbound) - unsurprisingly this suppresses star formation further, but only by a factor of ∼ 3. In the 'mod. RP+SNe' choice we discretize the radiation pressure into smaller kicks ( = 5kms -1 ) and see this has little effect (as expected). In all cases the results lie within the (rather large) range allowed by observations.  \nIn a companion paper, Kereš et al. (2014) consider the detailed effects of substantial changes to each aspect of our numerical method described in Appendix B. Here, we simply show a few basic comparisons. Considerable attention has recently been paid to differences between the results of grid codes and older SPH methods (such as that in Springel & Hernquist 2002) for certain problems (especially sub-sonic fluid mixing instabilities; see Agertz et al. 2007; Kitsionas et al. 2009; Bauer & Springel 2012; Vogelsberger et al. 2012; Sijacki et al. 2012; Kereš et al. 2012). The numerical method used for our standard simulations has been specifically shown to resolve most of these discrepancies (giving results quite similar to grid codes in test problems); this is verified in Hopkins (2013) for standard test problems and Kereš et al. (2014) for cosmological simulations. However we have re-run some of our simulations using the Springel & Hernquist (2002) formulation of  \nSPH (described in Appendix B), which shows the most pronounced forms of these discrepancies. Despite the known differences between such methods for certain test problems, we find in Fig. 6 ( top panel) that it makes little difference for the predicted galaxy masses. The older SPH method gives slightly lower M ∗ ( M halo ) (by about ≈ 0 . 15dex), primarily because cooling of diffuse 'hot halo' gas is suppressed by less-efficient mixing. But for this specific question, the effect is quite small compared to the effects of including the appropriate stellar feedback physics. We also show an experiment where adopt an entirely distinct artificial viscosity prescription (see Appendix B for details), which produces negligible differences.  \nIt is important to stress that our conclusion here - that our results depend only weakly on the numerical details - applies to the galaxy stellar masses and other lowest-order, integrated quantities. In future work, we will study other properties of the simulations, such as the galaxy morphologies, which can (and do) depend on some parameters much more sensitively. For example, the modifications to the SNe coupling algorithm described above, which produce very little systematic change in our predicted stellar masses and star formation histories, produce surprisingly large changes to the angular momentum content and thickness of disks in the more massive galaxies.  \nFinally, we show these results in part to stress, emphatically, that while there are always numerical choices in any code, there has been no 'tuning' of these parameters for our study here. Certainly none of these has been 'fit' or 'adjusted' to match any observations, and all the choices above are held constant across our standard set of simulations, using values calibrated from simple test problems (e.g. Hopkins 2013).", '5.4 (Strong) Dependence on Feedback': "The lower panel in Fig. 6 shows the effect of varying the physics of feedback: now, we see dramatic differences in M ∗ ( M halo ) . Removing all feedback (every mechanism listed in § 3.3), gas cools and collapses on a dynamical time t dyn within the disk, forming stars at a rate ˙ M ∗ ∼ M gas / t dyn ∼ ˙ M gas where ˙ M gas is the inflow rate from the halo. Most of the baryons are turned into stars. 12  \nIf we turn off SNe feedback, but retain all other forms of feedback, the results are nearly as bad: again, M ∗ is severely overpredicted in both dwarfs and MW-mass systems. In the glyph[lessorsimilar] 10 10 halos, with no SNe, other forms of feedback may still suppress SF significantly (so M ∗ glyph[lessmuch] fb M halo), but the masses are still much too large relative to those observed by factors of ∼ 100. We also note that, as many previous studies have pointed out (Murray et al. 2005; McKee & Ostriker 2007; Shetty & Ostriker 2008; Faucher-Giguère et al. 2013), it is ultimately the momentum injected by SNe, not just the thermal energy, which regulates star formation. So as expected, if we artificially turn off the SNe momentum (coupling only thermal energy, as is common in many cosmological simulations), then in our simulations of massive ( > 10 12 M glyph[circledot] ) halos, this is nearly as bad as removing SNe entirely. In the lowest-mass dwarfs, the discrepancy is not so severe (factor glyph[lessorsimilar] 2 changes in the SFH), because the mass resolution ( ∼ 100 M glyph[circledot] ) is such that the early expansion phases of SNe remnants (in which the thermal energy begins to be converted into momentum) are well-resolved.  \nFigure 7. Comparison of the M ∗ ( M halo ) relation (as Fig. 4) predicted by other published simulations in the literature using sub-grid stellar feedback models. We compile these results (colored points), where available, at lowz ( z = 0 -0 . 5; top ) and highz ( z = 2 -3; bottom ). We compare against the simulations presented here (gray points) with explicit feedback, and observational constraints (lines). Even sub-grid models which are 'successful' near ∼ L ∗ at z ∼ 0 over-predict M ∗ ( M halo ) by an order-of-magnitude relative to our explicit feedback simulations and observations at both low masses ( M halo glyph[lessorsimilar] 10 10 M glyph[circledot] ) and/or high redshifts ( z glyph[greaterorsimilar] 2). The exceptions appear to be the newest generation of sub-grid models which have been explicitly adjusted to mimic the effects of radiative feedback as well as SNe, seen in our explicit feedback models: this includes Stinson et al. (2013); Aumer et al. (2013); Ceverino et al. (2013); Trujillo-Gomez et al. (2013).  \n<!-- image -->  \nIf we remove radiative feedback entirely (both radiation pressure and local photo-ionization and photo-electric heating, as described in § 3.3), but retain SNe (and stellar winds), we see a nearly identical failure (to the no-SNe case) in both dwarfs and massive galaxies: while M ∗ < fb M halo, far too many stars form. As we showed in Paper II, these mechanisms are critical to disrupt the dense regions of GMCs in which young stars are born, before SNe explode, and thus allowing the SNe to heat larger, lower-density volumes of gas (which can both avoid over-cooling and feel the collective effects of many SNe rather than just one), and therefore actually generate significant galactic outflows. The same result is found (on smaller scales) in much higher-resolution simulations of either single star clusters or the first stars, which directly treat the radiation-hydrodynamics with each single star as a source (e.g.  \nOffner et al. 2009; Krumholz et al. 2011; Tasker 2011; Wise et al. 2012).  \nInterestingly, in the dwarfs, if we turn off only radiation pressure, or only photo-ionization heating, the effect is much less severe: the predicted stellar mass is still significantly larger, but it is > 100 times larger when both are removed. Radiation pressure can, to some extent, 'make up for' the loss of photo-heating, and vice versa. This should actually not be surprising: the most massive GMCs in dwarf galaxies have local characteristic velocities < 10kms -1 , thus either HII heating or UV radiation pressure alone can disrupt them (though we expect, under these conditions, HII heating should dominate), and this is completely consistent with both observations of star-forming regions (e.g. Lopez et al. 2011) and numerical radiation-hydrodynamic simulations of low-density, low-velocity clouds (Harper-Clark & Murray 2011; Sales et al. 2013). And indeed this tradeoff between photo-heating and radiation pressure in small clouds is exactly what we saw in our ultrahigh resolution simulations of isolated dwarfs of the same mass (discussed extensively in Paper II; see Figs. 7, 9, and 14-19 therein).  \nIn the massive systems, on the other hand, the radiation pressure term becomes more important than the HII heating. We see this in tests with both m12q and m12v . Even when the difference in stellar mass is not large (e.g. the m12v case), the lack of radiation pressure feedback is particularly evident in the dense, earlyforming center of the galaxy, where in the runs without radiation pressure feedback an enormous central density 'spike' appears, leading to a very large circular velocity of ∼ 400kms -1 in the central regions of these systems. At these densities, HII photo-heating is dynamically insignificant.  \nIf we disable stellar wind feedback (specifically, retaining stellar winds as a source of mass and metals, but associating no energy or momentum with that mass), and retain all other feedback, we see relatively weak effects. This is not surprising: their momentum flux is comparable to but not larger than other sources, and their energetics are much less than SNe. But they are obviously an extremely important source of mass and metals in the ISM.", '5.5 Comparison to Previous Work': 'In Fig. 7, we compare our results (grey points) at low and high redshifts, to those from previous simulations spanning a wide range of galaxy properties and numerical methods (Pelupessy et al. 2004; Stinson et al. 2007, 2013; Mashchenko et al. 2008; Valcke et al. 2008; Governato et al. 2010, 2012; Oser et al. 2010; Feldmann et al. 2009; Brooks et al. 2011; Guedes et al. 2011; Sawala et al. 2011; Scannapieco et al. 2011; De Rossi et al. 2013; Okamoto 2013; Kannan et al. 2013). All of these simulations include some form of sub-grid model designed to mimic the ultimate effects of stellar feedback, although the prescriptions adopted differ substantially between each. Most of these models are specifically tuned to reproduce reasonable scaling for MW-mass systems at z ∼ 0. However, two discrepancies are immediately evident. First, nearly all the previous models predict much larger stellar masses in dwarf galaxies with M halo glyph[lessorsimilar] 10 11 . 5 M glyph[circledot] , compared to either our simulations or the observational constraints. Second, even simulations which produce excellent agreement with the observations at z = 0 tend to predict far too much star formation at high redshift (take e.g. the simulation in Guedes et al. 2011, which produces a MW-like system with many properties consistent with observations at z = 0, but has turned nearly all its baryons into stars at z glyph[greaterorsimilar] 2).  \nThese are similar to the discrepancies that appear when we re-run our simulations excluding radiative feedback. And indeed, nearly all of the models from the literature in Fig. 7, even given  \nvarious freely adjustable parameters, are designed and motivated only to reproduce the effects of supernova feedback, which we have shown is insufficient to explain the observations.  \nIn fact, the only sub-grid models, to our knowledge, which currently do not produce such discrepancies (and agree broadly with our simulations both at low masses and high redshifts) are the recent generation of models in Stinson et al. (2013); Aumer et al. (2013); Ceverino et al. (2013); Trujillo-Gomez et al. (2013) (for some additional results from these see Kannan et al. 2013). These new models (all of which have been developed recently) are specifically designed/tuned to mimic the effects of radiative feedback (albeit indirectly), and to reproduce via simple sub-resolution prescriptions (including turning off cooling) some of the most important effects of radiation pressure and photo-heating which were studied in our previous work (Hopkins et al. 2012d). 13 Whether this is unique or not remains to be tested; the phase structure and other properties of outflows and the CGM in such models can be very different from those predicted here, even for the same massloading efficiencies (discussed further below). It will be particularly interesting to see whether other recently-developed sub-grid models such as that in Agertz et al. (2013), also incorporating the effects of radiative feedback but via very different prescriptions, will also agree well with observations at both low and high redshifts. In any case, these comparisons - and the results from this new generation of sub-grid models - highlight that some accounting for non-SNe feedback is critical.', '5.6 Instantaneous Suppression of Star Formation (at Fixed Gas Densities)': "We now examine galaxy star formation rates. In the previous section, we showed that the integrated SF is suppressed with feedback. But equally important is that feedback suppresses instantaneous SFRs in galaxies. This is manifest in the Kennicutt-Schmidt (KS) relation, shown in Fig. 8. 14 We plot the simulations at all redshifts (the redshift evolution is insignificant), and compare to observations at a range of redshifts (which also find little or no evolution). 15  \nThe predicted KS law agrees well with observations at all redshifts. As shown in Paper I-Paper III, this emerges naturally as a consequence of feedback, and is not put in by hand. Recall that the instantaneous SF efficiency (SF per dynamical time) in dense gas in the simulations is 100%; however the global SF efficiency is ∼ 2%. This difference arises because at ∼ 2% efficiency, feedback injects sufficient momentum to offset dissipation (indeed, given the same feedback, we obtain the identical KS law independent of the details of our small-scale SF law; see Hopkins et al. 2011, 2013d).  \n14 We define Σ SFR = ˙ M ∗ /π R 2 SFR (where ˙ M ∗ is the total SFR and R SFR is the half-SFR radius) and Σ gas = M gas /π R 2 SFR (where M gas is the gas mass within the 90% SFR radius). Defining both ˙ M ∗ and M gas within R SFR or the stellar effective radius shifts the points along the relation.  \n15 We compile the observed local galaxies in Kennicutt (1998) and Bigiel et al. (2008), and high-redshift galaxies in Genzel et al. (2010) and Daddi et al. (2010); shaded region shows the 90% inclusion range at each Σ gas from the compilation. As discussed in those papers, there is no significant offset between the high and low-redshift systems at fixed Σ gas.  \n4  \n<!-- image -->  \n/circledot  \nFigure 8. Kennicutt-Schmidt law, observed (Kennicutt 1998; Bigiel et al. 2008; Genzel et al. 2010; Daddi et al. 2010, yellow shaded range) and simulated (points as Fig. 4). We emphasize that this is a prediction: the instantaneous SF efficiency per dynamical time in dense gas is 100% in the simulations, but the emergent KS-law, as a consequence of feedback, has an efficiency a factor ∼ 50 lower. As shown in Paper I and Hopkins et al. (2013d), this is insensitive, with resolved feedback models, to the smallscale star formation law, and entirely determined by stellar feedback. 'No feedback' models lie a factor ∼ 50 above the observations; 'no radiation' and 'no SNe' models (Fig. 6) lie a factor ∼ 10 above observations.Figure 9. Example star formation history (SFH) for the m12v simulation, in our standard (explicit feedback) model compared to different sub-grid feedback treatments. We show the formation history of all stars in the simulated box at z = 0 (smoothed in 10 7 yr bins); this is not qualitatively different from the SFR versus time of the largest 'main' galaxy in the box at each time. 'No feedback' models force the galaxy to lie on the KS-law (SF is 'slow') but do not expel gas; gas piles up until the SFR balances the halo accretion rate, with a broad peak from z ∼ 2 -6. 'Sub-grid wind' models 'kick' gas at a rate proportional to the SFR; we show examples with different efficiencies and implementations. By design, model '2' produces a nearly identical z = 0 stellar mass M ∗ to our explicit feedback model. These sub-grid models (by construction) lower the SFR, but in both cases leave the qualitative behavior of the SFH identical. Explicit feedback models not only suppress the total M ∗ formed, but change the shape of the SFH. SF is more 'bursty' on small timescales, and the SFR is flatter in time (more biased to late times, without the broad highz peak).  \n<!-- image -->  \nFigure 10. SFH for each 'main' (largest) z = 0 galaxy in our standard (explicit-feedback) simulations. Lines show the mean SFR averaged on timescales of 10 8 yr (black solid) and 10 9 yr (red dashed). With explicit feedback, SFRs are highly variable below the galaxy dynamical time. Moreover, the (averaged) SFRs tend to be flat and/or rising with time. In contrast, with no feedback, the SFH has a sharp rise and fall peaking at z ∼ 2 -6. In the least massive dwarfs ( m09 ; M ∗ < 10 6 M glyph[circledot] and V vir ( z = 0 ) < 20kms -1 ), the SFR is strongly suppressed after reionization once a combination of the ionizing background and some small amount of feedback from the early star formation is able to expel most of the halo gas and prevent new cooling. We compare our m11 , m12 , and m13 runs to the observationally inferred 'mean' tracks (colored lines) for the main galaxies in halos of the same z = 0 mass, from Moster et al. (2013, magenta) and Behroozi et al. (2012, cyan). In each case the lines bracket the 1 σ range/scatter in the observed galaxy population. Our m11 and m12 runs agree very well with these constraints; however, in the most massive systems ( m13 ), the galaxy never 'quenches,' and the SFR remains high in conflict with observations below z ∼ 1.  \n<!-- image -->  \nFigure 11. Top: SFR versus galaxy stellar mass at different redshifts. We compare the observed (best-fit) relations from the compilation in Behroozi et al. (2012, black dashed) and Zahid et al. (2012, blue dashed) (the systematic offset is typical of different calibrations). Allowing for the typical factor ∼ 2 systematic observational calibration uncertainty, the agreement is good at all z . However, magenta + 's compare low-resolution (100 pc softening) runs of some massive halos which produce too-massive galaxies at z = 0: there is little offset between these simulations and our fiducial models. The observed relation is simply a consequences of galaxies having relatively 'flat' star formation histories. Bottom: Specific SFR of galaxies with different M ∗ , versus redshift. Observations are compiled in Behroozi et al. (2012, Table 5) and Torrey et al. (2013). The dynamic range here is smaller so the plot appears noisier, but the information is identical to that at top . SSFRs at z glyph[greaterorsimilar] 2 are relatively flat, indicating rising SF histories at highz .  \n<!-- image -->  \nIf we instead consider simulations with weak/no feedback, the global KS relation is severely over-predicted (efficient cooling leads to global efficiencies ∼ 100%). In most cosmological simulations, this is offset 'by hand' by simply enforcing a large-scale SF law that is sufficiently 'slow' that it agrees with the observations; however we see that this is already (implicitly) a sub-grid feedback model. Including explicit feedback obviates the need for these prescriptions, meaning that the instantaneous SF properties are truly predictive, and not simply a consequence of our chosen small-scale SF law.", '5.7 Global Star Formation Histories': "In Fig. 9, we examine the SFH of one MW-mass galaxy. We compare this to common sub-grid models. First, a 'no-feedback' model following Springel & Hernquist (2003a); this includes only a subgrid model for the effects of stellar feedback on the ISM structure (an 'effective equation of state') which ensures, by design (via tuned parameters) that the galaxy lies on the Kennicutt-Schmidt relation and has reasonable gas densities. However, without galactic winds, gas from inflows quickly builds up and the SFR rises until ˙ M ∗ ≈ ˙ M inflow, and nearly all the baryons are turned into stars. The galaxy at z = 0 is far too massive, and most of its stars are old (formed at z glyph[greaterorsimilar] 2, with the SFR peaking at z ∼ 5). 16 We then add a sub-grid wind model, in which gas is 'kicked' out of the galaxy (forced to free-stream to ensure it escapes the disk) at a rate proportional to the SFR: here the mass-loading is equal to the SFR ('subgrid wind 1'). This suppresses the SFR (as it is intended to do), by about a uniform factor ∼ 2, as expected. However this still leaves a too-massive galaxy, with most of its stars formed very early. Next, we consider a stronger wind model ('sub-grid wind 2'): the massloading is doubled, with the free-steaming length fixed. This further suppresses the SFR - in this model the final stellar mass agrees reasonably well with our explicit-feedback simulation. However, the sub-grid model still produces a SFR which peaks at very high redshifts z ∼ 2 -6. The problem is that in all the sub-grid models regardless of the absolute suppression of the integrated SFR or position on the Kennicutt-Schmidt relation - the shape of the galaxy SFH still closely resembles the shape of the halo inflow rate vs. time (for examples of this with other sub-grid models, see Oppenheimer & Davé 2008; Scannapieco et al. 2011; Stinson et al. 2013; Puchwein & Springel 2013).  \nThese broadly peaked SF histories are disfavored by a variety of observations. They produce too-massive galaxies at high redshift, as discussed above. But they also produce galaxies with SF histories at highz that disagree with direct observational constraints (see Papovich et al. 2005; Reddy et al. 2006; Stark et al. 2009).  \nWith our full, explicit feedback model included, we see that the shape of the SF history is qualitatively changed, and is more consistent with observations. At all times, SFRs are much more time-variable (this is discussed below). At the highest z glyph[greaterorsimilar] 6, halo and stellar masses both grow efficiently (albeit with some offset). 17  \n16 As noted in the many previous simulations with this sub-grid model (e.g. Springel & Hernquist 2003b), the 'effective equation of state' approach does not allow cooling below ∼ 10 4 K, so if those simulations properly included molecular or fine-structure cooling (with the appropriate high resolution), the SFH might peak even earlier.  \nThis is the 'rapid assembly' phase, before/during reionization, in which feedback - while able to eject some gas from the galaxy and provide some overall suppression and variability of ˙ M ∗ - does not appear to dominate the gas dynamics (the central potential and mass of the halo grow on timescales comparable to the galaxy dynamical time; so ˙ M ∗ ∝ ˙ M halo). But from z ∼ 2 -6, feedback acts strongly, and there appears to be a maximum, steady-state SFR which is constant or slowly increasing with time at which the galaxy is able to cycle new material into a fountain and so maintain equilibrium. This 'quasi-equilibrium' SFR scales with the central potential of the galaxy (see Paper III), as traced by quantities such as the central halo density or V max (the maximum circular velocity), not the halo mass or virial velocity. The central potential depth increases only weakly over this time as halos accrete material on their outskirts. Below z ∼ 2, a competition ensues between slowing halo accretion rates and more highly-enriched halo gas raising cooling rates. Individual mergers also have a more dramatic effect on SF histories.  \nIn Fig. 10, we show the SFH for each main z = 0 galaxy in our simulations, 18 and see that all cases with 10 9 glyph[lessorsimilar] M halo glyph[lessorsimilar] 10 13 M glyph[circledot] exhibit similar (relatively flat or slowly rising) SFHs. 19 In the most massive halos, some decline occurs when M halo glyph[greaterorsimilar] 10 12 M glyph[circledot] , as the cooling time of virialized gas becomes longer relative to the dynamical time (the system transitions to 'hot mode' accretion and filamentary infall is suppressed). However, we stress that the galaxies are clearly not 'quenched' - every system we simulate is still very much a star-forming, blue galaxy at z ∼ 0 (our m13 simulation would need a SFR glyph[lessmuch] 1 M glyph[circledot] yr -1 at z = 0 to be 'red and dead' by most definitions, but its SFR is ∼ 5 M glyph[circledot] yr -1 ). In very low mass halos (e.g. our m09 , with V vir ( z = 0 ) < 20kms -1 ), cooling is strongly suppressed after reionization.", "5.8 Specific SFRs and the SF 'Main Sequence'": "Fig. 11 compares the galaxy-integrated SFRs in all our simulated systems (including non-main halos) with observations of the SFR or specific SFR (SFR / M ∗ ) as a function of galaxy stellar mass, at various redshifts. The simulations agree well with the SFR 'main sequence' (SFR -M ∗ relation) observed at all z (observations plotted include compilations from Erb et al. 2006; Noeske et al. 2007; Daddi et al. 2007; Elbaz et al. 2007; Stark et al. 2009 and others in Behroozi et al. 2012 (see Table 5 therein) and Zahid et al. 2012. The scatter is also similar to that observed. There may be some slight tension (the predictions being slightly high at z = 0 and low at z = 2), but these are well within the range of systematic uncertainties owing to different SFR calibrations (we show a couple such examples). By extension, the simulations similarly agree with the evolution in specific SFRs of galaxies as a function of mass.  \nEvolution in specific SFRs and SFR -M ∗ towards higher  \nthe details of how feedback is coupled, even though the later-time SFRs and final stellar masses are robust to these variations. See Appendix A for details.  \n19 Interestingly, the m12q simulation shows a much higher highz SFRthan m12v or m12i . This is in part because the particular choice of a 'quiescent' halo led, in this case, to a halo with relatively little growth at late times ( z glyph[lessorsimilar] 3), hence a particularly early 'formation time.'  \nFigure 12. Variability of the SFHs shown in Fig. 10, quantified versus timescales ∆ t avg. For each 'main' galaxy in each simulation, we show the logarithmic dispersion in the SFR σ SFR about its mean on longer timescales, when the SFR is time-averaged over the timescale ∆ t avg. The variability rises substantially on timescales ∼ 10 7 -10 8 yr (galaxy dynamical times), owing to a combination of fountain dynamics, local structure in the galaxies, and stochastic effects from individual star forming regions. The shorttimescale variability is a factor of ∼ 2 -3 in ∼ L ∗ galaxies, but rises to order-of-magnitude level in dwarfs (where individual star clusters and bursts have a more dramatic effect).  \n<!-- image -->  \nSSFR at highz simply reflects rising gas fractions (as it must, since the simulations lie on the same KS-law in Fig. 8). The 'flattening' of SSFR at highz implies SF histories of individual galaxies are rising with time (as we see directly); physically it follows from the saturation of gas fractions at large values, and rapid growth of halo mass at these times. The SFR -M ∗ relation is, to lowest order, just ˙ M ∗ ∼ M ∗ / t Hubble ( z ) - this must be trivially true in any scenario where SFRs are relatively flat and/or rising with time (typical of star-forming galaxies). For this reason we see the same relation even in our simulations without feedback, as have other simulations with different feedback prescriptions (see Kereš et al. 2009; Davé et al. 2011). And we see that even the very massive halos (which produce 'too large' an M ∗ at low redshifts) lie on the extension of the observed relation (the problem is that they continue on the relation, rather than 'quenching' and moving below it, as observed at high masses).", '5.9 Quantifying Burstiness/Variability in SFRs': "In Fig. 9, we showed that the SFRs are significantly more timevariable in models with explicit/resolved feedback as compared to sub-grid feedback models. We quantify this in Fig. 12. We measure the dispersion in the SFR smoothed over various time intervals. Unsurprisingly, the scatter is larger on small timescales. On glyph[greatermuch] 10 8 yr timescales, the variability is always small (SFHs are 'smooth') this is more a function of the evolution of the halo over a Hubble time. Some such long-timescale variability is driven by mergers and global gravitational instabilities, but much of the short-timescale variability is not connected to these phenomena. Rather, on smaller timescales (comparable to the galaxy dynamical time) the dynamics of fountains, feedback, and individual giant molecular clouds and star clusters becomes important, so the scatter increases down to timescales ∼ 10 6 yr (comparable to the massive stellar evolu-  \non timescale). 20 The short-timescale scatter is modest ( ∼ 0 . 3dex) for massive systems ( M halo glyph[greaterorsimilar] 10 12 M glyph[circledot] ), but rises in smaller halos (where even single star clusters can have large effects) to ∼ 1dex at M halo glyph[lessorsimilar] 10 10 M glyph[circledot] . 21", '6.1 Key Results and Predictions': "Wepresent a series of cosmological zoom-in simulations 1 of galaxies with M halo ∼ 10 9 -10 13 M glyph[circledot] and M ∗ ∼ 10 4 -10 11 M glyph[circledot] . At this time, several of these runs represent the highest-resolution in both mass and force resolution of any fully cosmological runs to z = 0. But the most important improvement, compared to previous simulations, is that we for the first time include a fully explicit treatment of both the multi-phase (cold molecular through atomic, ionized, and hot diffuse) ISM and stellar feedback. Our treatment of the ISM is enabled both by our resolution and improved treatment of cooling and heating physics (e.g. molecular and metal line cooling, photoionization and photo-electric heating with self-shielding). Our stellar feedback model utilizes explicit time-dependent energy, momentum, mass and metal fluxes taken directly from stellar population models, without free/adjustable parameters. As such, the SFRs in our simulations, the resulting outflows, and galaxy stellar masses are not the result of tuning or 'by hand' adjusting feedback efficiencies. Critically, we include not just thermal energy from SNe, but the momentum and energy associated with SNe Types Ia & II, stellar winds (young star & AGB), local photo-ionization and photo-electric heating, and radiation pressure from UV and IR photons. In addition, our formulation of SPH resolves the historical numerical problems with this method, especially important for cooling in hot halo atmospheres (see Appendix B; Hopkins 2013; Keres et al. in prep).", 'Our key conclusions include:': "- · Stellar feedback - from known sources including SNe (energy and momentum), stellar winds, radiation pressure (primarily optical/UV), and photo-heating - is both necessary and sufficient to explain the observed relation between galaxy stellar mass and halo mass, and by extension the shape of the galaxy mass function and clustering, at stellar masses M ∗ glyph[lessorsimilar] 10 11 M glyph[circledot] . This appears to be true at all redshifts.\n- · No one feedback mechanism alone is sufficient: the effects add non-linearly, and the common approximation in simulations of including only SNe feedback severely over-predicts galaxy masses (especially at low masses and/or high redshifts). The effects are even worse if the feedback momentum is ignored (if only thermal energy is considered).\n- · The M ∗ -M halo relation evolves very weakly with redshift (because outflow efficiencies depend mostly on the central binding energy within the galaxy). At z glyph[greaterorsimilar] 2, weak evolution towards higher M ∗ ( M halo ) at low masses is equivalent to a steepening faint-end  \nslope of the galaxy luminosity function, similar to what is inferred observationally (Bouwens et al. 2007; Stark et al. 2009).  \n- · Stellar feedback and standard cooling physics explain low galaxy stellar masses, but do not appear sufficient to explain 'quenching' (late time suppression of star formation in massive halos) - none of our massive systems are 'red and dead.'\n- · Our simulations reproduce the observed Kennicutt-Schmidt relation. This is despite the fact that we assume a small-scale SF efficiency of 100% in self-gravitating dense gas. As such, the KS law and instantaneous SFRs are truly predicted, not simply a consequence of our sub-grid SF law. The low star formation efficiency we find is a consequence of stellar feedback, not the microphysics of how stars form in dense gas. Absent feedback, efficient cooling leads to a global SF efficiency of ∼ 100% per dynamical time. With feedback - from the same mechanisms that produce largescale outflows and regulate galaxy formation - the SF efficiency self-regulates at ∼ 2%, the level where feedback injects sufficient momentum to offset dissipation.\n- · Realistic feedback changes the shape of galaxy star formation histories. In particular, feedback from stellar radiation (both photoheating and radiation pressure) is critical for disrupting dense, cold gas, and so is especially important for suppressing star formation in high redshift galaxies. This leads to much flatter, or gently rising, star formation histories in subL ∗ galaxies. Most previous sub-grid models give qualitatively different results, in conflict with observations.\n- · The observed star formation 'main sequence' and specific SFRs emerge naturally from the shape of the galaxies' star formation histories (from M ∗ ∼ 10 8 -10 11 and z ∼ 0 -6). This includes 'flat' SSFR evolution at z ∼ 2 -6. However these are relatively insensitive to feedback, since any broadly flat or rising SF history predicts M ∗ ( z ) ∼ t Hubble ( z ) 〈 ˙ M ∗ ( z ) 〉 , consistent with the observations.\n- · Dwarf galaxies exhibit much more 'bursty' SF histories, with large variability in their SFRs on short timescales ( ∼ 1dex scatter on glyph[lessorsimilar] 10 7 yr timescales). This is because star formation and star cluster formation, and their associated feedback, are stochastic. The variability is not driven by mergers or global gravitational instabilities. Massive ( ∼ L ∗ ) galaxies are much less variable ( ∼ 0 . 3dex scatter in SFRs). This may translate into significantly larger scatter in M ∗ ( M halo ) at dwarf masses compared to ∼ L ∗ galaxies.", '6.2 Numerical Methods': "Wesee relatively weak dependence on simulation resolution, which is perhaps surprising given the small-scale structure present in the ISM. However, in Papers I-III & Appendix C, we presented extensive resolution studies of isolated disk galaxies simulated using the same prescriptions but numerical resolution varied from values comparable to those here, to order-of-magnitude superior mass and spatial resolution. We showed that the galaxy-averaged SFR is one of the very first quantities to converge, and is consistent to within factor ∼ 2 even for relatively poor resolution: this is because it traces the integral effect of feedback balancing turbulent dissipation. That said, we do see qualitative changes in behavior if the resolution falls below that needed to meaningfully resolve ISM phase structure, at which point 'self-regulation' by feedback loses meaning. As a rule of thumb, the simulations must at least resolve the Toomre/Jeans length and mass (the size of the largest GMCs) in each star-forming disk, and the results are especially numerically stable if the mass resolution can be ideally pushed to < 10 4 M glyph[circledot] . Even in this regime, quantities such as the phase structure of dense  \ngas and outflows are much more sensitive to resolution, and will be discussed in more detail in future work.  \nGiven the same feedback model, we also see little difference in the stellar mass buildup between our standard simulations, run with a numerical algorithm designed and shown to eliminate essentially all major differences between grid (Eulerian) and smoothedparticle (Lagrangian) hydrodynamics methods, and an older version of SPH that exhibits large differences in test problems. Thus we expect little or no difference between the results here and those from adaptive-grid or moving-mesh codes, if the same feedback and ISM physics could be included. This owes to two key points: first, the differences between numerical methods in this respect, even where significant, are generally much smaller than the ordersof-magnitude differences owing to the inclusion or exclusion of the relevant physics. The stellar mass content of galaxies is set by the total amount of feedback injected, and so it is unsurprising that the time-averaged star formation rate is insensitive to changes in the detailed phase structure of the gas around galaxies. Second, the numerical differences primarily affect mixing instabilities in multiphase, sub-sonic, pressure-dominated gas. As such, many comparison studies have shown that while the numerical differences can be important for details of the structure of hot halos in massive galaxies, they are generally unimportant inside cold star-forming gas, or in subL ∗ halos, where the flows of interest tend to be highly supersonic and gravity-dominated (see e.g. Kitsionas et al. 2009; Price &Federrath 2010; Bauer & Springel 2012; Sijacki et al. 2012). As one considers more detailed galactic properties, we expect the differences between numerical methods to manifest as discrepancies in the cooling properties, phase structure, or distribution of heavy elements in the CGM, and to impact the way in which both inflowing cool gas and feedback-driven outflows interact with gas in galactic halos. For these reasons, an accurate numerical scheme is critical if one hopes to study the detailed structure of both gas in and around galaxies with realistic feedback. A much more extensive comparison of numerical methods is presented in a companion paper (Kereš et al. 2014).", '6.3 Future Work': "This is a first exploration of cosmological simulations with explicit stellar feedback models, and many open questions remain. We have studied the effects of realistic stellar feedback on galaxy star formation histories and stellar masses; however, a complete understanding of this self-regulation requires a much more detailed examination of the dynamics of galactic outflows. In companion papers, we will study how outflows are generated, and how these interact with the circum-galactic and inter-galactic medium. It will be particularly important to build new observational diagnostics and explore whether or not different feedback mechanisms lead to different observable properties in the ISM, CGM, and IGM (Faucher-Giguere et al. 2014). Complementary questions regarding the morphology of galaxies - how the sizes, bulge-to-disk ratios, kinematics, and other properties of the simulated systems here depend on different feedback mechanisms - will be developed as well. The resolution and explicit treatment of the ISM in these simulations make possible many additional studies.  \nGoing forward, it will also be important to examine the role of additional physics. Some other physics is probably needed to explain the 'quenching' of star formation in massive systems ( M halo glyph[greatermuch] 10 12 M glyph[circledot] ). AGN feedback is a plausible candidate, which we have studied in previous work using idealized sub-grid models for the ISM. But the consequences could easily be completely different in a resolved multi-phase medium. Other physics such as  \nmagnetic fields, anisotropic conduction, and cosmic rays may be important as well, and their consequences are just beginning to be explored (e.g. Jubelgas et al. 2008; Hanasz et al. 2013; Salem & Bryan 2013).", 'ACKNOWLEDGMENTS': 'We thank the many friends and peers who discussed this work in progress and sent helpful suggestions after the first draft was posted to the arXiv. The simulations here used computational resources granted by the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number OCI-1053575; specifically allocations TG-AST120025 (PI Keres), TG-AST130039 (PI Hopkins), and TG-TG-AST090039 (PI Quataert). Collaboration between institutions for this work was partially supported by a workshop grant from UC-HiPACC. Partial support for PFH was provided by NASA through Einstein Postdoctoral Fellowship Award Number PF1-120083 issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the NASA under contract NAS8-03060, and by the Gordon and Betty Moore Foundation through Grant #776 to the Caltech Moore Center for Theoretical Cosmology and Physics. JO also thanks the financial support of the Fulbright/MICINN Program and NASA Grant NNX09AG01G. DK acknowledges support from the Hellman Fellowship fund at the UC San Diego and NASA ATP grant NNX11AI97G. CAFG is supported by a fellowship from the Miller Institute for Basic Research in Science and by NASA through Einstein Postdoctoral Fellowship Award number PF3-140106 and grant number 10-ATP10-0187. EQ is supported in part by NASA ATP Grant 12-ATP12-0183, a Simons Investigator award from the Simons Foundation, the David and Lucile Packard Foundation, and the Thomas Alison Schneider Chair in Physics.', 'REFERENCES': 'Zahid, H. J., Dima, G. I., Kewley, L. J., Erb, D. K., & Davé, R. 2012, ApJ, 757, 54  \nZuckerman, B., & Evans, II, N. J. 1974, ApJL, 192, L149', 'A1 Cooling': 'Gas cooling is solved implicitly each timestep (using the iterative algorithm from GADGET-3). Heating/cooling rates are computed including free-free, photo-ionization/recombination, Compton, photo-electric, metal-line, molecular, and fine-structure processes. We follow 11 separately-tracked species (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe, each with its own yield tables associated directly with the different mass return mechanisms below; see Wiersma et al. 2009b). The appropriate ionization states and cooling rates are tabulated from a compilation of CLOUDY runs, including the effect of a uniform but redshift-dependent photoionizing background computed in Faucher-Giguère et al. (2009), together with local sources of photo-ionizing and photo-electric heating (described below with the relevant feedback mechanisms). Self-shielding is accounted for with a local Sobolev/Jeans-length approximation (integrating the local density at a given particle out to a Jeans length to determine a surface density Σ , then attenuating the background seen at that point by exp ( -κ ν Σ) ). Confirmation of the accuracy of this approximation in radiative transfer experiments can be found in Faucher-Giguère et al. (2010) and Rahmati et al. (2013). With this accounting, metal-line cooling follows the rate tables from Wiersma et al. (2009a), free-free rates follow Katz et al. (1996), and photo-electric rates follow Wolfire et al. (1995). Compton heating/cooling is included both from the CMB and local sources, accounting as in Faucher-Giguere & Quataert (2012) for possible two-temperature plasma effects at very high temperatures by limiting the Compton rates by the Coulomb energy exchange rates (though in practice this is only relevant at much higher temperatures than are seen in these simulations). Fine-structure and molecular cooling at low temperatures ( T < 10 5 K) is tracked using an interpolation table for a compilation of CLOUDY runs as a function of the density, temperature, metallicity, and local ionizing background, as in Robertson & Kravtsov (2008). A temperature floor is included at the maximum of either 10K or the CMB temperature at the given redshift.', 'A2 Star Formation': 'Star formation occurs probabilistically. At each timestep d t , a gas particle has a probability of turning into a star particle p = 1 -exp ( -˙ m i ∗ d t / m i gas ) , where ˙ m i ∗ is the SFR integrated over the particle, and m i gas is the particle gas mass. The SFR is non-zero only for particles with density above n > n crit (generally n crit = 100cm -3 ), which are also locally self-gravitating using the criteria developed in Hopkins et al. (2013d) ( α ≡ δ v 2 δ r / Gm gas ( < δ r ) ≈ β ( |∇ × v | 2 + |∇· v | 2 ) / G ρ < 1, with β ≈ 0 . 25), and which have a non-zero molecular fraction f mol > 0. The molecular fraction is determined following Krumholz & Gnedin (2011), using the local Sobolev approximation and metallicity to estimate the integrated column to dissociating radiation ( τ ≈ κ 〈 Σ 〉 with 〈 Σ 〉 = ρ [ h sml +( ∇ ln ρ ) -1 ] and κ = κ gas + κ dust , MW ( Z / Z glyph[circledot] ) ). The SFR per unit volume for gas that meets all of these criteria is then 100% per free-fall time, ˙ ρ ∗ = ρ mol / t ff = f mol ρ gas / t ff. When a gas particle becomes a star particle, the star particle inherits the metallicity of each followed  \nspecies from its parent, and the conversion/formation time of the particle is used to determine its age in subsequent timesteps. The star particles also inherit their mass and gravitational softening from their parent gas particles, so that they are not over or underresolved relative to the medium from which they form. 22', 'A3 Stellar Feedback': "The stellar feedback algorithms follow those developed in Paper IPaper II, with some modifications necessary for cosmological simulations.  \nRadiation Pressure: Gas surrounding stars (see below) receives a direct momentum flux ˙ P rad ≈ ( 1 -exp ( -τ UV / optical ))( 1 + τ IR ) L incident / c where 1 + τ IR = 1 +Σ gas κ IR accounts for the absorption of the initial UV/optical flux and multiple scatterings of the re-emitted IR flux if the region between star and gas particle is optically thick in the IR (assuming the opacities scale linearly with gas metallicity). At each timestep we evaluate the optical depth in a smoothing kernel around the star particle (whose pre-absorption stellar spectrum L ν is tabulated as a function of age and metallicity). The UV/optical absorption ( τ UV / optical = ∫ κρ d glyph[lscript] ) is estimated via the Sobolev approximation (as above) in multiple frequency bins (see Paper II for details). The absorbed fraction of L is then distributed within the SPH smoothing kernel according to each particles' relative contribution to the optical depth. This absorbed luminosity is assumed to re-radiate isotropically in the IR (and while it can be re-absorbed, it is again re-radiated in the IR), so for an (assumed) gray-body opacity it imparts an acceleration a IR = κ IR F IR / c to all gas in the kernel. 23  \nPhotons which are not absorbed in the UV/optical, and the reemitted IR flux, define an effective 'emergent spectrum' for each kernel. This is propagated to large distances in the gravity tree, where it is used to calculate the local incident flux on all gas particles from stars outside the smoothing kernel; the same frequencydependent opacities are used to calculate the local absorption and momentum flux ( ˙ P = L abs / c ). The momentum flux is imparted in every timestep along the direction determined by the flux-limited diffusion approximation (along the local optical depth gradient). 24  \n22 As discussed in the main text, we have confirmed (as seen in previous simulations of isolated galaxies in Hopkins et al. 2011 and Hopkins et al. 2013d) that with feedback active, the star formation prescription makes little difference to our results. We have re-run our m09 , m10 , and m12v runs removing the virial criterion and/or molecular criterion from the star formation law, changing the SF density threshold from ∼ 10 -1000cm -3 , and changing the SFR per free-fall time from ∼ 10% -200%. These changes yield only small (factor < 2), random (non-systematic) changes to the star formation history and resulting stellar mass.  \n- 23 If the gas is optically thick in the IR out to the edge of the smoothing kernel, the kernel is iteratively expanded so this region is treated explicitly. But this is almost never the case at our resolution.  \n24 Note that this avoids the 'clump detection' algorithm described in Paper I. That was important there, where regions which were optically thick in the IR (e.g. cores within GMCs) could be well-resolved, so coherent radiative transfer effects required a means to estimate the clump 'membership' and so avoid artificial numerical cancellation of momentum (which would occur if momentum fluxes were randomly oriented within the clump, as shown in Fig. A1 therein). The importance is decreased here because such regions are always at most marginally resolved. But in either case, computing the optical depth gradient allows us to recover a very close approximation to the previous information without the added neighbor search (which is very expensive in cosmological simulations), and is more faithful to full radiative transfer in the diffusion limit. In Appendix C we show both algorithms give very similar results in the isolated galaxy simulations from Paper I. We have also re-run our m12v simulation with the older 'clump detection' algorithm  \nPhoto-Ionization and Photo-Electric Heating: Here the algorithm is identical to that in Paper II. We first tabulate the rate of production of ionizing photons for each star particle (as a function of age and metallicity); moving radially outwards from the star, we then ionize each neutral gas particle until the photon budget is exhausted (using the gas density, metallicity, and ionization state to determine the necessary photon number). Note that - unlike the purely local Stromgren sphere approximation sometimes used in the literature - this accounts for whether each particle is already ionized, and (if so) allows the ionized region to continue to expand (thus accounting for large coherent/overlapping HII regions). 25 In the cooling routine, ionized gas is flagged as having a sufficiently strong local ionization field to keep it fully ionized for the duration of the timestep; this local ionizing field information together with the escaped UV flux defined above is also used to determine the photo-electric heating rates in the cooling routine.  \nSupernovae and Stellar Winds: The SNe Type-I and TypeII rates are tabulated from Mannucci et al. (2006) and STARBURST99, respectively, as a function of age and metallicity for all star particles; this determines a probability per unit time d p = ( d N SNe / d M ∗ d t ) m ∗ , i d t which we use to determine whether a SNe occurs in a given particle each timestep drawing from a Poisson distribution. 26 If so, the appropriate ejecta mass, metal yields (for all followed species), energy, and momentum are tabulated and directed radially from the star, and we assume the ejecta shocks within the gas in the smoothing kernel h sml (appropriately weighted) 27 around the star.  \nHowever, coupling this appropriately requires knowing whether the shock is energy or momentum conserving up to the scales we resolve. 28 To estimate this, consider the cooling radii calculated in high-resolution simulations of individual blastwaves: R cool ≈ 28pc E 0 . 29 51 〈 n cgs 〉 -0 . 43 ( Z / Z glyph[circledot] + 0 . 01 ) -0 . 18 (where E 51 ≈ 1 is the ejecta energy in units of 10 51 ergs -1 ; 〈 n cgs 〉 is the local density in cm -3 ; and Z is the local gas metallicity; see Cioffi et al.  \nand find it has little systematic effect on any quantities considered here. We also note that since the radiative acceleration is implemented as a continuous force term (as recommended in Paper I), there is no need to estimate an 'escape' or 'kick' velocity in the calculation. And we emphasize that there is no 'boost factor' in the equations above, only resolved absorption.  \n- 25 For numerical convenience, we cut off this walk at the minimum of ∼ 5 kpc or 50 times the particle smoothing length. Remaining photons are assumed to have escaped.\n- 26 It is possible, in principle, to have multiple SNe in a single particle and timestep; however, particle masses and timesteps in the dense regions where young stars live are sufficiently small that this is quite rare (occuring < 1% of the time when there is a SNe in the particle). To be avoid rare but severe exceptions, we also enforce a timestep limiter of 10 5 yr for all stars < 30Myr in age.  \n27 A particularly useful test problem is a series of SNe exploding at a constant rate from a fixed point in a thin disk (modeled as an infinite plane) with a low-density atmosphere out of plane. At high resolution, the correct solution shows the energy and momentum 'venting' out of plane. In the case where the disk is poorly resolved (less than one smoothing length in height), if the SNe energy/momentum coupling among the neighbor particles/cells is weighted by the standard SPH smoothing kernel (effectively mass-weighted coupling), then almost all the momentum goes into the disk plane and artificially drives an expanding ring. We therefore chose to assign the weights according to the fraction of the surface area subtended by each particle as seen by the SNe; this recovers the correct behavior even in the limit of low resolution.  \n28 Unlike the ultra-high resolution isolated galaxy simulations in Paper II, we cannot be sure we always resolve the energy-conserving phase of SNe expansion in these simulations.  \n1988). When h sml glyph[lessmuch] R cool, the full ejecta (shocked) kinetic energy is coupled as thermal energy (with the ejecta momentum included as a momentum flux, and the ejecta mass and metals as mass/metal fluxes into the gas). Otherwise, at coupling, a fraction of the initial ejecta energy is instead converted from energy into momentum as would occur within the un-resolved cooling radius so the coupled momentum is p = p ej √ [ M ejecta + M enc ( < R cool ) ] / M ejecta, and an additional fraction of the shocked thermal energy is cooled away before being allowed to artificially do any work (according to E thermal , shocked ∝ ( R / R cool ) -6 . 5 ; see Thornton et al. 1998). We stress that this is very different from assuming the SNe energy goes directly into a wind, or from turning off cooling in the gas. We are simply accounting for the possibility of an unresolved SedovTaylor phase and depositing the appropriate momentum , not just thermal energy, in the ambient medium in that case.  \nAs noted in the text, at high redshifts, the progenitors of massive galaxies ( glyph[greaterorsimilar] 10 12 M glyph[circledot] ) are not as well mass-resolved, given our particle masses of ∼ 10 4 M glyph[circledot] . At z ∼ 6, the progenitor of a z = 0, 10 12 M glyph[circledot] halo has M halo ∼ 10 10 M glyph[circledot] and so should have M ∗ ∼ 10 7 M glyph[circledot] , just ∼ 1000 particles. As such, the details of how momentum from SNe is coupled into a kernel of ∼ 100 particles can have significant effects on the entire baryonic galaxy. For example, we see significant changes to the SFR at z glyph[greaterorsimilar] 4 if we 'cap' the momentum input at a modest value p ≈ 20 p ej. However at lower redshifts (or in simulations with better mass resolution, such as our dwarf galaxies) this has no systematic effect (as shown in Fig. 6, where we substantially vary both the power-law scaling of the SNe momentum above and the coupling weights within the SPH kernel, and see only small effects on the final stellar mass). Because of the scaling of remnant momentum with entrained mass, the key criterion is whether the mass resolution of the simulation is glyph[greatermuch] M ejecta. In practice we find that this explicit accounting for the SNe remnant momentum has little effect when particle masses are glyph[lessorsimilar] 1000 M glyph[circledot] . We have also re-run all our simulations using the alternative cooling radius estimate from Chevalier (1974): R cool ≈ 58pc E 0 . 32 51 〈 n cgs 〉 -0 . 16 ( P gas / 10 -12 dyncm -2 ) -0 . 2 . This makes little difference to the stellar masses and star formation histories here (there is no significant difference in their consistency with observations); the primary effect is from dropping the metallicity dependence (leading to slightly less efficient feedback in low-metallicity, poorly-resolved regions).  \nStellar winds are algorithmically nearly identical to SNe, except they occur continuously. We tabulate the wind mass, metal, energy, and momentum fluxes (as a function of stellar age and metallicity), and inject these into the neighboring gas identically to the SNe. 29  \nIn both cases, we include the relative gas-star particle velocities added to the wind/ejecta velocity centered on the star in calculating the initial ejecta momentum and energy fluxes, but this has very little effect in star-forming systems (since massive star winds and SNe ejecta are fast compared to the relative velocity of stars). This can, however, be significant in old stellar populations when AGBejecta (with wind launching velocities glyph[lessorsimilar] 10kms -1 ) dominate the mass loss.  \n29 Purely for numerical convenience, we find it useful to limit the timesteps on which this occurs, so that a fixed fraction (say, ∼ 1%) of the particle mass is lost per 'feedback step' (which may be longer or shorter than the timesteps on which the star particle dynamics are evolved). Testing this we see it has no effect on our results, but is significantly less expensive computationally than invoking the feedback routine every dynamical timestep.", 'APPENDIX B: SIMULATION NUMERICAL DETAILS': 'As noted in the text, these runs adopt the P-SPH formulation of TreeSPH in GIZMO, which features many improvements to SPH and has been tested in a wide range of problems (listed in the main text). We describe the most important differences between our numerical method and previous widely-used algorithms below.', 'B1 SPH Formulation': "The simulations use the Lagrangian 'pressure-entropy' formulation of the SPH equations developed in Hopkins (2013). This formulation derives the SPH equations exactly from the particle Lagrangian and manifestly conserves momentum, energy, angular momentum, and entropy (in the absence of sources/sinks), and also eliminates the artificial 'surface tension' error term which appears at contact discontinuities in previous ('density-energy' or 'densityentropy') formulations of SPH (see also Saitoh & Makino 2013). 30 The pressure-entropy formulation dramatically improves the behavior of fluid mixing instabilities (e.g. the Kelvin-Helmholtz and Rayleigh-Taylor instabilities), and eliminates most of the known differences between the results of grid and SPH methods for the problems of interest here (for discussion of these in historical SPH implementations, see Agertz et al. 2007; Price & Federrath 2010; Bauer & Springel 2012). For extensive numerical tests, see Hopkins (2013). 31", 'B2 Artificial Viscosity & Entropy': "In SPH, the algorithm is inherently inviscid and some artificial viscosity is necessary to capture shocks. We adopt the 'inviscid SPH' viscosity prescription with higher-order switches from Cullen & Dehnen (2010). This allows for excellent shock-capturing, while reducing the viscosity to identically zero away from shocks. A wide range of tests of this algorithm are presented therein. This viscosity treatment allows accurate treatment of sub-sonic turbulence down to Mach numbers glyph[lessorsimilar] 0 . 1 while simultaneously accurately capturing shocks with Mach numbers over ∼ 10 4 . 32  \nSimilarly, SPH is inherently dissipationless, so a mechanism is also needed to generate mixing entropy in shocks. We implement this following Price (2008), using the same higher-order dissipation switch from Cullen & Dehnen (2010): this ensures that entropy exchange only occurs in crossing, sub resolution-scale flows with  \n30 As noted in the text, the choice of 'density-entropy' or 'pressureentropy' SPH makes relatively small differences to the predictions here. We have also re-run a few simulations using the 'pressure-energy' form of SPH (in which the internal energy, rather than entropy, is the explicitly followed variable). In adiabatic flows with a constant timestep the pressure-entropy and pressure-energy forms are identical to machine accuracy; with adaptive timesteps, the error reduction in the latter formulation is slightly better (poorer) in cooling (adiabatic) steps. We confirm this makes little difference to our results.  \n31 Note that non-Lagrangian schemes, in particular, have severe difficulties accurately propagating strong blastwaves and can lead to 'selfacceleration' of particles in some regimes (see Hopkins 2013). The 'traditional' method in GADGET, for example, (by which we refer to the implementation in Springel & Hernquist 2002) is Lagrangian, but adopts the 'density' formulation of SPH, which introduces the problems with fluid mixing and contact discontinuities noted above.  \n32 Wehave compared (Kereš et al., in prep.) a number of simulations adopting the simpler, time-dependent prescription from Morris & Monaghan (1997). In most respects this gives very similar results, but gives higher viscosity in sub-sonic turbulence, and produces some 'particle noise' (from interpenetration) in extremely strong shocks which can lead to unphysically high temperatures in particles leading the shock front.  \ndiscontinuous entropies (i.e. prevents artificially multi-valued entropies). For comparison, 'traditional' GADGET (Springel & Hernquist 2002) adopts a constant artificial viscosity following Gingold & Monaghan (1983) with a Balsara (1989) switch, and no entropy diffusion. This is far more dissipative, and smears out structure in sub-sonic turbulence as well as producing large (artificial) shear viscosities, which can lead to significant angular momentum transfer and eliminate structure in sub-sonic turbulence. 33  \nWehave also implemented the identical artificial viscosity and dissipation switches, and higher-order spatial gradient estimators, from Read & Hayfield (2012). We have re-run a couple of our simulations using these methods, making our hydrodynamic solver essentially identical to that in SPHS. This produces only small differences, with higher mean artificial viscosity in turbulent regions, and lower in 'smooth' spatially extended shocks (e.g. large-scale structure; for a comparison showing this method gives similar results to grid codes for cosmological accretion/halo gas, see Power et al. 2013).", 'B3 Thermodynamic Evolution & Timestep Criteria': "We employ a standard adaptive timestep algorithm and limiter. As shown in Saitoh & Makino (2009) and Durier & Dalla Vecchia (2012), in problems with high Mach number flows, adaptive timesteps (without a limiter) can lead to errors if particles with long timesteps interact suddenly mid-timestep with those on much shorter timesteps. Fortunately this is easily remedied by our timestep limiter, identical to that in Durier & Dalla Vecchia (2012). At all times, any active particle informs its neighbors of its timesteps and none are allowed to have a timestep > 4 times that of a neighbor. Whenever a timestep is shortened (or energy is injected in feedback) particles are forced to return to the timestep calculation. The limiter is not included in the 'traditional' (Springel & Hernquist 2002) GADGET.  \nWe have confirmed the importance of this limiter in our simulations: if it is absent, a small number of particles in explosive blastwaves generate large energy conservation errors which can artificially over-heat under-dense regions of the IGM. Provided the switch is included it makes no difference if we restrict the timestep ratio between neighbors to 2 , 4 , 6 or 8.", 'B4 Smoothing Kernel': "We adopt a quintic (fifth-order) spline kernel, with neighbor number designed to optimally resolve sound waves down to a wavelength ≈ h , the 'core radius' of the kernel. 34 The kernel size is  \n33 The first draft of this paper used the artificial viscosity prescription exactly from Cullen & Dehnen (2010), which produces slightly higher than desired viscosities in cosmological accretion shocks, and as implemented in our leapfrog time-integration scheme could lead to incorrect applications in the rare timesteps where the artificial viscosity term was larger than the cooling term. In all subsequent drafts this is fixed (all simulations were rerun) with an additional limiter similar to that suggested in Hu et al. (2014), and an explicit requirement that the viscosity only grows when ∇· v < 0. This does not qualitatively or systematically change any of our conclusions, but does result in small differences in detail.  \n34 This is calculated explicitly in Dehnen & Aly 2012. Note that it is meaningless to speak of the resolution scale h as the 'full size' of the kernel out to some fixed number of neighbors for higher-order kernels, since this will change with the exact functional form of the kernel adopted even while the actual resolution scale is identical. The maximum kernel radius of compact support, and equivalently radius at which gravity reduces to exactly that of a point mass, is ≈ 3 h . For more discussion see Price (2012).  \nadaptive (following the approximate 'fixed mass in kernel' prescription in Springel & Hernquist 2002), enclosing ∼ 64 neighbors. This choice is the 'optimal' spline kernel suggested by a wide range of tests in Hongbin & Xin (2005), Dehnen & Aly (2012), and Hopkins (2013). However we have also experimented with the Wendland kernels in Dehnen & Aly (2012) and triangular kernels in Read et al. (2010) and see no significant improvements up to neighbor numbers ∼ 500. For comparison, the traditional GADGET kernel is a cubic spline. This becomes unstable outside the range ∼ 30 -50 neighbors; within this range, the 'effective resolution' of the kernel is identical to that adopted here, but kernel errors are larger by nearly an order of magnitude.", 'B5 Gravity': "The gravity solver follows the GADGET-3 hybrid tree-particle mesh (Tree-PM) method. However, we have modified this to allow for adaptive gravitational softenings 35 and to more accurately symmetrize the force between interacting particles with different softening, following the fully Lagrangian method in Price & Monaghan (2007); this manifestly maintains conservation of energy, momentum, and angular momentum. In 'traditional' GADGET, softenings are not adaptive, and pairwise interactions are simply smoothed by the larger of the two particle softenings.  \nWe have also modified the softening kernel as described therein (see also Barnes 2012) to represent the exact solution for the potential of the SPH smoothing kernel. With this change, the softening no longer represents non-Newtonian gravity; rather, the gravitational force is exactly Newtonian on small scales, but for particles which are not point masses but represent the extended mass distribution represented by the SPH kernel (matching the assumption made in the hydrodynamic equations). 36  \nWe have also tested most of our 'standard' simulations with fixed gravitational softening lengths for the baryons equal to the minimum values in Table 1, and for the high-resolution darkmatter particles 37 equal to ( 20 , 20 , 50 , 150 , 150 , 150 , 200 ) pc for  \n35 With this choice, the gravitational softening of the gas follows the SPH smoothing length (so there are always ∼ 100 particles inside the full softening extent), with the minimum softening/smoothing lengths set to the values in Table 1 (fixed in physical units). These minima are chosen for numerical efficiency and the appropriate mass resolution, but also ensure that the 'hardest' scatterings of baryons off dark matter particles are weaker than the numerical errors of the long-range gravitational forces (much smaller than the real rms accelerations from irregular structure in the disks). If adaptive softening is used for other particle types, an equivalent smoothing length is computed separately for each type: each dark matter particle is assigned a smoothing length enclosing ∼ 100 dark matter neighbors, and so on. Since stars (unlike gas and dark matter) are not actually volume-filling, it is not clear if it is more physically accurate to use adaptive softenings or fixed softenings (representing e.g. open cluster extents); in these simulations (which resolve ISM structure) we prefer the latter, but find that the choice makes very little difference in practice.  \n36 This makes the force softening length in Table 1 similar to a Plummer equivalent softening: the force becomes exactly that of a point mass (as opposed to an extended mass distribution) at ≈ 2 . 8 glyph[epsilon1] and deviates by ∼ 10% from a point mass at ≈ 1 . 5 glyph[epsilon1] .  \n37 In all our runs, the high-resolution region zoom-in region is sampled with high-resolution dark matter particles, whose properties are shown in Table 1. Outside of this region there is a hierarchy of progressively more massive dark matter particles (generally in 'steps' of factors of 8 in particle mass). These particles always have much larger softenings than the highresolution particles, generally > kpc, chosen to approximately match the rms particle separation in the 'boundary' region where these particles meet their higher-resolution counterparts.  \nruns ( m09, m10, m11, m12v, m12q, m12i, m13 ). The latter corresponds very roughly to the dark matter softening kernel extent matching the mean inter-particle spacing at ∼ R vir / 4. Of course, with fixed softenings matched to Table 1, the baryons will be undersoftened (i.e. can undergo hard scatterings) when they are very diffuse (in e.g. the inter-galactic medium), but we are most concerned with accuracy in the dense regions in the actual galaxy we are simulating, where smaller softenings are more accurate. In any case this produces only small differences ( glyph[lessorsimilar] 10% in stellar mass at z = 0). Re-running a couple simulations with fixed softenings but dark matter softenings equal to the smaller minimum values in Table 1 (chosen to match the kernel extent and inter-particle spacing at ∼ 10 kpc in each simulation) similarly produces little effect. We have also re-run several of our simulations using identical softenings for both the baryonic and dark matter particles (with that softening taken to be about the geometric mean of the two values in Table 1). The differences are again small, comparable to our slightly lower-resolution runs; if we compare to a run with different baryon/dark matter softenings but the baryonic softening matched to these new runs, the differences are almost completely eliminated. There is almost no difference, between the models above, in the predicted central dark matter profiles (Oñorbe et al., in prep).", 'B6 Domain Decomposition & Parallelization': 'The simulation architecture has been heavily optimized from previous versions of the code. Gravity is still solved with a TreePM algorithm, with nested PM grids solving the large-scale forces while the tree is used for small-scale interactions. But the tree walk, domain decomposition, feedback routines, and SPH density and hydrodynamic force calculation have all been optimized substantially relative to the GADGET-3 implementation (increasing the memory requirement by a factor of ∼ 3 -5, but decreasing run-times and load-imbalances by a factor ∼ 5). The code has also been optimized for hybrid OpenMP+MPI application, allowing near-ideal strong scaling to ≈ 256 cores and modest gains to 512 cores in pure-MPI mode, and positive strong scaling (a decrease in wall-clock time for runs of fixed size, albeit less than ideal) to > 1000 cores in hybrid mode for runs presented in this work.', 'APPENDIX C: TESTING IN IDEALIZED, HIGH-RESOLUTION GALAXY SIMULATIONS': "In a series of papers (Paper I-Paper III), we have extensively studied and tested the feedback models used here in even higherresolution simulations of idealized (non-cosmological) individual model galaxies. Our dwarf galaxy simulations here are run at essentially the same resolution as the 'ultra-high resolution' dwarf galaxy models in these earlier papers, so we can safely apply the same feedback models. However, as noted in the text, in the simulations of massive halos ( glyph[greaterorsimilar] 10 12 ) here, we are forced to lower resolution (comparable to the 'high-resolution,' but not 'ultra-high resolution' simulations therein); it is therefore important to check the results of the isolated galaxy simulations at lower resolution as well.  \nTo this end, we have re-run the 'HiZ' simulation from Paper I-Paper III, with the identical code used for the cosmological simulations here. This was the most massive system considered therein, a disk with properties typical of star-forming galaxies at z ∼ 2 -4. The halo, stellar bulge, stellar disk, and gas disk have masses M halo = 1 . 4 × 10 12 M glyph[circledot] , M bulge = 0 . 7 × 10 10 M glyph[circledot] , M ∗ , disk = 3 × 10 10 M glyph[circledot] , and M gas = 7 × 10 10 M glyph[circledot] , with scale-lengths for the gas and stellar disk h ∗ , disk = 1 . 6kpc and h gas = 3 . 2kpc. We  \n-1  \nFigure C1. Tests of our methodology in idealized simulations of a single, isolated (non-cosmological), gas-rich massive disk galaxy. Top: SFR versus time. Middle: Wind mass-loading (unbound mass M wind versus total mass in new stars formed since the beginning of the simulation). Bottom: Distribution of gas densities at fixed time ≈ 200Myr after the beginning of the simulation. We separate different phases by temperature: cold gas ( T < 5000K), warm gas (5000K < T < 10 5 K) and hot gas ( T > 10 5 K; the bimodal distribution reflects low-density, high volume-filling factor material which has escaped the disk, and hot bubbles actively heated by SNe within it). In each, we compare four runs with identical initial conditions (different linestyles, as labeled). (1) A model run with our standard numerical method, at resolution about equal to our zoom-in simulations of massive galaxies. (2) A run using the density-entropy form of SPH, with our improvements to the numerical method from § B removed. (3) A run using the exact same feedback algorithms used in Paper II-Paper III, without the optimizations for cosmological runs described in § A. (4) A run at the ultra-high resolution from those papers, with much better spatial and mass resolution than can be achieved in cosmological runs. We see very little difference between the runs, suggesting both that our results are stable with respect to resolution and that the changes to the code and resolution for cosmological simulations do not fundamentally alter our conclusions from the previous work.  \n<!-- image -->  \nre-run this at two resolutions: first, the 'ultra-high resolution' level used in our earlier papers, with force softening ≈ 3 pc, and particle mass ≈ 800 h -1 M glyph[circledot] . Second, at resolution about equal to our massive cosmological simulations, force softening ≈ 10pc and particle mass ≈ 5 × 10 4 M glyph[circledot] .  \nThe results are shown in Fig. C1. We specifically compare the SFR and wind mass-loading versus time, and the phase distribution of the gas at a time t ≈ 200Myr when the galaxy has reached a quasi-steady state. The mass loading is defined as the instantaneous ratio of total wind mass (defined as the gas mass which has positive Bernoulli parameter - i.e. would escape in the absence of additional forces or cooling - with outward radial velocity > 100kms -1 ) to total stellar mass formed since the beginning of the simulation. In all cases, the steady-state SFR, wind mass-loading, and gas phase distribution are similar (to within a factor < 2). The largest differences appear at very early times, and are generally artifacts of the initial conditions (in which the disk is supported by thermal pressure, before a turbulent cascade is established); these should not appear in our cosmological simulations.  \nWe have also re-run this second 'cosmological resolution' simulation with the identical feedback implementation used in Paper II (as opposed to the updated algorithms described above), to examine whether the improvements to the code would substantially change our previous results. We see very little difference. Since the total energy and momentum inputs are fixed to the same levels by the initial mass function, and no fundamentally new source of feedback has been introduced, this is a reassuring indication that the results are not especially sensitive to the purely numerical details of the algorithmic implementation of feedback.  \nAs a final comparison, we re-run the simulation using the identical feedback implementation, but adopting the 'density-entropy' form of SPH, with a constant artificial viscosity, no entropy diffusion, a cubic spline kernel, and constant (non-adaptive) gravitational softenings. This is the form of SPH which exhibits the largest and most severe differences from grid codes in test problems (especially involving fluid mixing). But consistent with our comparison in Fig. 6, we see almost no difference between this simulation and our new SPH run. This is consistent also with our comparison in Hopkins (2013), where we perform a similar comparison with a model dwarf (SMC-mass) galaxy. At those masses we saw almost no difference in SFR, and ∼ 50% differences in the wind massloading owing to purely numerical effects, mostly the effect of fluid mixing instabilities altering the cooling time of hot SNe 'bubbles.' In the HiZ galaxy model considered here, these differences are further minimized because the hot gas is not the dominant driver of outflows."}