The Feedback In Realistic Environments FIRE project explores feedback in cosmological galaxy formation simulations. Previous FIRE simulations used an identical source code FIRE1 for consistency. Motivated by the development of more accurate numerics  including hydrodynamic solvers gravitational softening and supernova coupling algorithms  and exploration of new physics e.g. magnetic fields we introduce FIRE2 an updated numerical implementation of FIRE physics for the GIZMO code. We run a suite of simulations and compare against FIRE1 overall FIRE2 improvements do not qualitatively change galaxyscale properties. We pursue an extensive study of numerics versus physics. Details of the star formation algorithm cooling physics and chemistry have weak effects provided that we include metalline cooling and star formation occurs at higherthanmean densities. We present new resolution criteria for highresolution galaxy simulations. Most galaxyscale properties are robust to numerics we test provided 1 Toomre masses are resolved 2 feedback coupling ensures conservation and 3 individual supernovae are timeresolved. Stellar masses and profiles are most robust to resolution followed by metal abundances and morphologies followed by properties of winds and circumgalactic media. Central kpc mass concentrations in massive gtLSUBSUB galaxies are sensitive to numerics via trappingrecycling of winds in hot haloes. Multiple feedback mechanisms play key roles supernovae regulate stellar masseswinds stellar massloss fuels late star formation radiative feedback suppresses accretion on to dwarfs and instantaneous star formation in discs. We provide all initial conditions and numerical algorithms used.
['methods: numerical', 'stars: formation', 'galaxies: active', 'galaxies: evolution', 'galaxies: formation', 'cosmology: theory', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'Astrophysics - Instrumentation and Methods for Astrophysics']
FIRE2 simulations physics versus numerics in galaxy formation
{'FIRE-2 Simulations: Physics versus Numerics in Galaxy Formation': "Philip F. Hopkins glyph[star] 1 , Andrew Wetzel 1 , 2 , 3 † , Dušan Kereš 4 , Claude-André FaucherGiguère 5 , Eliot Quataert 6 , Michael Boylan-Kolchin 7 , Norman Murray 8 , Christopher C. Hayward 9 , Shea Garrison-Kimmel 1 , Cameron Hummels 1 , Robert Feldmann 6 , 10 , Paul Torrey 11 , Xiangcheng Ma 1 , Daniel Anglés-Alcázar 5 , Kung-Yi Su 1 , Matthew Orr 1 , Denise Schmitz 1 , Ivanna Escala 1 , Robyn Sanderson 1 , Michael Y. Grudi'c 1 , Zachary Hafen 5 , Ji-Hoon Kim 12 , Alex Fitts 7 , James S. Bullock 13 , Coral Wheeler 1 , T. K. Chan 4 , Oliver D. Elbert 13 , Desika Narayanan 14  \n1 TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA  \n- 2 Carnegie Observatories, Pasadena, CA 91101, USA\n- 3 Department of Physics, University of California, Davis, CA 95616, USA\n- 4 Department of Physics, Center for Astrophysics and Space Science, University of California at San Diego, 9500 Gilman Drive, La Jolla, CA 92093, USA\n- 5 Department of Physics and Astronomy and CIERA, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA\n- 6 Department of Astronomy and Theoretical Astrophysics Center, University of California Berkeley, Berkeley, CA 94720\n- 7 Department of Astronomy, The University of Texas at Austin, 2515 Speedway, Stop C1400, Austin, TX 78712, USA  \n8  \nCanadian Institute for Theoretical Astrophysics, 60 St. George Street, University of Toronto, ON M5S 3H8, Canada  \n- 9 Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA\n- 10 Institute for Computational Science, University of Zurich, Zurich CH-8057, Switzerland\n- 11 Department of Physics, MIT, 77 Massachusetts Avenue, Cambridge, MA 02139, USA\n- 12 Kavli Institute for Particle Astrophysics and Cosmology, Department of Physics, Stanford University, Stanford, CA, USA\n- 13 Department of Physics and Astronomy, University of California, Irvine, CA 92697, USA\n- 14 Department of Astronomy, University of Florida, Gainesville, FL 32611, USA  \nSubmitted to MNRAS, February 2017", 'ABSTRACT': "The Feedback In Realistic Environments (FIRE) project explores feedback in cosmological galaxy formation simulations. Previous FIRE simulations used an identical source code ('FIRE-1') for consistency. Motivated by the development of more accurate numerics - including hydrodynamic solvers, gravitational softening, and supernova coupling algorithms - and exploration of new physics (e.g. magnetic fields), we introduce 'FIRE-2', an updated numerical implementation of FIRE physics for the GIZMO code. We run a suite of simulations and compare against FIRE-1: overall, FIRE-2 improvements do not qualitatively change galaxy-scale properties. We pursue an extensive study of numerics versus physics. Details of the star-formation algorithm, cooling physics, and chemistry have weak effects, provided that we include metal-line cooling and star formation occurs at higher-than-mean densities. We present new resolution criteria for high-resolution galaxy simulations. Most galaxy-scale properties are robust to numerics we test, provided: (1) Toomre masses are resolved; (2) feedback coupling ensures conservation, and (3) individual supernovae are time-resolved. Stellar masses and profiles are most robust to resolution, followed by metal abundances and morphologies, followed by properties of winds and circum-galactic media (CGM). Central ( ∼ kpc) mass concentrations in massive ( > L ∗ ) galaxies are sensitive to numerics (via trapping/recycling of winds in hot halos). Multiple feedback mechanisms play key roles: supernovae regulate stellar masses/winds; stellar mass-loss fuels late star formation; radiative feedback suppresses accretion onto dwarfs and instantaneous star formation in disks. We provide all initial conditions and numerical algorithms used.  \nKey words: galaxies: formation - galaxies: evolution - galaxies: active - stars: formation - cosmology: theory methods: numerical", '1 INTRODUCTION': "Feedback from stars is an essential and still poorly-understood component of galaxy formation. In the absence of stellar feedback, most gas accreted into galaxies should cool rapidly on a timescale much shorter than the dynamical time, collapse, fragment, and turn into stars (Bournaud et al. 2010; Tasker 2011; Hopkins et al. 2011; Dobbs et al. 2011; Krumholz et al. 2011; Harper-Clark & Murray 2011). Cosmologically, efficient cooling inevitably results in most  \n- glyph[star] E-mail:phopkins@caltech.edu\n- † Caltech-Carnegie Fellow  \nof the baryons turning into stars, producing galaxies much more massive than observed (Katz et al. 1996; Somerville & Primack 1999; Cole et al. 2000; Springel & Hernquist 2003b; Kereš et al. 2009), regardless of the details of star formation in the simulation (White & Frenk 1991; Kereš et al. 2009).  \nHowever, the observed Kennicutt-Schmidt (KS) relation implies that gas consumption timescales are long ( ∼ 50 dynamical times; Kennicutt 1998), and giant molecular clouds (GMCs) appear to turn just a few percent of their mass into stars before they are disrupted (Zuckerman & Evans 1974; Williams & McKee 1997; Evans 1999; Evans et al. 2009). Observed galaxy mass functions and the halo mass-galaxy mass relation require that galaxies in-  \nTable 1. Initial suite of FIRE-2 simulations run to redshift z = 0  \nParameters describing the initial suite of FIRE-2 simulations in this paper. Each simulation contains several (in some, several dozen) galaxies in the high-resolution zoom-in region; halo and stellar properties listed refer only to the main 'target' halo around which the high-resolution zoom-in region is centered. All properties refer to our highest-resolution simulation using each initial condition. All units are physical. (1) Name of simulation.  \n- (2) M vir halo : Virial mass (following Bryan & Norman 1998) of the 'target' halo at z = 0 in simulation with baryons.  \n(3)  \nR  \nvir: Virial radius (spherical) of the main halo at  \nz  \n=  \n0.  \n- (4) M ∗ : Stellar mass (within ≤ 3 R 1 / 2 ) of the central galaxy in the target halo at z = 0.\n- (5) R 1 / 2 : Half-mass radius of stars in the central galaxy at z = 0 (see § 3).\n- (6) mi , 1000: the mass of baryonic (gas or star) particles, in units of 1000 M glyph[circledot] . Dark matter (DM) particle masses are ≈ 5 × larger, according to the universal baryon fraction.  \n(8) r conv DM : Radius of convergence in DM properties (based on the Power et al. 2003 criterion, with our calculation of where DM-only simulations converge from § 4). This is approximately the radius enclosing > 200 DM particles. We show below that convergence in DM profiles can in fact extend to much smaller radii in runs with baryons. § 4 shows that the DM force softening glyph[epsilon1] DM is much less important, as long as it is glyph[lessorsimilar] r conv DM . In our default runs, glyph[epsilon1] DM is fixed with values = 40pc for our dwarfs and = 30pc in our MW-mass (m12) systems.  \n(9)  \nNotes: Additional information on each simulation.  \nTable 2. Physics & Numerics Explored in This Paper (and Papers II & III)  \nA cursory outline of the physics and numerics explored in this paper. All 'standard' FIRE-2 simulations, including all in Table 1, are run with the identical simulation code and physics. However, to understand how physical and numerical changes influence our results, we systematically 'turn off' different physics and vary the numerical method and/or resolution in the sections listed here.  \nFigure 1. Mock images of two Milky Way (MW)-mass galaxies at z = 0 simulated using FIRE-2: ( m12i and m12f ). Each image is a mock Hubble Space Telescope u / g / r composite with a logarithmic stretch, using STARBURST99 to determine the SED of each star based on its age and metallicity and raytracing following Hopkins et al. (2005) with attenuation using a MW-like reddening curve with a dust-to-metals ratio of 0 . 4. We show face-on ( top ) and edge-on ( bottom ) images. Both form thin disks, with clear spiral structure, clear dust lanes, and visibly resolved star-forming regions. Properties of each galaxy are in Table 1.  \n<!-- image -->  \ncorporate or retain only a small fraction of the universal baryon fraction in stars and the ISM (Conroy et al. 2006; Behroozi et al. 2010; Moster et al. 2010). Observations of the intergalactic medium (IGM) and circum-galactic medium (CGM) require that many of those baryons must have been accreted into galaxies, enriched, and then expelled in galactic super-winds with mass loading ˙ M wind many times larger than the galaxy SFR (Aguirre et al. 2001; Pettini et al. 2003; Songaila 2005; Martin et al. 2010; Oppenheimer & Davé 2006), and indeed such winds are ubiquitously observed  \n(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).  \nUntil recently, numerical simulations treated stellar feedback in highly-simplified fashion and have had difficulty reproducing these observations. This is especially true of models which invoke only energetic feedback (thermal injection) via supernovae (SNe), which typically found the energy was efficiently radiated away (Katz 1992; Guo et al. 2010; Powell et al. 2011; Brook et al.  \nFigure 2. Mock galactic map, similar to Fig. 1, but as seen from within the galaxy, for m12i ( top ) and m12f ( bottom ). We ray-trace a Galactic (Aitoff) projection, as seen from a random star ∼ 10kpc from the galactic center. Individual, filamentary giant molecular cloud (GMC) complexes and young star clusters are visible, and both galaxies have a clear thin disk plus bulge morphology.  \n<!-- image -->  \n2011; Nagamine 2010; Bournaud et al. 2011). By 'turning off cooling' for some adjusted duration, as in Stinson et al. (2006); Governato et al. (2010); Macciò et al. (2012); Teyssier et al. (2013); Stinson et al. (2013); Crain et al. (2015), or directly putting in winds 'by hand' as in Springel & Hernquist (2003a); Davé et al. (2006); Anglés-Alcázar et al. (2014); Vogelsberger et al. (2014), it is possible to reproduce some of the observed galaxy properties. But this obviously does not prove that known stellar feedback mechanisms actually act in this way, nor can it predict many ISM and CGMscale properties that depend explicitly on e.g. the phase-structure of feedback-driven outflows (see Hummels et al. 2013).  \nAccurate treatment of star formation and galactic winds ul-  \ntimately requires realistic treatment of the stellar feedback processes that maintain the multi-phase ISM. Observationally, many stellar feedback processes - SNe, protostellar jets, photo-heating, stellar mass loss (O-star and AGB winds), and radiation pressure - act efficiently on the ISM (see Evans et al. 2009; Lopez et al. 2011, and references above). Simulations of either single molecular clouds/star clusters or the 'first stars,' which resolve individual stars and can treat these microphysics in detail, have universally found that the non-linear interaction of these feedback mechanisms successfully suppresses star formation, pre-process giant molecular clouds before SNe explosions (so that SNe occur in rarified environments), and generate galactic chimneys and super-bubbles  \nFigure 3. Mock images, as in Fig. 1, but for a subset of dwarf galaxies in our sample: m09 , an ultra-faint with M ∗ ∼ 10 4 M glyph[circledot] (similar to Coma Berenices, Leo IV, or Canes Venatici II); m10v , a faint dwarf with M ∗ ∼ 10 5 M glyph[circledot] (similar to Hercules or Leo T); m10q , an intermediate-mass dwarf with M ∗ ∼ 10 6 M glyph[circledot] (similar to Sextans, Carina, or Leo II); m11q , an SMCmass dwarf with M ∗ ∼ 10 9 M glyph[circledot] ; and m11v , an LMC-mass galaxy with M ∗ ∼ 2 × 10 9 M glyph[circledot] . Most have spheroidal morphologies, as is observed and as was seen in FIRE-1 (Wheeler et al. 2017). We show, m11v , the LMC-mass galaxy, both face-on and edge-on to illustrate the dramatic bar and elongated/flattened structure, similar to the actual LMC. Note that the surface-brightness scale is not the same in each image (an ultra-faint has ∼ 1000 × lower mean surface brightness than the LMC, so it would be invisible on the same scale).  \n<!-- image -->  \nthat generate fountains and super-winds (e.g. Krumholz et al. 2007, 2011; Offner et al. 2009, 2011; Harper-Clark & Murray 2011; Bate 2012; Wise et al. 2012; Pawlik et al. 2013; Muratov et al. 2013). A new generation of high-resolution galaxy-scale simulations has since emerged, which reach resolution sufficient to begin directly incorporating these physics, and to begin to resolve the multi-phase structure of the ISM (Tasker 2011; Hopkins et al. 2011, 2012c; Kannan et al. 2014; Agertz et al. 2013). For example, in these works and a series of related papers focused on isolated galaxy simulations (Narayanan & Hopkins 2013; Hopkins et al. 2013e, 2012a, 2013b,d,a,c), the authors showed in isolated galaxy simulations that the combination of multiple feedback processes together produce a quasi-steady state ISM, in which GMCs form and disperse rapidly, with turbulence, phase structure, GMC properties, a KS law, and galactic winds in reasonable agreement with observations.  \nMotivated by the success and predictive power of these simulations, in Hopkins et al. (2014) we introduced the Feedback In Realistic Environments (FIRE) project. 1 The FIRE code synthesized the physics and numerical methods developed in the previous work (with relevant improvements) into a single code suitable for high-resolution cosmological simulations of galaxy formation. These simulations explicitly treat the multi-phase ISM with heating and cooling physics from gas at a range of temperatures T ∼ 10 -10 10 K, star formation restricted only to self-gravitating, self-shielding, molecular, high density ( n H glyph[greaterorsimilar] 5 -50cm -3 ) gas, resolution reaching ∼ 250 M glyph[circledot] or ∼ 0 . 5pc, and (most importantly) explicit treatment of stellar feedback including the energy, momentum, mass, and metal fluxes from SNe Types Ia & II, stellar massloss (O-star and AGB), radiation pressure (UV and IR), and photoionization and photo-electric heating. All stellar evolution and feedback inputs are taken directly from stellar evolution models, without subsequent 'parameter tuning.'  \nIn a series of papers, we have subsequently shown that cosmological zoom-in simulations incorporating these physics can reproduce a diverse range of galaxy properties at a wide range of redshifts, including stellar masses, star formation histories (SFHs) and the galactic 'main sequence' (Hopkins et al. 2014; Sparre et al. 2017; Feldmann et al. 2016); metallicities and metal abundance ratios in both 'standard' and r-process elements (Ma et al. 2016a; van de Voort et al. 2015); detailed morphological and kinematic structure of thin and thick disks (Ma et al. 2017b); rotation curves and morphologies of Milky Way-mass galaxies (Wetzel et al. 2016); observed satellite mass functions, rotation curves/kinematics, and cusp/core structure of dwarfs (Oñorbe et al. 2015; Chan et al. 2015; Wheeler et al. 2015, 2017; Wetzel et al. 2016); abundance gradients (El-Badry et al. 2016; Ma et al. 2017a); neutral hydrogen absorption in the CGM (Faucher-Giguere et al. 2015; Faucher-Giguère et al. 2016; Hafen et al. 2017); galactic outflows (Muratov et al. 2015, 2017; Anglés-Alcázar et al. 2017a); star-formation properties of galactic nuclei (Torrey et al. 2017); escape fractions of ionizing photons needed for reionization (Ma et al. 2015, 2016b); and at least some of the diversity of star-forming and quiescent massive galaxies at high redshifts (Narayanan et al. 2015; Feldmann et al. 2016). There are of course a number of areas where the simulations fail to reproduce the observations: most notably, the bimodality of galaxy colors at both z = 0 (Hopkins et al. 2014) and z = 2 (Feldmann et al. 2017) - these are likely clues to important physics missing from the FIRE-1 simulations.  \nFor the sake of consistency and clarity, all FIRE simulations have used an identical source code - what we will now refer to as 'FIRE-1.' This ensured 100% identical physics and numerical choices (up to the simulation resolution and choice of the specific halo simulated) in all runs, necessary for simulation comparisons. Unfortunately, this ignores development of new, more accurate hydrodynamic solvers and gravitational force softening algorithms (see e.g. Hopkins 2015), improvements to the numerical accuracy of feedback coupling algorithms (i.e. ways to ensure machine-accurate momentum conservation in SNe coupling to gas), code optimizations that would allow higher-resolution simulations, and physics neglected in FIRE-1 such as magnetic fields,  \n1  \nSee the FIRE project website:  \nhttp://fire.northwestern.edu  \nFor additional movies and images of FIRE simulations, see:  \nhttp://www.tapir.caltech.edu/~phopkins/Site/ animations  \ncosmic rays, conduction, viscosity, optically thick radiative cooling, and more. These effects could, in principle, have large consequences for galaxy formation. For example, FIRE-1 used an improved version of the smoothed-particle hydrodynamics (SPH) method to solve the hydrodynamic equations; but it is well-known that SPH has certain low-order errors that do not converge accurately, add noise, and artificially suppress phenomena such as fluid mixing and sub-sonic turbulence (Agertz et al. 2007; Bauer & Springel 2012; Sijacki et al. 2012), potentially leading directly to large differences in cooling in massive galaxies (Kereš et al. 2012; Torrey et al. 2012). There has been considerable effort to 'fix' SPH, and FIRE-1 used the improved P-SPH methods from Saitoh & Makino (2013) and Hopkins (2013a) to reduce these errors, but some (e.g. the zeroth-order errors) cannot be entirely eliminated in SPH without de-stabilizing the method (see Price 2012). As a result, especially for fluid mixing problems, newer moving-mesh or mesh-free Godunov methods provide still greater accuracy and more rapid convergence (see Springel 2010; Hopkins 2015).  \nWe therefore introduce the 'FIRE-2' simulations: an update of the FIRE physics modules in the code GIZMO. This includes a new, more accurate hydrodynamics solver that resolves the main known issues of SPH, as well as more accurate treatments of cooling and recombination rates, gravitational force softening, and numerical feedback coupling. In this paper we present a large suite of cosmological zoom-in simulations, and compare these to our FIRE-1 results and to some basic observed galaxy properties. We find that the qualitative results from the FIRE-1 simulations are reproduced in FIRE-2. We then use these simulations to extensively explore numerical and algorithmic choices in the simulation setup, and whether these have any effect on the predictions. Some first science results from these FIRE-2 simulations have been presented in Wetzel et al. (2016); Su et al. (2017); Fitts et al. (2017).  \nThe goals of this paper are twofold. First, this is a methods and numerical/physical tests paper for the FIRE-2 simulations; we present the simulations, extensive tests of the methods, and explicitly detail all aspects of the numerical methods and algorithms. Second, we survey numerical and physical effects, e.g.: resolution (mass, spatial, and temporal), hydrodynamic solvers (SPH vs. modern Godunov methods), criteria for star formation, details of the cooling physics, and stellar feedback from radiation, winds, and SNe. For each effect we present an extensive study in simulations to understand which effects are physical, and which numerical, and where our simulations should and should not be trusted. Because we will show that feedback is the most important property determining the galaxy's formation history, a pair of companion papers will separately explore the details of the numerical implementation and physics of mechanical/SNe feedback (Hopkins et al. 2018, henceforth Paper II; this paper is Paper I) and radiative feedback (Hopkins et al., in prep., henceforth Paper III).  \nTable 1 presents the initial suite of FIRE-2 simulations studied here; Table 2 provides an 'executive summary' of our study and key conclusions; Table 3 provides a high-level overview of what 'numerical resolution' actually means in our simulations. § 2 summarizes our methods: we direct the reader to the appropriate appendices where the complete algorithmic details are presented in detail. § 3 presents a basic overview of the resulting simulations and specifically examines any differences between FIRE-1 and FIRE-2 predictions. § 4 extensively studies the effects of resolution, in mass (§ 4.1), space (§ 4.2), and time (§ 4.3). § 5 examines the effects of the hydrodynamic methods, including SPH vs. finite-volume methods and § 6 considers the effects of so-called 'artificial pressure' terms used in some (non-FIRE) simulations. § 7 studies the details  \nof cooling, chemical yields, and numerical metal-mixing terms; and § 8 considers the star formation algorithm. § 9 considers the effects of different stellar feedback physics, turned on and off in turn, to provide an indication of which feedback processes dominate, and to provide a way of quantifying the relative importance of numerical and physical (feedback) uncertainties for our results. The Appendices present additional cooling/feedback tables and algorithmic information necessary for implementing the FIRE-2 simulations.", '2 METHODS': "Here we describe our numerical methods in the FIRE-2 simulations. For further details, at the end of each sub-section below, we direct interested readers to the appendix, paper, or public code where an exact algorithmic breakdown is provided. We will study each aspect in more detail below.  \nBefore FIRE-1, a series of papers developed the numerical methods, and tested each physical addition individually using higher-resolution ISM-scale simulations and analytic solutions where possible (we refer the interested reader to Hopkins et al. 2011, 2012c,b,a, 2013a,b,c,e,d). These developments and improvements were then synthesized into the physics implemented in FIRE-1, described in Hopkins et al. (2014).  \nIn our FIRE-2 runs, the 'core' or 'baseline' physics is the same as in FIRE-1: we simply seek to improve the numerical accuracy with which we solve the relevant equations. However, the switch to a new hydrodynamics method in FIRE-2 also makes it possible in principle to include new physics such as magnetic fields: these are described here as 'additional' physics, and will be studied in separate work.  \nAs with FIRE-1, all runs denoted as 'FIRE-2' here or in any other papers (e.g. Wetzel et al. 2016; Su et al. 2017; Fitts et al. 2017) use the identical source code and physical parameters, unless explicitly labeled otherwise for comparison purposes. Of course certain numerical parameters (e.g. force softening) scale explicitly with resolution; these are provided for each simulation.", '2.1 Hydrodynamics': "A major motivation for our introduction of 'FIRE-2' is to take advantage of a new generation of accurate, mesh-free Godunov hydrodynamics methods that have been recently developed (see Gaburov & Nitadori 2011; Hopkins 2015). Because we enforced the strict requirement that all FIRE-1 simulations use the identical source code, all FIRE-1 runs used the older 'P-SPH' method (Hopkins 2013a), an improved 'pressure-energy' variant of smoothedparticle hydrodynamics (SPH). We will explore the effects of the hydrodynamic solver in our simulations in § 5. More importantly, however, the new Godunov methods allow us to accurately include more complicated plasma physics such as magnetic fields and anisotropic diffusion, which were not possible to solve accurately in P-SPH (see Hopkins & Raives 2016; Hopkins 2017).  \nWe therefore employ the meshless finite-mass (MFM) magneto-hydrodynamics solver in GIZMO. 2 This is a mesh-free, Lagrangian finite-volume Godunov method designed to capture advantages of both grid-based and particle-based methods, built on the gravity solver and domain decomposition algorithms of GADGET-3 (Springel 2005). In a series of methods papers (Hopkins 2015; Hopkins & Raives 2016; Hopkins 2016, 2017), GIZMO  \nhas been tested extensively (involving ∼ 100 distinct test problems) compared to state-of-the-art fixed grid Godunov codes (e.g. ATHENA and RAMSES; Teyssier 2002; Stone et al. 2008), movingmesh codes (e.g. AREPO; Springel 2010), and 'modern' SPH methods (e.g. P-SPH; Hopkins 2013a; Rosswog 2015; Hu et al. 2014).  \nWe emphasize that in essentially every test problem we find MFM gives more accurate results (at fixed particle number or CPU time) and faster convergence compared to state-of-the-art SPH methods, and demonstrates accuracy and convergence in good agreement with well-studied fixed-grid and moving-mesh codes. Most importantly, this includes areas where SPH has had historical difficulty, including sharp shock-capturing, fluid-mixing instabilities, magneto-rotational instabilities, and anisotropic diffusion (Ritchie & Thomas 2001; Agertz et al. 2007; Price 2008; Wadsley et al. 2008; Read & Hayfield 2012; Saitoh & Makino 2013). For some problems relevant in cosmological simulations, e.g. those with moving contact discontinuities, orbiting thin disks, supersonically shearing fluid-mixing instabilities, poorly resolved explosions, hydrostatic gravitational equilibrium or gravitational collapse, the Lagrangian nature of the method here also allows us to converge at much lower resolution compared to fixed-grid methods (Müller & Steinmetz 1995; Zingale et al. 2002; O'Shea et al. 2005; Heitmann et al. 2008; Hopkins 2015) and provides excellent angular momentum conservation (avoiding 'grid alignment' and spurious torques common in grid-based codes; Hahn et al. 2010; Byerly et al. 2014; Hopkins 2015).  \nAs discussed in Hopkins (2015), this increased accuracy and convergence rate effectively makes our simulations effectively higher-resolution (at least in terms of the spatial resolution of the hydrodynamics and its convergence), compared to FIRE-1 simulations at the same particle number.  \nFor reasons discussed in § 6, we do not adopt an artificial 'pressure floor' of any kind for hydrodynamics; unresolved fragmentation is instead explicitly treated via our star formation model.  \nReaders interested in further details of the hydrodynamic solver should consult Hopkins (2015) and the public GIZMO source code. Tests and comparisons of different hydrodynamic methods are in § 5.", '2.2 Gravity': "The N -body gravity solver is extensively detailed in Hopkins (2015); this is an improved version of the GADGET-3 Tree-PM solver and additional details can be found in Springel (2005). This solver is well-tested under a huge range of applications; we will focus here only on how this relates to the spatial or force resolution, and what this means for 'resolved scales' in FIRE. Detailed discussion and tests of resolution and force softening are presented in § 4.1-4.3.  \nBy default, as described in Hopkins (2015), the resolution of gravity and hydrodynamics are equal and the two use the same , consistent assumptions about the gas mass distribution in the simulations. Specifically, we follow Price & Monaghan (2007) and compute gravitational forces from gas particles by assuming the gas is in an extended mass distribution with the same functional form as the interaction kernel centered at the particle. 3 This means that,  \nfor gas, the gravitational force softening resolution glyph[epsilon1] gas i identically follows the inter-particle/cell separation h gas i = ∆ x , where ∆ x is the equivalent cell-length in a fixed, Cartesian mesh (i.e. the average distance between particle centers, around element i , is h gas i ), so glyph[epsilon1] gas i ≡ h gas i = 16pc m 1 / 3 i , 1000 ( n H / 10cm -3 ) -1 / 3 (where mi , 1000 is the mass resolution in units of 1000 M glyph[circledot] ). Note this definition is independent of the exact softening kernel shape or 'neighbor number.'  \nAs discussed in § 4.2, a number of studies have shown this produces 'optimal' softening in terms of (i) physical consistency with the hydrodynamics, (ii) maximizing accuracy and convergence rates, and (iii) reducing N -body integration and force errors (e.g. Merritt 1996; Bate & Burkert 1997; Romeo 1998; Athanassoula et al. 2000; Dehnen 2001; Rodionov & Sotnikova 2005; Price & Monaghan 2007; Barnes 2012; Hubber et al. 2013). Explicit tests validating the accuracy and convergence of our adaptive selfgravity implementation in GIZMO are presented in Hopkins (2015), including the Evrard (1988) polytropic collapse test, the gas (and gas+DM) Zel'dovich (1970) pancake collapse, the 'Santa Barbara Cluster' adiabatic zoom-in simulation (Frenk et al. 1999), and rotating steady-state stable disk tests.  \nFor dark matter (DM) and stars, the collisionless nature of the fluid makes the 'correct' softening ambiguous. In § 4.2 we therefore explore and compare a range of choices (both adaptive and constant). In our default simulations we set glyph[epsilon1] DM i to a constant chosen to give optimal convergence and integration accuracy based on the tests therein (essentially the largest possible value before noticeable 'over-softening' effects appear), and set glyph[epsilon1] ∗ i to a constant matched to the gas softening at the mean gas density of star formation. However we will show that these choices have little effect on any predictions here, consistent with previous studies (Bagla & Khandai 2009; Iannuzzi & Dolag 2011; Iannuzzi & Athanassoula 2013).  \nReaders interested in further details of the gravity solver and adaptive softening scheme should consult Hopkins (2015) and the public GIZMO source code. Extensive tests of force softening algorithms and choices are presented in § 4.", '2.3 Cooling': "Gas cooling is solved using a standard implicit algorithm, described in Hopkins et al. (2014), and all details are given in Appendix B. To summarize: heating/cooling rates are computed including free-free, photo-ionization/recombination, Compton, photo-electric, metalline, molecular, fine-structure, dust collisional, and cosmic ray processes, from T = 10 -10 10 K. 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). The relevant ionization states are tabulated from CLOUDY simulations including the effects of uniform but redshift-dependent background (from Faucher-Giguère et al. 2009) together with local radiation sources (every star particle is treated as a source; see the feedback description below). We account for self-shielding with a local Sobolev/Jeans-length approximation (Appendix B), calibrated in full radiative transfer experiments in Faucher-Giguère et al. (2010) and Rahmati et al. (2013).  \nWith this, high-temperature ( > 10 4 K ) metal-line excitation, ionization, and recombination rates then follow Wiersma et al. (2009a). Free-free, bound-free, and bound-bound collisional and  \nderived from the softened N -body particle Lagrangian and the expressions used in GIZMO are presented in Hopkins (2015), Appendix H2.  \n1  \nFigure 4. Basic galaxy-scale properties in a subset of FIRE-2 simulations. Each column shows a different galaxy from Fig. 1 and Table 1. Top: Star formation rate (averaged in 100Myr intervals) of the primary (most massive) galaxy in each simulation versus cosmic time. Ultra-faint galaxies ( Vc glyph[lessorsimilar] 20kms -1 ) are quenched after reionization. In more massive dwarf galaxies, SF is highly 'bursty', but it becomes less so at even higher masses and at later times (as seen in our FIRE-1 runs; Sparre et al. 2017). Second from Top: Total stellar mass versus scale factor, a = 1 / ( 1 + z ) , in the zoom-in region (we show against scale factor so that the rate of growth at early times is more clear). The value at z = 0 appears on the plot. We show all stellar mass within R vir, but this is dominated by the main galaxy, so it evolves similarly. Growth occurs rapidly at high redshifts then settles into a more steady state at late times, allowing galaxy structure to relax (El-Badry et al. 2016). Middle: Stellar Mass-Weighted Average Metallicity versus scale factor (value at z = 0 also shown). This rise is similar to the stellar mass, because these galaxies evolve on a redshift-dependent stellar mass-metallicity relation; the metallicities at each mass and redshift are nearly identical to FIRE-1 galaxies (see Ma et al. 2016a). Second from Bottom: Baryonic ( thick ) and total ( thin ) mass density profiles as a function of radius around the most massive galaxy at z = 0. Profiles are averaged in spherical shells. The dwarfs with stellar masses M ∗ ∼ 10 6 -10 9 . 5 M glyph[circledot] exhibit central 'cores' in their mass profiles, in both stars and dark matter, most prominent at M ∗ ∼ 10 9 M glyph[circledot] , where the cores extend to ∼ kpc scales, consistent with FIRE-1 results (Oñorbe et al. 2015; Wheeler et al. 2015; Chan et al. 2015). Bottom: Rotation curves of circular velocity, Vc , versus radius around the most massive galaxy at z = 0. Dwarfs exhibit slowly-rising rotation curves, while MW-mass systems have flat rotation curves with small-to-modest bulges, as in FIRE-1 (Chan et al. 2015).  \n<!-- image -->  \nradiative rates for H and He follow Katz et al. (1996) with the updated fitting functions in Verner & Ferland (1996). Photo-electric rates follow Wolfire et al. (2003), accounting for PAHs and local variations in the dust abundance. Compton heating/cooling (off the combination of the CMB and local sources) follows FaucherGiguère & Quataert (2012). Fine-structure and molecular cooling at low temperatures (5 -10 4 K) follows a pre-computed tabulation of CLOUDY runs as a function of density, temperature, metallicity, and local radiation background (see Robertson & Kravtsov 2008). Collisional dust heating/cooling rates follow Goldsmith & Langer (1978) with updated coefficients from Meijerink & Spaans (2005) assuming a minimum grain size of 10Å, and a dust temperature of 30K. Cosmic ray heating follows Guo & Oh (2008) accounting for both hadronic and Compton interactions, with a uniform cosmic ray background of ∼ 5eVcm -3 . At very high densities ( ∼ 10 10 cm -3 ), gas can become optically thick to its own cooling radiation; this is treated self-consistently following Rafikov (2007), but this is irrelevant for the simulations here (because they do not reach sufficiently high densities). Hydrodynamic heating/cooling rates (including shocks, adiabatic work, reconnection, resistivity, etc.) are computed in standard fashion by the hydro solver, then included  \ndirectly into our fully-implicit solution. A 10K temperature floor is enforced, but has no detectable effect on our results.  \nIn Appendix B, we provide fitting functions to the complete set of cooling physics above.", '2.4 Star Formation': "Gas which is locally self-gravitating, self-shielding, Jeans unstable, and above some minimum density is turned into stars using a sinkparticle approach. All details are given in Appendix C. Briefly, gas is eligible to turn into stars if and only if it meets the following criteria:  \n(i) Self-Gravitating: We require the potential energy be larger than the thermal plus kinetic energy within the resolution scale: specifically we use the sink-particle criterion developed in Hopkins et al. (2013d), α ≡ ( δ v 2 + c 2 s ) δ r / Gm gas ( <δ r ) = [ ‖∇⊗ v ‖ 2 i + ( cs , i / hi ) 2 ] / ( 8 π G ρ i ) < 1, where δ v = ‖∇⊗ v ‖ i hi and cs give the kinetic and thermal energy, respectively, within the resolution scale δ r → hi around the particle ( ⊗ is the outer product). Hopkins et al. (2013d) and many other studies (Li et al. 2005; Federrath et al. 2010; Padoan & Nordlund 2011; Padoan et al. 2012) have shown  \nthis is more useful than a density criterion (as well as more accurate in converging to the results of higher-resolution simulations), since it actually identifies gas which is collapsing under self-gravity at the resolution scale (i.e. the gas which should, physically form stars), independent of the exact spatial, mass, or density scale. This also does not allow unbound material to form stars (e.g. tidally unbound gas, or gas in strong shocks and winds which is dense, but not selfgravitating owing to large internal motions). The exact order-unity coefficient is calibrated from higher-resolution simulations of collapsing clouds/cores (see Padoan et al. 2012; Federrath & Klessen 2012), but our results are insensitive to variations.  \n- (ii) Self-Shielding (Molecular): We estimate the self-shielded ('molecular') fraction of each gas element following Krumholz & Gnedin (2011), using the local Sobolev approximation and metallicity to estimate the integrated column to dissociating radiation, and allow star formation only from the molecular component ( ρ mol = f mol ρ ). This is specifically a requirement that the gas be self-shielding, and therefore able to cool to low temperatures. Given the high n crit (see below), this criterion typically has negligible impact, since the high-density gas is typically all shielded and molecular anyways.  \n(iii) Jeans Unstable: We require a thermal Jeans mass below the maximum of the particle mass or 10 3 M glyph[circledot] in the element. This is done to ensure that any resolved, massive self-gravitating objects which should collapse coherently (as opposed to fragmenting into stellar mass-scale objects) are followed self-consistently and not simply assigned to stars (the choice of 10 3 M glyph[circledot] is designed to separate massive 'clumps' from normal very massive stars, but is not important). In practice this criterion is always easily met when other criteria are met.  \n(iv) Dense: To prevent spurious application of the criteria above, we also check that n H > n crit = 1000cm -3 (much larger than the mean galaxy density 〈 n 〉 ). This restricts star formation to dense molecular clouds fragmenting out of the background disk.  \nIf gas meets all the criteria above, we assume it turns into stars at a rate ˙ ρ ∗ = ρ mol / t ff where t ff is its free-fall time. This also comes directly from higher-resolution simulations of turbulent clouds (Padoan & Nordlund 2011; Padoan et al. 2012), as well as analytic models for star formation via turbulent fragmentation (Hopkins 2012b,a; Guszejnov & Hopkins 2015, 2016; Guszejnov et al. 2016). We stress that this is an assumption about the rate at which small, locally self-gravitating clumps (which may only be a small fraction of the dense gas mass) fragment; it does not imply the global efficiency of star formation (either on galaxy or GMC scales) is necessarily high - we find that it is self-regulated by feedback at ∼ 1 -10% per free-fall time (see Hopkins et al. 2011, 2013e, 2014) even in gas with densities glyph[greaterorsimilar] 100cm -3 , in agreement with observations (Lee et al. 2016). Recently, similar implementations to ours have also found consistent results on large scales (Semenov et al. 2016; Agertz & Kravtsov 2016), consistent with analytic expectations (Ostriker & Shetty 2011; Faucher-Giguère et al. 2013). If particles do not meet all of the criteria above, their SFR is zero. Gas particles which turn into star particles begin life as zero-age main sequence populations, with abundances and total mass inherited from their progenitor gas particle.  \nWe provide the complete set of formulae and detailed algorithmic implementation of star formation in Appendix C.", '2.5 Stellar Feedback': "Once a star particle forms, it is treated as a single stellar population, with known age t ∗ = t -t form, metallicity (inherited from its progen-  \nor gas particle), and mass (equal to its progenitor gas particle). All feedback quantities are tabulated directly without subsequent adjustment or fine-tuning - from standard stellar population models (STARBURST99; Leitherer et al. 1999) assuming a Kroupa (2001) IMF (the same as FIRE-1). 4  \nHere, we briefly summarize the feedback mechanisms. Because these are the most important and novel aspect of the FIRE simulations, we discuss the exact physics and algorithmic implementation in much greater detail in the companion papers, Paper II & Paper III. These papers present extensive tests of the algorithms, with idealized simulations of e.g. individual SNe explosions reaching resolution < 0 . 01 M glyph[circledot] and experiments using detailed radiation-hydrodynamics simulations, used to validate the exact implementations here. But for the sake of completeness, we summarize them here and provide the complete algorithms in the Appendices.  \nThe physics of stellar feedback in FIRE-2 are the same as in FIRE-1 (Hopkins et al. 2014), and the algorithms are identical up to improvements in accuracy which we explicitly detail below.  \n(i) Supernovae (Ia & II): (For details, see Appendix D.) Every timestep ∆ t , for each star particle, the tabulated SN rate as a function of star particle mass, age, and metallicity is used to determine the probability p of an event (Type-Ia and/or Type-II) occurring within the particle within ∆ t ; our mass and time resolution is such that p glyph[lessmuch] 1, i.e. we explicitly treat individual SNe explosions, rather than model their collective effects indirectly. We determine probabilistically if an event occurs within ∆ t ; if so, the appropriate ejecta mass, metal yields, energy, and momentum (also determined from the stellar evolution tables) are deposited directly in the surrounding gas around the star particle. The algorithm for deposition is constructed to ensure manifest, machine-accurate conservation of mass, metal mass, energy, and momentum, while also ensuring that the ejecta are distributed isotropically in the rest frame of the star. In Paper II we show that this is non-trivial in Lagrangian codes such as ours, where highly anisotropic gas distributions around a star particle can easily bias the momentum distribution and even violate linear momentum conservation, if the algorithm is not carefully designed to prevent this. We properly account for the relative star-gas motion (so e.g. the exact shock solution includes the initial stellar motion through the background gas). We determine the coupled momentum by computing the exact Sedov-Taylor solution for an energy-conserving spherical shock, at the coupling location (resolved separation between gas and star); if the resulting momentum exceeds the terminal momentum at which point the shock should have become radiative (equivalently, if the resolved coupling radius is larger than the cooling radius), we deposit only the momentum which would have been present when it reached that cooling radius. Paper II shows that this ensures our simulations exactly reproduce the fully-converged solutions (with resolution < 0 . 01 M glyph[circledot] ) for individual SNe explosions in high-resolution ISM simulations (once they reach the same radius as our coupling radius), independent of our resolution, for the same ambient density. We stress that we do not turn off cooling or otherwise impose any assumption about 'galactic wind driving.'", '(ii) Continuous Stellar Mass-Loss (OB/AGB-star Winds):': "4 Of course, alternative stellar evolution/IMF models may predict different feedback properties, but we will not investigate this here. In general the predicted variation is small, but for some quantities, e.g. the escape fraction of ionizing photons at high redshift, we have shown it can be important (Ma et al. 2016b).  \n(Details in Appendix D.) Similarly, stellar mass-loss is injected continuously in the gas surrounding each star particle as a function of stellar age and metallicity, with the appropriate mechanical luminosity, momentum, mass, and metal content, including both fast (O/B-star) and slow (AGB) winds, calculated directly from the stellar evolution models. The algorithm for deposition is exactly the same as for SNe, except there is an 'event' every timestep with associated ejecta mass =∆ t ˙ M wind.  \n(iii) Photo-Ionization and Photo-Electric Heating: (Details in Appendix E.) For computing radiative feedback properties, each star particle is treated as a source with an appropriate age and metallicity-dependent, IMF-averaged spectrum. We approximate the complete spectrum with a five-band treatment that includes ionizing, far-UV (relevant for photo-electric heating), near-UV, optical/near-IR, and mid/far-IR (re-radiated dust emission) photons. The background radiation owing to these sources is then locally extincted by the gas immediately surrounding the star (using a Sobolev approximation to estimate the column integrated to infinity, and extincting each band accordingly), with frequency and metallicity-dependent opacities from dust and neutral gas. The luminosity absorbed by dust (non-ionizing bands) is assumed to reradiate in the mid/far-IR band. The resulting, post-extinction luminosities are then propagated to long-range distances through an optically thin transport network (using a tree structure), to calculate an incident flux in each band at all positions. We therefore refer to this transport algorithm as the 'Locally Extincted Background Radiation in Optically-thin Networks,' or LEBRON, approximation. Since we simulate only a single small region of the Universe surrounding one galaxy in our 'zoom-in' simulations, we add to the diffuse ionizing-band flux a uniform but redshift-dependent meta-galactic background tabulated from Faucher-Giguère et al. (2009). The fluxes are then corrected for self-shielding using the same local-extinction Sobolev approximation, at the location of the gas. The resulting incident ionizing and FUV fluxes are then used to self-consistently compute the gas ionization states and radiative heating/cooling rates in our standard cooling algorithms described in Appendix B.  \n(iv) Radiation Pressure: (Details in Appendix E.) As photons are tracked according to the algorithm above, each explicitly resolved absorption transfers the appropriate photon momentum = L abs ˆ n / c (where ˆ n is the direction of ray propagation) to the gas. This automatically accounts for both direct UV/optical singlescattering, and indirect re-radiated IR photons (which can, in principle, be multiply-scattered, although this rarely occurs on the relatively coarse scales we resolve in the FIRE simulations). We stress that we do not assume any 'sub-grid' photon coupling, multiplescattering, or radiation pressure - there is no 'boost factor' anywhere in the model: in both FIRE-2 and FIRE-1, the only radiation pressure in the simulations is from explicitly resolved photon absorption. In Paper III, we show that only ∼ 1 / 2 of the total bolometric luminosity of stellar populations is absorbed at all, and, given our numerical resolution (which prevents us from resolving e.g. proto-stellar cores), the multiple-scattering IR term accounts for < 10% of the galaxy-averaged radiation pressure (it may be important, however, in dense galactic nuclei corresponding to observed systems like Arp 220).  \nWe emphasize that while quantities like SNe rates and stellar spectra are IMF-averaged, individual SNe are always discrete events (not continuous energy injection). In future work we will consider the effects of explicitly sampling the spectrum of stellar masses from the IMF (Su et al. 2018); however our preliminary re-  \nsults indicate the effects on large scales are (unsurprisingly) small compared to our IMF-averaged approach.  \nFor readers interested in reproducing our results, we provide simple fitting functions to all of our stellar evolution tabulations (and yield tables) needed for the feedback mechanisms above, in Appendix A. All details of the algorithmic implementation of mechanical feedback (SNe and stellar mass-loss) are given in Appendix D, and all details of the algorithmic implementation of radiative feedback (radiation pressure, photo-ionization, and photoelectric heating) are given in Appendix E.", "2.6 'Additional' Physics: Magnetic Fields, Conduction, Viscosity, Diffusion, Cosmic Rays, Black Holes, and more": "As noted above, a major motivation for our migration to FIRE-2, using the new MFM hydrodynamic solver, is to compare simulations including more complicated plasma physics, e.g. magnetic fields. However, for the sake of clarity and direct comparison with FIRE-1, in this paper we will focus on simulations that include our 'core' set of FIRE physics (gravity, hydrodynamics, cooling, star formation, and stellar feedback, as described above). This means that our 'default' or 'core physics only' FIRE-2 simulations use the same physics as FIRE-1, just more accurate numerical integration of those physics.  \nThe effects of additional physics will of course be the subject of their own studies. Some examples include (i) magnetic fields and (ii) anisotropic (Braginskii) conduction & viscosity (both studied in Su et al. 2017); (iii) passive-scalar turbulent eddy diffusion (e.g. metal diffusion), discussed briefly here in § 7.2 and in more detail in Escala et al. (2018); (iv) cosmic rays (Chan et al., in prep.); (v) alternative radiation-hydrodynamics (using e.g. alternative RHD solvers such as the M1 method as implemented in Hopkins & Grudic 2018 or direct integration following Jiang et al. 2014), discussed in detail in Paper III; (vi) super-massive black hole formation, accretion, and feedback (see e.g. Anglés-Alcázar et al. 2017b, for a preliminary exploration).", '2.7 Timesteps & Integration': "Our time integration scheme is discussed in detail in Hopkins (2015). Following Springel (2005) we use an adaptive powerof-two hierarchy for assigning individual timesteps for particles. As shown in Saitoh & Makino (2009) and Durier & Dalla Vecchia (2012), in problems with high Mach number flows, adaptive timesteps can lead to errors if particles with long timesteps interact suddenly mid-timestep with those on much shorter timesteps; this is remedied by requiring that, 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 of any sort) particles are forced to return to the timestep calculation. This has been tested extensively in Hopkins (2013a) and Hopkins (2015).  \nThe timestep is set by the minimum of various criteria. All particles obey limits ∆ t < 0 . 2 ( hi / | a i | ) 1 / 2 (see Power et al. 2003; here hi is the minimum of the inter-particle separation or Plummerequivalent force softening) 5 and ∆ t < 0 . 25 / |∇· v i | , where v i and a i are the total velocity and acceleration of particle i (including all sources of acceleration: e.g. feedback and hydrodynamic forces, for gas). For further safety, we always enforce a maximum timestep of ∆ a < 10 -4 a (where a is the scale-factor), but this is almost never  \nimportant. Gas elements must also obey the Courant (CFL) condition: ∆ t < 0 . 4 hi / v max sig , i , where v max sig , i is the usual maximum signal velocity between all particles interacting with i (see Hopkins 2015 for tests and details). Additional timestep criteria apply if additional fluid physics (magnetic fields, diffusion, cosmic rays, radiation) are included (see Hopkins 2017). 6 In the above equations, note hi is defined by the inter-particle separation as defined in § 4.2, which is the appropriate value for the prefactors here (see Hopkins 2015; the pre-factors would need to decrease by a factor ∼ 2 if we replaced hi with the maximum allowed neighbor distance, for example). If adaptive gravitational softening is used for collisionless (star and dark matter) particles, they must obey additional Courantlike timestep criteria given in § 4.3.3.  \nFor star particles, we additionally impose a restriction ∆ t < MAX ( 10 4 yr , t ∗ / 300 ) , where t ∗ is the age of the star; this prevents the code from 'skipping' any significant portion of stellar evolution if, somehow, a star formed in a region where the other timestep criteria allowed long timesteps (although this is very rare), and also ensures that the expectation value of the number of SNe per particle per timestep is always < 1 at our production resolution.  \nReaders interested in more details of the time integration scheme should consult the public GIZMO source code.", '2.8 Initial Conditions': "All simulations in this paper are fully cosmological 'zoom-in' simulations: a large box is simulated at low resolution to z = 0, and then the mass within and around the halo(s) of interest is identified, traced back to the starting redshift, and the Lagrangian region containing this mass is re-initialized at much higher resolution for the ultimate simulation (Porter 1985; Katz & White 1993). The initial conditions are generated with the MUSIC code (Hahn & Abel 2011), using second-order Lagrangian perturbation theory to evolve the initial conditions to redshift z ∼ 100, at which point the GIZMO simulation begins. In the simulations here, the Lagrangian high-resolution regions are defined by a convex hull including all particles within ∼ 5 R vir of the final ( z = 0) 'primary' galaxy (most massive galaxy within the high-resolution region); we have used a series of re-simulations with progressively higher resolution, including baryons, to refine the Lagrangian regions more accurately, with a target of zero low-resolution DM particles contaminating the region within ∼ 2 R vir (following Oñorbe et al. 2014). Typically, these regions include a number of smaller galaxies; however, in this paper, we exclude any galaxy with > 1% contamination (from low-resolution particles) by mass within R vir.  \nTable 1 describes the initial conditions for the initial set of halos we have simulated to z = 0. We consider a series of halos with different masses; many of these are chosen to match the halos from our FIRE-1 studies (specifically simulations first presented in Hopkins et al. 2014 and Chan et al. 2015). In all cases, the ICs are re-generated if needed to meet our strict contamination standard above. A couple of FIRE-1 ICs are not re-simulated here, because they were not generated from the MUSIC code (they were taken from older work); for consistency and clarity we will only include ICs generated in a consistent manner here. We have added new simulations here to increase our statistical sampling of halo growth  \nhistories and mass. The 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. Simulations labeled 'q' (e.g. m10q ) have more 'quiescent' halo growth histories at late times (i.e. tend to form earlier) while those labeled 'v' have more 'violent' latetime histories (tend to form later); however we stress that these all lie well within the typical scatter in such histories at each mass (for example, each 'q' history has several major mergers at high redshifts). Other labels ('i', 'f') are purely for bookkeeping. Several ICs ( m10q , m10v , m11q , m11v , m12q , m12i ) are taken from the AGORA project (Kim et al. 2014, 2016), to enable easy comparisons with a wide range of different codes. We adopt a standard, flat Λ CDM cosmology with h ≈ 0 . 70, Ω M = 1 -Ω Λ ≈ 0 . 27, and Ω b ≈ 0 . 045 (consistent with current constraints; see Planck Collaboration et al. 2014). 7  \nWe scale the resolution with simulation mass, to achieve the optimal possible mass and force resolution for each halo; we study both mass and force resolution extensively in § 4.  \nFor readers interested in more details, or reproducing our results, all initial conditions used here are publicly available. 8", '2.9 Parallelization & Runtime Requirements': "The simulations here use a hybrid OpenMP-MPI parallelization scheme with a number of optimizations specific for 'zoom-in' simulations. These are listed in Appendix G, together with explicit strong and weak scaling tests. Our improvements allow us to extend good weak scaling on production zoom-in simulations to at least > 16 , 000 CPU cores (and > 10 6 cores for large-volume simulations).  \nWith these optimizations, each high-resolution, productionquality simulations of a MW-mass galaxy with particle masses ∼ 7000 M glyph[circledot] (a few × 10 8 total particles) typically requires ∼ 10 6 cpuhrs; for our smallest dwarfs (particle masses ∼ 250 M glyph[circledot] ), their much lower star formation efficiencies and baryonic densities reduce this to ∼ 10 4 cpu-hrs. All simulations here were run on the XSEDE Stampede, Comet or NASA Pleiades machines.  \nDetails of our code optimizations and scaling tests are presented in Appendix G.", '2.10 A Complete List of Differences Between FIRE-1 & FIRE-2': "Although they are discussed in great detail throughout this paper, for the sake of clarity we here summarize the differences between the FIRE-1 and FIRE-2 simulations, in order from most to least important.  \n- (i) More Accurate Hydrodynamic Solver: As described in § 2.1, FIRE-2 uses the newer, more accurate mesh-free finitevolume Godunov-type MFM method to solve the hydrodynamic equations. FIRE-1 used the older 'pressure-energy' SPH ('PSPH') method. In § 5 we show how this affects our results; while the differences are generally second-order, this appears to be the single change with the largest effects on our predictions.  \n1  \nFigure 5. Comparison of galaxy properties and formation histories in FIRE-1 versus FIRE-2, as in Fig. 4. We show galaxies for which the identical halo is a member of the 'core set' of both FIRE-1 and FIRE-2 simulations. FIRE-2 combines a more accurate hydrodynamic method, higher resolution, a more accurate numerical algorithm for depositing supernova ejecta into gas around explosions, and updated cooling tables (for a complete list of changes, see § 2.10). Nevertheless the results are qualitatively similar in every property that we examine here. We do see some quantitative differences. For dwarf galaxies, we find slightly lower stellar masses, because of the updated photo-heating tables. Massive galaxies show somewhat higher masses and central rotation velocities, because of enhanced mixing, which occurs because our more accurate hydrodynamic method changes the cooling and efficiency of wind escape in 'hot halos' at late times. The enhanced 'burstiness' in FIRE-1 m11v occurs because it was run with ∼ 10 × lower resolution as compared to FIRE-2. We examine each of the numerical aspects of the method in detail below.  \n<!-- image -->  \n(ii) Manifestly-Conservative Supernovae Ejecta Distribution: In Appendix D and Paper II, we describe in explicit detail how, algorithmically, we distribute the products (mass, metals, momentum, and energy) of mechanical feedback (SNe and continuous stellar mass loss) from star particles into the surrounding gas particles. As discussed there, the FIRE-1 runs used a simpler algorithm, which can, in situations where the gas elements surrounding a star are highly disordered, produce a distorted (anisotropic) deposition (e.g. biasing the momentum deposition so it is not deposited symmetrically in the rest frame of the star, violating linear momentum conservation). We stress that the FIRE-1 algorithm still ensured the mass and energy of ejecta were exactly conserved; the issue comes with the spatial/vector distribution of the ejecta (momentum conservation). We have developed a novel scheme in FIRE-2 which eliminates this numerical bias and ensures manifest conservation. This difference generally has small effects, but does appear to influence the central stellar masses/densities of massive galaxies, and because the error term in the older implementation was resolutionindependent, it actually can influence galaxy morphologies more dramatically at the highest resolutions.  \n(iii) More Accurate Photo-Ionization Heating: In Appendix E and Paper III, we describe our treatment of radiation transport in explicit detail, including UV/optical/IR radiation pressure, photoelectric heating, and photo-ionization heating. The method and source terms are almost entirely identical between FIRE-1 and FIRE-2. In our default treatment of photo-ionization heating specif-  \nically, we conduct a gas neighbor search around each star particle, consuming ionizing photons (as we move outward) using a local Stromgren approximation until the photon budget is exhausted. However in FIRE-1, the search was simply terminated at the boundary of the local computational domain - any remaining photons were lost. In FIRE-2, any remaining photons are propagated via the long-range tree-based radiative transfer method in Appendix E. The fraction of photons affected is small since the vast majority are absorbed locally, and so this produces weak or negligible differences on galactic scales (nearly undetectable except in small dwarfs), but it eliminates the explicit domain-dependence of local HII regions.  \n(iv) Removal of 'Artificial Pressure' Terms: In FIRE-1, we included an artificial numerical 'pressure floor' for cold gas in the ISM, following the approach in e.g. Robertson & Kravtsov (2008) (designed to suppress collapse of any gas resolved with < 4 thermal Jeans lengths). As discussed in § 6, this term is now (a) redundant with our star formation prescription, and (b) potentially unphysical, as it fails to conserve energy and can introduce noise or suppress real fragmentation. We therefore include no such terms in FIRE2, but instead follow standard practice in the star formation community and rely on the sink-particle (star formation) criterion to treat unresolved fragmentation (see Federrath et al. 2010, for discussion). In § 6 we show the removal of these terms has no effect except to eliminate some obviously unphysical resolution-scale artifacts in the cold gas, as expected.  \nFigure 6. Mock images, as in Fig. 1, comparing FIRE-1 ( left ) and FIRE-2 ( right ) versions of the same galaxy on the same scale. Top: Dwarf galaxy ( m10q ). Because their morphologies are irregular or spheroidal, they are similar independent of numerical details. Middle & Bottom: MW-mass galaxy ( m12i ) seen face-on and edge-on. Qualitatively, the morphologies are similar. The FIRE-2 run is higher-resolution, which translates to a slightly thinner thin disk and a more extended, low surface-brightness outer disk.  \n<!-- image -->  \n(v) Updated Cooling Tables & SNeII Yields: The physical mechanisms of stellar feedback, and assumptions about stellar evolution, are the same between FIRE-1 and FIRE-2. This means SNe rates (Ia and II), wind mass loss rates and kinetic luminosities, bolometric luminosities and luminosities in different bands, yields, etc., are the same. We have made one minor update: in FIRE-1, we used the SNe II yields of Woosley & Weaver (1995); however, it is widely known that these older models significantly under-predict the observed yields in Mg and Ne, and we confirmed this in Ma et al. (2016a). We have therefore updated this to the more recent Nomoto et al. (2006) yields, which remedies this issue. We stress, though, that for all other species (especially C and O, which constitute most of the metal mass and are the dominant coolants), the IMF-averaged yield is within ∼ 10% of Woosley & Weaver (1995). Since Mg and Ne are negligible coolants, this has no detectable effect on our main results. Similarly, the cooling physics is the same in FIRE-1 and FIRE-2. However we have updated some of the actual fitting functions used to compute the cooling functions (specifically for the recombination rates, photo-electric heating including PAHs, optically-thick cooling, and dust cooling), to match more accurate cooling tables made public since FIRE-1 was developed.  \nFigure 7. Stellar mass-halo mass relation for FIRE-2 simulations (colored points) at z = 0. Stellar masses and halo virial masses are defined as in Table 1, for all resolved, uncontaminated halos (116 galaxies total; see text, § 3). Large points show the 'primary' (most massive) galaxy within the zoom-in region, in each simulation (different point styles). Grey triangles show FIRE-1 simulations. While individual galaxies may differ in mass, the effects are primarily stochastic: the two agree well on average. We compare observational estimates as labeled; black dotted lines show the observationally estimated ∼ 95% intrinsic scatter (see text). Within the scatter and systematic variations between fits, the simulations agree well with the observations at all masses.  \n<!-- image -->  \nFor the sake of transparency and clarity, a complete set of fits to the FIRE-2 stellar evolution, yield, and cooling tabulations are presented in Appendices A-B.  \n(vi) Code Optimization, Higher Resolution: For FIRE-2, we have made a number of purely numerical optimizations to the GIZMO code, to improve speed and parallelization efficiency (for details, see Appendix G). We have also re-compiled some lookup tables and re-fit cooling functions for greater precision. This has no effect on our results, of course, but it has allowed us to run new simulations at even higher resolution compared to FIRE-1.", '3 BASIC RESULTS & COMPARISON BETWEEN FIRE-1 &FIRE-2': "Table 1 summarizes all the production-quality FIRE-2 simulations run as of writing this paper. For each, we give the ( z = 0) halo virial mass, virial radius, stellar mass of the 'target' galaxy (the galaxy used to identify the initial zoom-in region), half-mass radius of the target galaxy, mass resolution of the simulation, and some values describing the 'spatial resolution' (because our simulations are Lagrangian, mass resolution is well-defined, but 'spatial resolution' is inherently variable: we discuss this in detail in § 4.2). We have considered simulations spanning a z = 0 halo mass range from M halo ∼ 10 9 -10 12 , similar to our FIRE-1 simulations. All the simulations here have been run to redshift z = 0.  \nFig. 1 shows both face-on and edge-on images of two of our FIRE-2 MW-mass systems ( m12i and m12f ), at the highest resolution we have studied ( mi , 1000 = 7). These use STARBURST99 (the same assumptions used in-code ) to compute the stellar spectra as a function of age and metallicity for each star particle, and then ray-trace through the ISM assuming a constant dust-to-metals ratio and physical dust opacities to volume-render the observed images in each band, which we use to construct a mock HST u / g / r composite image as seen by a distant observer. Fig. 2 shows images from within the galaxy: we select a random star ∼ 10kpc from  \nthe galactic center and construct a Galactic Aitoff projection of the ray-traced image from all stars in the galaxy to the mock observer. Fig. 3 shows images of several dwarf galaxies from the ultra-faint through LMC mass scales.  \nFig. 4 shows several properties of a representative subset of our simulations: the star formation rate and stellar mass versus time (archeological formation history of stars within the z = 0 galaxy); the stellar mass-weighted mean metallicity of those stars versus time; the z = 0 baryonic and total mass profiles; and the z = 0 circular velocity curve. Each property is measured for the 'target' galaxy in the simulation. Essentially all of our high-resolution simulations show qualitative behavior broadly similar to one of the galaxies plotted.  \nIt is not our intention in this paper to explore a quantitative comparison of the simulations and observations: this will be the subject of future work. For example, detailed comparison of the scaling relations of galaxy angular momentum, rotation curves, sizes and the Tully-Fisher relation can be found in El-Badry et al. (2018b,a), while Fitts et al. (2017) and Chan et al. (2018) compare the size, structural properties, and surface brightness distribution functions of dwarfs, and Garrison-Kimmel et al. (2017) the sizes of Milky Way-mass systems at z = 0 (and Ma et al. 2018b at z glyph[greaterorsimilar] 5), as a function of halo properties. We will show mass profiles, however, so that one can infer where they are sensitive to the numerical choices explored in this paper.  \nFig. 5 compares the galaxies for which we have both production-quality FIRE-1 and FIRE-2 simulations. Fig. 6 compares the visual morphology of the same galaxies. Here we can directly compare formation histories and morphologies of the same galaxy, with our improved numerical methods.  \nIn Fig. 7, we plot the stellar mass-halo mass relation for our FIRE-1 and FIRE-2 simulations, compared to observations. We identify all resolved, un-contaminated halos in the high-resolution region and plot their virial masses and the stellar mass of the central galaxy in each halo, 9 which is (as expected) always smaller than the total mass in R vir plotted in Fig. 4. The FIRE-1 results here are taken directly from Hopkins et al. (2014). We compare these predictions to recent observational constraints, from a combination of abundance matching and mass modeling. The observational fit from Moster et al. (2013) only includes galaxies with M ∗ glyph[greaterorsimilar] 10 9 M glyph[circledot] (so it is extrapolated here), but this extrapolation has been shown to reproduce well the observed local group dwarf luminosity func-  \n9 We use the HOP halo finder (Eisenstein & Hut 1998) to identify halos in Fig. 7, for the sake of consistency with the FIRE-1 results published in Hopkins et al. (2014). This combines an iterative overdensity identification with a saddle density threshold criterion to merge subhalos and overlapping halos. We define halo mass M vir and radius R vir as the Bryan & Norman (1998) virial mass/radius. We discard any halo outside the fully high-resolution region ( > 1% contamination by mass, from low-resolution particles), as well as unresolved halos (with < 5 × 10 4 DM, < 100 baryonic, or < 10 star particles), and subhalos (any halo within < 2 R vir of a more massive halo center). The exact value of these cuts makes no difference to our conclusions. We define central stellar mass as in Table 1 iteratively by first measuring the half-mass radius of all stars within a large cut (inside 15% of R vir), then taking all stars within 3 × this radius (and then re-defining the half-stellar mass radius on these stars). This eliminates all satellites we identify by visual inspection and gives results reasonably close to fitting mass profiles of the central system (a detailed mock observational study is presented in Price et al. 2017; but for an exponential disk this recovers 97% of the mass). Using a simpler cut of all stars at < 0 . 1 R vir gives similar results, except one case with a z ≈ 0 merger; but for massive galaxies the 0 . 1 R vir cut includes stars that are clearly part of the halo.  \ntions to ∼ 10 4 M glyph[circledot] (Garrison-Kimmel et al. 2016), so we consider it as well. Brook et al. (2014) combine the constraints from Behroozi et al. (2012) at high masses ( M ∗ > 10 9 M glyph[circledot] ) with local group and field dwarf constraints. Garrison-Kimmel et al. (2016) perform a similar exercise, allowing the scatter below M halo < 10 11 M glyph[circledot] to vary; we plot their best-fit median relations for a constant scatter below this mass ≈ 1dex (but note that for any scatter in the range 0 . 5 -2dex, the results are similar). We expect the scatter to vary continuously with mass, so we show the 95% inclusion contour if we take the model from Garrison-Kimmel et al. (2016) where the scatter is constant at 0 . 2dex above M halo > 10 11 . 5 M glyph[circledot] (the value favored by Behroozi et al. 2012 and Moster et al. 2013), and varies linearly with halo mass as σ dex = 0 . 2 -0 . 2 log 10 ( M halo / 10 11 . 5 M glyph[circledot] ) at lower masses (rising to ≈ 0 . 5dex at M halo ≈ 10 10 M glyph[circledot] ). 10  \nWe emphasize that although matching the full observed (2D) distribution of galaxies in stellar mass-halo mass space is equivalent to matching the observed stellar mass function (SMF), with the limited sample here, we can only test whether our simulations are consistent with being drawn from this distribution (we do not have a sufficiently large ensemble of halos to forward-model the entire SMF). However, in Wetzel et al. (2016) and with a much larger sample in Garrison-Kimmel et al. (in prep), we have sufficient statistics to forward-model the observed dwarf galaxy SMF and compare to the local (Milky Way and M31) observed satellite and local field SMFs, at stellar masses ∼ 10 4 -10 9 M glyph[circledot] . And in Ma et al. (2018a) we use an expanded sample of simulations run to high-redshift to compare the SMF (and luminosity functions) at ∼ 10 6 -10 11 M glyph[circledot] to observations at z glyph[greaterorsimilar] 5. However more quantitative comparison to the SMF of massive galaxies at lower redshifts ( z ∼ 0 -2) will require larger statistical volumes.  \nIn any case, in Figs. 1-7, we confirm the conclusions of our previous FIRE-1 studies. Although it is impossible to be exhaustive at this point, we have yet to identify any area in which the FIRE2 predictions differ at the order-of-magnitude level from FIRE-1 predictions. In future work, we will examine detailed properties of the CGM (e.g. column density distributions of different absorbers) where the hydrodynamic solver could, in principle, have a larger effect.  \nAs in FIRE-1, in FIRE-2 the dwarfs tend to have spheroidal morphologies, with relatively little rotation (see Wheeler et al. 2017). Especially around M halo ∼ 10 11 M glyph[circledot] , repeated bursts of star formation driving cycles of outflow, subsequent infall, and repeated star formation leads to 'puffing up' of the dark matter and stellar orbits, generating large cores in the dark matter profiles (Oñorbe et al. 2015; Chan et al. 2015). This also leads to expansion of the galaxy size and low surface brightness in their stellar distribution (El-Badry et al. 2016). 'Bursty' star formation predominates in low-mass dwarfs and high-redshift progenitors of massive galaxies (Sparre et al. 2017; Faucher-Giguère 2018), and in galactic nuclei (Torrey et al. 2017). In FIRE-1 and FIRE-2, MW-mass galaxies develop a clear disk+bulge morphology, with large thin stellar disks (despite the presence of strong feedback; Ma et al. 2017b), with clear spiral structure, pronounced radial metallicity and age gradients (Ma et al. 2017a). The stellar (and gas phase) metallicities agree well with the observed mass-metallicity relation both as a  \nfunction of stellar mass at z = 0 and as a function of redshift (Ma et al. 2016a); to the extent that galaxies differ in metallicity between FIRE-1 and FIRE-2, they primarily move along the stellar mass-metallicity relationship. In both sets of simulations, galaxies drive strong winds, with higher mass-loading in low mass galaxies, sufficient to place galaxies on the observed relationship between stellar mass and halo mass (Hopkins et al. 2014; Muratov et al. 2015; Hayward & Hopkins 2017). Initial examination of our FIRE2 runs shows these winds produce galactic outflow rates and covering factors of neutral hydrogen similar to our FIRE-1 simulations, similar to observations at a wide range of galaxy masses (FaucherGiguere et al. 2015; Faucher-Giguère et al. 2016), although this will be studied in more detail in the future.  \nThere are some modest quantitative differences between FIRE-2 and FIRE-1; most of this manuscript will explore the origin of these differences, but they appear to primarily owe to the change in the hydrodynamic solver. On average, dwarfs are slightly lowermass in FIRE-2; this owes to both hydrodynamics and the fact that (unlike in FIRE-1) we do not artificially 'ignore' ionizing photons when they pass outside numerical domains (hence they can still heat gas). The difference in Fig. 5 can be as large as a factor of ∼ 3 for a single galaxy, but this is largely stochastic (since the star formation histories are dominated by a few bursts, a small perturbation to the formation history or feedback strength can lead to factor ∼ 2 changes in mass) - Fig. 7 makes this clear, as the systematic offset between FIRE-1 and FIRE-2 appears to be a factor < 2, well within the systematic uncertainties in the M ∗ -M halo relation at M halo ≤ 10 11 M glyph[circledot] . We have also examined this relation at z =( 0 . 5 , 1 , 2 , 4 , 6 ) and find similar agreement, so we refer to Hopkins et al. (2014) for an extensive analysis. Metallicities at the same mass are slightly higher, by < 30%; this owes primarily to updated yield tables; this is far smaller than the factor of ∼ 2 -5 systematic uncertainty in the calibration of observed galaxy metallicities. For MW-mass galaxies, the stellar masses are slightly higher in FIRE2, and the bulges slightly more concentrated (the rotation curves have stronger peaks, by a modest amount). We will show that this is a direct consequence of the hydrodynamic treatment, but is also sensitive to simulation resolution and the SNe coupling algorithm.", '4 RESOLUTION IN THE FIRE-2 SIMULATIONS': "We now discuss our mass, spatial, and time resolution. For each, we will present a series of tests, and summarize these with a set of resolution criteria.  \nThroughout this paper, when we refer to 'convergence' of some property, we more precisely mean to test whether the property is strongly, systematically sensitive to resolution (at our highest resolution). Because even pure N -body integration (let alone galaxy formation, with explicitly stochastic effects such as SNe included) is a fundamentally chaotic problem, and has no known exact solution, we of course cannot define convergence in the more formal sense (of e.g. some error norm relative to said solution). Some properties in nature (e.g. halo mass functions, or turbulence in the ISM) extend down to scales vastly smaller than any simulation could resolve - for these, convergence must be defined 'down to' some minimum resolvable scale (e.g. the mass function for objects larger than some minimum number of particles). This also means that 'convergence' in one quantity should not be taken to necessarily imply convergence in another.", '4.1 Mass Resolution': "In Lagrangian or N -body methods such as ours, there is a well defined mass-resolution given by the particle mass mi in the highresolution Lagrangian region. 11  \nTo maintain a well-defined mass resolution scale and minimize N -body integration errors, single gas particles are converted into single star particles with the same masses (rather than, for example, spawning star particles with much smaller/larger masses). However, mass loss from stars to gas in O-star/AGB winds and SNe means that particles will not have perfectly equal masses, so to prevent pathological behavior in very rare circumstances (if e.g. a single gas particle sees many SNe that increase its mass) we use the standard particle splitting and merging routine from Hopkins (2015) to ensure no particle ever deviates by more than a factor of 3 from the median particle mass. This affects only a tiny number of particles (one in ∼ 10 6 ). Averaged over our runs at z = 0, ∼ 90% (99%, 99 . 99%) of all gas particles are within < 0 . 005 (0 . 1, 0 . 2) dex of the median particle mass. We will therefore refer only to the median baryonic particle mass mi in this paper.  \nDark matter particles are always more massive by the universal baryon fraction, m dm ≈ 5 mi ; this ensures halos and baryonic galaxies are comparably resolved at initial collapse. Of course, since many galaxies retain only a small fraction of their baryons, and dark matter does not cluster on small scales, dark matter structures tend to be vastly larger than baryonic structures and so are far better mass-resolved.  \nWedefine, for convenience, the baryonic particle mass in units of 1000 M glyph[circledot] , and note that this specifies both baryonic and dark matter particle masses:  \nmi , 1000 ≡ m baryonic i 1000 M glyph[circledot] (1)  \nm DM = Ω m -Ω b Ω b m baryonic i ≈ 5000 M glyph[circledot] mi , 1000 (2)  \n4.1.1 Requirements for 'Resolved' Self-Gravitating Structures  \nThe mass resolution 'required' to accurately model different phenomena depends, of course, on the question being asked and desired level of accuracy. A wide range of studies have shown that structures with masses of ∼ 5 -100 times the element mass are 'believable' in the sense that they exist and can be identified as self-gravitating in higher-resolution re-simulations (Klypin et al. 1999; Springel et al. 2001; Kravtsov et al. 2004; Nurmi et al. 2006; Boylan-Kolchin et al. 2009; Wetzel & White 2010). This corresponds to halos, stellar galaxies, or gas clouds of mass:  \nM min halo ∼ 0 . 6 × 10 5 M glyph[circledot] N 10 mi , 1000 (3)  \nM min ∗ = M min cloud ∼ 1 × 10 4 M glyph[circledot] N 10 mi , 1000 (4)  \nwhere mi , 1000 ≡ mi / 1000 M glyph[circledot] is the baryonic particle mass, to which we reference all quantities, and N 10 = N desired / 10 reflects the 'desired' number of particles ( N desired).  \nWe will show below for DM halos and subhalos, stellar galaxies, gas clumps and GMCs within galaxies, only a few particles are sufficient for robust prediction of masses and mass functions.  \n11 The low-resolution regions of the box, outside several virial radii of the main galaxies, are populated by lower-resolution collisionless particles, stepping up in a powers-of-eight hierarchy. This is sufficient for resolution of long-range tidal forces from these regions but we do not consider any halo contaminated ( > 1% by mass) by these particles to be 'resolved.'  \nTable 3. Illustrative Examples of 'Resolution' Scales in a MW-Mass Halo ( m12i )  \nSeveral 'spatial resolution' and 'time resolution' properties of the simulations discussed in § 4. We focus on our resolution study of the m12i system, but results for other MW-mass galaxies are nearly identical, and results for dwarf galaxies are qualitatively similar at similar mass resolution. Because our simulations are Lagrangian, only mass resolution is truly 'fixed.' Spatial (both hydrodynamic and gravitational, which are the same always) and time resolution are both adaptive, and in principle can reach arbitrarily small values, but in practice reach minimum values based on the densest mass-resolved structures in the simulation.  \nOf course, the presence of even a few baryonic particles might require a quite massive halo. The observed relationship between galaxy stellar mass and halo mass implies that galaxies with M ∗ ∼ 10 4 M glyph[circledot] typically live in halos of mass M halo ∼ 10 9 M glyph[circledot] ; using the M ∗ -M halo relation at low masses from Moster et al. (2013), we estimate the minimum halo mass with a 'believable' baryonic relic is ∼ 10 9 M glyph[circledot] at mi , 1000 ∼ 1. This is actually quite well-resolved in dark matter, with glyph[greaterorsimilar] 10 5 particles. This also means that glyph[greaterorsimilar] 10 5 baryonic particles have participated in the formation history of the halo and have 'cycled through it' (assuming something like the Universal baryon fraction is associated with the halo). It is just that the star formation efficiency is so low that the residual mass in stars is small. There is no question, then, that such objects are 'real' in the simulations and, if they have such small stellar masses, that  \nfeedback had a real effect (it had to prevent > 10 5 particles from cooling and forming stars!).  \nMoreover, because dark matter dominates the gravitational forces in these small galaxies, the orbital dynamics of the surviving stars (determined by the dark matter potential, not the negligible self-gravity) can be believed so long as the dark-matter structure is resolved (i.e. the stars are just tracer particles). This is easily satisfied in halos with M halo glyph[greaterorsimilar] 10 9 M glyph[circledot] , independent of the number of stars (given a realistic M ∗ -M halo relation).  \nWe will show below that even the internal dynamics, evolution over a Hubble time, and mass profile of DM structures, are robust to resolution down to radii enclosing just ∼ 200 elements. This is much more demanding than 'existence' of halos or gaseous structures (which requires just a few elements), as expected.  \nFigure 8. Resolution study, as in Fig. 4. Each column is one galaxy, and each line shows a different mass resolution ( mi , 1000) up to our production resolution. Top: Star formation rate versus cosmic time. In dwarfs ( m10v, q ), low resolution leads to artificially large 'burstiness': a few bursts dominate the history and eject the remaining baryons. At higher resolution, star formation is still 'bursty', but it becomes robust to changes in resolution. In MW-mass galaxies, higher resolution shifts star formation to slighly later times by more efficiently regulating the low-mass projenigtor at high redshifts. Second from Top: Stellar mass versus scale factor. Note that the mi , 1000 = 16 (130) simulations of m10v have only 10 (2) star particles in the main halo, and yet they are within a factor ∼ 3 of the highest-resolution stellar mass; by the time ∼ 100 particles are in the main galaxy, the mass is stable to within ∼ 20%. Middle: Average stellar metallicity versus scale factor. This converges more slowly than stellar mass; with < 100 star particles, dwarf galaxies show artificially suppressed metallicity, because low resolution under-samples the enrichment history and leads to artificially bursty SFH that blows out metals completely. Massive galaxies show smaller differences that match differences in their SFHs. Second from Bottom: Mass density profile at z = 0: results are robust down to radii enclosing ∼ 100 -200 particles of the 'type' of interest (details below). Bottom: Circular velocity profile at z = 0. For dwarfs, this is dark-matter dominated and therefore underresolved only at our lowest resolution (where the DM 'convergence radius' discussed below is > kpc). For the MW-mass galaxy, the modest difference in the SFH at z ∼ 0 . 5 -2 in m12i translates to a factor ∼ 2 difference in the central bulge mass, which leads to a more strongly peaked central Vc at low resolution. m12f (not shown) shows nearly identical behavior to m12i ; m12m shows better agreement at all resolution levels.  \n<!-- image -->", '4.1.2 Requirements for Well-Behaved Internal Hydrodynamics': 'In Hopkins (2015), we show that, for the MFM method here, a few hundred resolution elements are sufficient to capture the orbital dynamics of thin Keplerian disks for ∼ 10 -100 orbital times (for a galactic disk, roughly equivalent to a Hubble time), the shock structure of strong blastwaves and/or implosions (the Sedov and Noh tests), and all the correct qualitative behaviors of self-gravitating polytropic sphere collapse (the Evrard test) and cosmological structure formation (the Zeldovich test), all to within a factor much better than ∼ 2 of the exact solution. This means that a single star particle generating SNe in a halo and blowing out ∼ 100 -1000 surviving gas particles (total baryonic mass glyph[greaterorsimilar] 10 6 M glyph[circledot] mi , 1000 or, at the Universal baryon fraction, halo mass glyph[greaterorsimilar] 10 7 M glyph[circledot] mi , 1000) is  \nat least hydrodynamically well-behaved, if not converged. 12 This is consistent, roughly, with the resolution dependence of properties such as the enrichment history presented below, which require similar particle numbers in the relic baryonic galaxy to resolve the entrainment/recycling of SNe ejecta and its re-incorporation into subsequent generations of star formation.', '4.1.3 Requirements for Capturing ISM Properties, Star Formation, and Stellar Feedback': "The most demanding mass resolution requirements come from the physics of the ISM and CGM.  \nFigure 9. Mock images, as in Fig. 3, showing the effects of resolution on the morphology of the simulated galaxies at z = 0, for a subset of the resolution series from Fig. 8. Left: Dwarf galaxy ( m10q ): because dwarf morphologies are disordered, their qualitative morphology does not depend sensitively on resolution (this is true for all our dwarfs). Center: MW-mass galaxy m12i : In m12i we see a trend towards a more extended thin disk component dominating as we increase resolution. The effective radius also increases, but much less dramatically (just by ∼ 40% from mi , 1000 = 56 to mi , 1000 = 7 . 0); the bulge and central ∼ 2 -3 kpc remain similar in each case, but the extended, gas-rich disk is much more prominent. At the lowest resolution ( mi , 1000 = 450), we see little disklike structure at all (although there is a rotating gas+stellar disk ∼ 1kpc in size). At intermediate resolution, a clear disk with ∼ 5kpc radius appears - there is also an extended, smooth-light component out to ∼ 10kpc. At our highest resolution, this extended component exhibits spiral structure and the gas+young stellar disk extends to > 10kpc from the galaxy center (albeit at low surface brightness: all figures here use an 8-magnitude stretch). Right: MW-mass galaxy m12f : In this slightly more-massive galaxy, an extended disk is present at all resolution levels (the trends seen in m12i are still present, but much weaker). m12m (not shown) similarly shows a disk at every resolution level.  \n<!-- image -->  \n(i) The Toomre Mass & ISM Structure: Many previous studies have shown that reliably capturing star formation (both rates and spatial distribution) requires the ability to resolve at least the existence and self-gravity of the largest self-gravitating gas structures (i.e. fragmentation) in a galactic disk (see Saitoh et al. 2008; Hopkins et al. 2011, 2012c,b, 2014; Kim & Ostriker 2017). These studies showed that once this criterion is met, the SFR predicted in the simulations becomes independent of the numerical star formation model (see § 8), because it is regulated by feedback. This means we do not strictly need to resolve the internal dynamics of such clouds to capture galaxy-scale dynamics. And although there will always be clouds and sub-structure below our resolution limits, both observations and simulations have repeatedly shown that almost all star formation occurs in the largest GMC complexes in galaxies (see Evans 1999; Rosolowsky et al. 2003; Blitz & Rosolowsky 2005; Bolatto et al. 2008; Elmegreen 2007; Tasker & Tan 2009; Feldmann & Gnedin 2011; Murray 2011; Harper-Clark & Murray 2011; Hopkins 2012a). This mass M max GMC is set by the Toomre mass, M Toomre ∼ f 2 gas M gas where f gas is the gas fraction inside the relevant effective radius of the baryonic galaxy. So, if we require M min cloud glyph[lessorsimilar] M Toomre, we obtain the desired resolution criterion  \nmi , 1000 glyph[lessorsimilar] m Toomre i , 1000 where  \nm Toomre i , 1000 glyph[lessorsimilar] 100 N 10 M Toomre 10 6 M glyph[circledot] ∼ 100 ( f 2 gas M gas 10 6 M glyph[circledot] N 10 ) (5)  \nThis defined the target resolution of the original FIRE-1 massivegalaxy simulations (10 12 M glyph[circledot] halos).  \n(ii) SNe Cooling Radii: A still more demanding criterion is set by SNe physics. As discussed in detail in Paper II (as well as Cioffi et al. 1988; Thornton et al. 1998; Martizzi et al. 2015; Walch & Naab 2015), a single SN remnant radiates its thermal energy rapidly upon reaching a nearly-invariant 'swept up' mass M cool ∼ 3000 M glyph[circledot] ; if a number N SNe occur before the cooling time expires, the 'cooling mass' simply scales M cool ∝ N SNe. If SN energy is injected in a kernel-weighted fashion over N NGB elements as we do here, then almost all of the energy goes into the nearest N 1 / 3 NGB neighbors; requiring they have a total mass < M cool in turn requires:  \nm SNe i , 1000 glyph[lessorsimilar] 0 . 9 ( N NGB / 32 ) -1 / 3 N SNe (6)  \nThis set the target mass resolution for the original FIRE-1 dwarf galaxies ( < 10 10 M glyph[circledot] halos). As we show in Paper II, unless mi , 1000 glyph[lessmuch] m SNe i , 1000 , it is necessary to properly account for the con-  \nFigure 10. Star formation rate versus lookback time, in our MW-mass m12i simulations, re-started at late times. In each case, we re-start using the snapshot at z = 0 . 06 from our run with mi , 1000 = 56 (in Fig. 8) as our initial condition, so the late-time ICs of the different realizations are identical. We then run from z ≈ 0 . 06 -0 ( ∼ 700Myr physical time), to study how the SFR varies in a massive galaxy given the same initial galaxy properties (that is, factoring out how variations affect the earlier phases when the progenitor was low-mass). Here, we use particle splitting/merging before running to vary the mass resolution. We see that for the same ICs, the SFR is almost completely insensitive to resolution; variations in the SFR with resolution in Fig. 8 are dominated by (1) early stages when the progenitor galaxy was much lower mass, hence much less well resolved, and (2) less prominently, weak resolution dependence of wind mixing versus escape from the outer halo, which changes the inflow rate back into the galaxy at later times.  \n<!-- image -->  \nversion of energy into momentum in the unresolved Sedov-Taylor (S-T) phases. Failure to do so will significantly under-estimate the effects of feedback. 13  \n(iii) Dwarf Galaxy 'Burstiness': In previous papers we have shown that for small dwarf galaxies, star formation is robustly 'bursty'; 14 however, at low resolution, the 'burstiness' is artificially enhanced by numerical effects (see e.g. Sparre et al. 2017). This owes to the fact that stars form in units of star particles. At sufficiently low resolution, a single star particle implies a massive, co-eval star cluster, therefore a large number of approximately coeval SNe. Although we know SNe are clustered in reality, this artificial numerical clustering could easily 'overshoot' reality, producing too many synchronized SNe in the same location, which in turn leads to an over-large super-bubble that could heat the entire galaxy gas supply to super-virial temperatures. If we consider the 'unit' of star formation (minimum resolved cloud) to be, say  \n13 Note that Eq. 6 does not mean we cannot resolve hot gas and/or overlapping SNe bubbles at lower resolution. If any mass of stars m 0 forms and the corresponding N SNe ∼ m 0 / 100 M glyph[circledot] go off in an overlapping resolution element within a cooling time, then they can heat a mass ∼ 60 m 0 to > 10 6 K. So because star formation is clustered, we can still resolve super-bubbles and galactic chimneys at relatively modest resolution. However, resolving the full hot gas content, venting, momentum contribution from confined blastwaves, and early evolution of SNe explosions requires a criterion like Eq. 6.  \n14 For the sake of quantitative comparison, we define a specific measure of 'burstiness' as the standard deviation in the quantity log ( ˙ M ∗ (∆ t 1 ) / ˙ M ∗ (∆ t 0 )) where ˙ M ∗ (∆ t ) is the SFR averaged over a time interval ∆ t , and we compare a short interval ∆ t 1 = 10Myr and longer interval ∆ t 0 = Gyr.  \n∼ 10 particles (since a single lone gas particle cannot be at much higher density than its neighbors), then taking the IMF-average N SNe ∼ mi / 100 M glyph[circledot] per particle, the entire baryonic mass in the galaxy can be heated to T glyph[greatermuch] 10 6 K ('super-bubble' temperatures) if that baryonic mass is below ∼ 600 mi . So resolving 'venting' of individual super-bubbles (much smaller than the entire galaxy) requires glyph[greaterorsimilar] 1000 gas particles in the galaxy. Assuming further a gas fraction of ∼ 1 during the gas-rich phases of star formation and/or typical dwarf galaxies, and using the stellar mass-halo mass relation above, this implies a particle mass target:  \nm burstiness i , 1000 glyph[lessorsimilar] 5 ( M halo 10 10 M glyph[circledot] ) (7)  \nEquivalently, we can say that for a simulation with fixed mass resolution, once halos exceed a mass M halo glyph[greaterorsimilar] 2 × 10 9 M glyph[circledot] mi , 1000 the numerical 'excess clustering/burstiness' is not important.  \n(iv) The Jeans Mass: A common mis-conception is that one needs to resolve the Jeans mass in order to capture basic fragmentation physics. This is incorrect, for two important reasons. First, the warm/cold ISM is super-sonically turbulent, so the 'fragmentation cascade' (hierarchical structure of fragmentation into GMCs, clumps, cores, etc.) is determined by the turbulent Jeans mass, not the thermal Jeans mass, at least down to the sonic scale ( R glyph[lessorsimilar] R sonic ∼ 0 . 1pc and mass scale ∼ 1 M glyph[circledot] ). We discuss this below, but the turbulent Jeans mass of a super-sonically turbulent GMC with virial parameter α ∼ 1 is of order the GMC mass itself ( ∼ 10 7 M glyph[circledot] ), whereas the thermal Jeans mass is ∼ 0 . 1 M glyph[circledot] . Second, in a homogeneous medium, the Jeans mass defines the smallest scales of the fragmentation cascade - by definition, all larger scales are also unstable. Therefore, failure to resolve the Jeans scale simply means all resolved scales fragment as they should - the fragmentation cascade is just truncated at the resolution scale, instead of the (smaller) Jeans scale. This is analogous to 'resolving turbulence' in the ISM; the actual Kolmogorov ('termination') scale of the turbulent cascade ( ∼ au) is far smaller than achievable resolution, but that does not mean turbulence cannot be captured (it simply limits the dynamic range of the cascade that can be followed). In both cases (for the same physical reason, in fact; see Hopkins 2012a), the power in the cascade is dominated by the largest-scale structures, so integrated quantities converge very quickly with increasing resolution (Padoan et al. 2012; Federrath & Klessen 2012; Hopkins 2013b; Guszejnov & Hopkins 2015, 2016). We show this explicitly below - the mass function of dense, cold clouds is well behaved regardless of resolution, it simply extends to smaller and smaller sub-structures as we increase the resolution. Thus, the full cascade is not followed, but this is not problematic, provided that the goal of our simulation is not to resolve individual brown dwarf or star formation (the actual 'end point' of the fragmentation cascade in the ISM). This is where the sub-grid model for star formation enters - it is, explicitly, a sub-grid model for the un-resolved fragmentation cascade from a resolved, self-gravitating cold gas clump/cloud all the way down to an aggregate of individual stars with some IMF. The Toomre scale, on the other hand, is important to resolve, because it is the largest unstable scale in a disk in other words, failure to resolve the Toomre scale means that no fragmentation will occur (when it physically should).", '4.1.4 Mass Resolution Tests': "We now explore in detail how mass resolution can alter our conclusions. Figs. 8-14 present a series of explicit mass resolution tests of our full-physics FIRE-2 simulations. Table 3 gives some typi-  \nFigure 11. Resolution studies in DM-only simulations of the dwarf and MW-mass halos from Fig. 8. We show the radially averaged DM mass profile at z = 0 around the main halo, varying the mass resolution. Inset shows the difference plot ( ∆ log ρ in dex) relative to the 'reference' highest-resolution run available. For each factor of 8 in mass resolution, we change the DM force softening by a factor of 2. With increasing mass resolution we converge more and more closely to an NFW-like profile at small radii (thin green line). At radii larger than the Power et al. (2003) radius, r 0 . 6. where the relaxation time t relax ≈ 0 . 6 t circ ( R 200 ) , which encloses approximately ≈ 2200 particles (dotted vertical lines), the agreement is near-perfect. However, even at radius, r 0 . 06, where the t relax ≈ 0 . 06 t circ ( R 200 ) , which encloses just ≈ 200 particles (dashed vertical lines), the agreement is still quite good: densities are under-estimated by at most ∼ 0 . 05 -0 . 15 dex. We define the latter as our DM 'convergence radius' throughout this paper.  \n<!-- image -->  \ncal values for mass, spatial, and time resolution for m12i runs at different resolution.  \nFig. 8 shows our full-physics cosmological simulations at varying mass resolution. We compare two low-mass dwarfs and two MW-mass galaxies, to bracket the range of behavior. Fig. 9 compares the visual morphologies of the galaxies.  \nFor the dwarfs, the total stellar mass (and circular velocity profile) becomes robust to better than a factor < 2, with just glyph[greaterorsimilar] 2 -16 star particles in the z = 0 galaxy. This is because the total stellar mass is set by an integral balance between feedback energetics and gravity binding the baryons. Also, given the low particle masses, even these low-resolution dwarf runs easily satisfy our Toomre-mass criterion (Eq. 5). However, at such low resolution ( mi , 1000 = 16 -130) they do not satisfy our 'numerical burstiness' or 'SNe cooling' criteria (Eq. 7 and 6, respectively). As a consequence, the SFHs are visibly more 'bursty' (dominated by just a couple large bursts), and the metallicities are systematically under-predicted (compared to our high-resolution runs). The latter occurs because the single bursts blow out nearly all the baryons (and metals) from the galaxy: they fail to capture partial entrainment/mixing/incomplete blowout that would keep the metals in the  \ngalaxy. The metallicity and burstiness of the SFH appear to robust to a factor ∼ 2 ( ∼ 10%) when the number of stars reaches > 100 ( > 500) at mi , 1000 = 2, as expected from our criterion in Eq. 7. More subtle properties, such as the internal SF structure of bursts and hot gas properties of the galaxy, and escape fraction of ionizing photons (Ma et al. 2015, 2016b), require still higher resolution, hence our highest-resolution runs, with mi , 1000 = 0 . 25, satisfying Eq. 6. These results are also supported by the analysis of dwarf satellite galaxies of our massive, high-resolution m12i halo in Wetzel et al. (2016); the convergence in the stellar mass function of dwarfs appears good down to ∼ 5 -10 star particles per galaxy; but the metallicities of the satellites with glyph[lessorsimilar] 100 star particles are suppressed relative to our higher-resolution isolated-dwarf simulations and observations (compare Ma et al. 2016a). Not surprisingly, because the gross morphology of the dwarfs is irregular, it is essentially resolution-independent. Galaxy sizes (effective radii) are also nearly independent of resolution (because of the more 'bursty' star formation, the sizes tend to be ∼ 10 -20% larger at the lowest resolutions here). In Paper II, we show the star formation histories of dwarf galaxy simulations reaching ∼ 30 M glyph[circledot] resolution (which will  \nFigure 12. Resolution studies of DM mass profiles in the dwarf and MW-mass halos from Figs. 8-11, in simulations with baryons (solid) compared to DMonly simulations (dotted). We re-normalize the simulations with baryons by Ω m / (Ω m -Ω b ) so that if the baryons traced the DM perfectly, the DM-only and baryonic simulations would agree exactly. For clarity, we label the 'convergence radius' ( r 0 . 06, enclosing ≈ 200 DM particles; vertical dashed lines) for just the simulation with baryons. In m10v , the galaxy is sufficiently low mass ( M ∗ ∼ 10 5 M glyph[circledot] ) that it is DM dominated and has little or no core, so DM and baryonic simulations agree well. We also see good agreement down to radii ∼ 30pc enclosing just ∼ 10 particles. In m10q , the more massive dwarf is still DM-dominated but forms a core in the baryonic runs (suppression of ρ DM to ∼ 600pc); the core is actually more robust to resolution than the 'cusps' in the DM-only runs, down to similar radii ∼ 30pc. In MW-mass runs, the baryons dominate the central mass and cause some contraction of the DM (although less than would be expected from pure adiabatic contraction given the galaxy mass). Differences in the central DM profile are, in these galaxies, dominated by differences in the baryonic mass (high-resolution runs can be less dense, owing to slightly smaller stellar masses), not by traditional DM resolution considerations such as N -body relaxation.  \n<!-- image -->  \nbe studied in more detail in Wheeler et al., in prep), and show that they agree well with the ∼ 250 M glyph[circledot] -resolution runs here.  \nFor the massive, MW-mass system, even the lowest-resolution ( mi 1000 = 450) run has glyph[greatermuch] 10 5 particles in the baryonic galaxy and easily satisfies the Toomre and 'burstiness' criteria (although its progenitor dwarf galaxies at high redshifts may not). As such the SFH and metallicity are more robust to resolution (compared to the dwarf runs). However we do see higher resolution systematically shifts the SFH from early to later times, as higher resolution allows better resolution of two key physics. First, the generation of winds via resolution of the hot gas 'channels' and their escape ('venting') from a multi-phase halo (our 'SNe' criterion; for explicit studies see Hopkins et al. 2012b, 2013c; Muratov et al. 2015; Martizzi et al. 2015; Fielding et al. 2017). At low resolution, hot gas is necessarily 'dragged' by a large mass of implicitly-coupled gas, because mass is locked into massive particles: at our lowest resolution, for example, a single SN cannot affect less mass than ∼ 10 6 M glyph[circledot] - this means much lower 'launch velocities' and tem-  \nratures, unless a huge number of SNe explode simultaneously (see Paper II). Second, high resolution allows better resolution of the 'burstiness' and SFHs within the smaller progenitor galaxies of the more massive z = 0 MW-mass system (note that the largest differences appear at early times, when the galaxy is a progenitor dwarf), which reduces their stellar and gas masses, leading to more gas expelled to large enough radii where its recycling times are long. The gas is still re-incorporated eventually, evident in the similar late-time masses and SFRs, but appears to re-accrete later.  \nAs a result of these effects, at higher resolution the final z = 0 MW-mass galaxy is slightly lower-mass, but more importantly, because it shifts star formation and re-accretion of recycled material to later times, that material carries larger angular momentum, and the galaxy is less compact. For m12i , the effective radius in our high-resolution mi , 1000 = 7 run is ∼ 1 . 4 times larger than the mi , 1000 = 450 run, which translates to a factor ∼ 1 . 8 lower mass inside < 5kpc; this in turn lowers the peak in the rotation curve from ∼ 370kms -1 to ∼ 270kms -1 . Thus the galaxy rotation curve is  \nFigure 13. Resolution study of ISM structure, in the runs from Fig. 10 (re-running the MW-mass m12i from z = 0 . 07 -0 with different mass resolution). Top: Time-averaged mass function (MF) of 'GMCs': dense, cold gas clouds (identified with a friends-of-friends algorithm; see § 4.1.4). We plot the total mass in clouds, per logarithmic interval in cloud mass ( M GMC dN GMC / d log M GMC), as a function of mass ( M GMC). All identified structures as small as 3 gas particles are shown. We compare the observed GMC MF in the MW (Rice et al. 2016; inside/outside the solar circle as solid/dotted), normalized to the same total mass. The MF has power-law slope dN / dM ∝ M -( 1 . 6 -1 . 8 ) and exponential cutoff around the Toomre mass ( ∼ 10 7 M glyph[circledot] ), similar to both observations and analytic predictions from turbulent fragmentation theories. Bottom: Linewidth (internal velocity dispersion) vs. size relation (for the same clouds; median and 5 -95% interval in thick/thin lines), compared to observations (labeled). Down to ∼ 3 -5 particles per cloud, the predictions (both MF and linewidth-size) agree well; higher resolution simply samples smaller clouds. Most of the GMC mass is around the Toomre mass, so the total cloud mass is within ∼ 30% in all runs shown. Only the lowest-resolution run cannot resolve the peak, biasing the mean M GMC higher by a factor ∼ 2.  \n<!-- image -->  \nnoticeably 'less bulgy' at high resolution. The differences in m12i are more dramatically evident in the low surface-brightness outer disk morphology, which goes from being entirely absent at very low resolution to quite prominent at high resolution. Interestingly, however, in m12m there is much weaker dependence on resolution (the effective radius changes by < 15%). Also in m12f , although the change in effective radius and the rotation curve with resolution is similar to m12i , the visual morphology changes much less dramatically - even at mi , 1000 = 450 there is still a prominent, extended thin disk. This may owe to the fact that m12m and m12f have a some-  \n3  \n4  \nSatellite Distance  \nr  \n[  \nkpc  \n]  \nFigure 14. Mass function ( top ), circular velocity distribution ( middle ), and spatial distribution ( bottom ) of satellites (DM subhalos) within the primary halo in DM-only runs of our m12i simulation (at z = 0), at varied mass resolution (as labeled). Top: Number of subhalos versus bound subhalo mass. Vertical lines show 32 DM particles; agreement is good down to structures containing ∼ 10 -30 particles. Middle: Number of subhalos versus maximum circular velocity. Resolution dependence is similar to the mass function (deviations occur only below ∼ 10 -15kms -1 , corresponding to the un-resolved subhalos above). Bottom: Distribution of subhalos in radial distance from the center of the primary halo. Up to shot noise in the exact position of the individual subhalos in their orbits, the distributions agree well. As the DM is approximately self-similar, we find qualitatively identical results comparing our other MW-mass ( m12f ) or dwarf mass ( m10q , m10v ) simulations.  \n<!-- image -->  \nwhat larger mass, but likely also owes to the specific fact that m12i has a series of mergers around z ∼ 1, which launch strong fountains and change the angular momentum of the gas that will form its disk (whereas m12f and m12m grow more smoothly). So it appears m12i is simply more sensitive to resolution effects. Clearly, it is important to push to even higher resolution (in progress), to test whether or not this is actually converged.  \nIf the dependence on resolution in our MW-mass runs owes to  \n3  \na combination of (1) better-resolving the progenitor (dwarf) galaxies, and (2) resolving channels by which hot gas outflows can escape, we should see similar SFRs independent of resolution at late times if we start from the same ICs . In other words, it is useful to test whether indeed the difference with resolution comes from progenitor influence on the late-time galaxy, or whether it is present for fixed conditions in a massive MW-mass galaxy. We therefore consider a series of simulations where we use our m12i simulation with mi , 1000 = 56, re-started at z ≈ 0 . 06 from a snapshot of the simulation in Fig. 8, but using our particle splitting/merging routine to first split/merge the ICs until they are re-sampled with a desired target resolution. Fig. 10 shows that when we do this, the SFR is nearly identical over ∼ 2 . 5 dex in mass resolution (small deviations at the highest resolution owe mostly to some artifacts of our very aggressive particle splitting/up-sampling routine applied for this specific test).  \nThis is consistent with our argument above, that we only need to marginally resolve the Toomre scale to achieve a robust prediction for the SFR, given a specific gas disk initial condition. Quantities such as the Kennicutt-Schmidt relation are thus extremely robust to mass resolution. However, non-linear, long-timescale effects in cosmological simulations (e.g. recycling) shift the SFR by changing the supply or loss rate of gas in the CGM .  \nTo better understand how much of the resolution dependence owes to purely gravitational physics, Figs. 11-12 consider how the dark matter mass profiles of the galaxies change with mass resolution. First we consider DM-only simulations (i.e. gravity-only simulations), in Fig. 11. We see near-ideal convergence in the mass profiles of the DM halos (as well as their substructure mass functions, shown below, and halo formation/growth histories). As expected, the z = 0 DM mass profiles are well-fit by a Navarro et al. (1996) (NFW)-like profile, down to some minimum scale where numerical effects flatten the profile. Because of the Lagrangian nature of our code, improving mass resolution also improves the effective spatial resolution/force resolution; in fact, we will show below (consistent with many previous studies) that the nominal spatial force softening is generally much less important than mass resolution. Power et al. (2003) argue that the central 'flattening' in DM profiles is dominated by N -body relaxation, and that robust results should be obtained outside a radius where the N -body relaxation time t relax ∼ 0 . 6 t 0 (where t 0 ≡ t circ ( R 200 ) = 2 π R 200 / V 200) is comparable to the Hubble time. Because N -body relaxation depends most strongly on N , this is effectively a requirement on the number of particles more exactly ( 5 / √ 8 )( N [ < r ] / ln N [ < r ] )( ρ crit / ¯ ρ [ < r ] ) 1 / 2 < 0 . 6, satisfied for the densities in Fig. 11 when N glyph[greaterorsimilar] 2200. Not surprisingly, we see almost excellent agreement at these radii. In fact, as others have shown, because of more accurate integration criteria, shorter timesteps, and a smoother spline for gravitational softening (which reduces the N -body relaxation time below the equation for strict point masses assumed by the argument in Power et al. 2003), if we are willing to tolerate slightly larger errors, within a factor ∼ 1 . 3 (0 . 1dex) of the converged (NFW) solution we find excellent agreement independent of resolution, down to radii containing just N glyph[greaterorsimilar] 200 particles.  \nFig. 12 compares the DM-only result to the full baryonic physics runs from Fig. 8. Following standard practice, we correct the profiles from baryonic runs by the cosmic mean baryon fraction, so that if the baryons behaved identically to the DM, the curves would lie exactly on top of one another. With or without baryons, in both our dwarf and MW-mass runs, agreement is good to glyph[lessorsimilar] 0 . 15dex outside the radii enclosing ∼ 200 particles. In m10v , the z = 0 galaxy is sufficiently low-mass that baryons have a neg-  \ngible effect on the DM profile (see Chan et al. 2015), so the two sets of runs track each other closely. In m10q , the galaxy starts to reach masses where stellar feedback generates a 'core': even though small (factor < 2) differences in the stellar mass have large effects on core creation at these masses (see Oñorbe et al. 2015), the resolution dependence is similar in baryonic and DM runs. In the MW-mass runs, we see that runs with baryons have higher central densities than DM-only runs, because the large galaxy baryonic mass has caused some halo contraction. This necessarily leads to better 'convergence' in an N -body sense, but it means that the DM profile is more sensitive to changes in the galaxy stellar mass (here, the highest-resolution run has a lower M ∗ and correspondingly lower central DM density). In all cases, our highest-resolution baryonic simulations reach an approximate 'convergence radii' of ∼ 50 -100 pc in the dark matter . We of course achieve much higher effective spatial resolution in the baryons dense enough to form stars.  \nNote that there is no analogous Power et al. (2003) criterion for baryons. However, Fig. 13 and various self-gravitating baryonic collapse tests in Hopkins (2015) show that self-gravitating structure can be captured across just a couple inter-particle separations. Larger numbers of particles are needed for convergence in the mass profile in dark matter because the N -body orbits must be integrated for a Hubble time. In the gas, however, self-gravitating sub-structures (e.g. GMCs) survive for only a few dynamical times. Thus the appropriate Power et al. (2003)-equivalent criterion is t relax glyph[greaterorsimilar] couple × t dyn, which is easily satisfied even for point masses (un-softened gravity) whenever an object is comparable in mass (or larger) than a kernel. With adaptive softening for gas (our default choice), this is by definition true even for 'structures' of ∼ 2 particles.  \nFig. 13 specifically considers the resolution dependence of structure within the ISM; we use the same simulations from Fig. 10 of the MW-mass system at z ∼ 0, but study the mass function of dense, cold gas structures within the ISM. 15 A more detailed discussion of the cloud properties in these simulations is presented in Guszejnov et al. (2017); for our purposes here we simply desire a proxy for the GMC mass function to understand the resolutiondependence of the simulations. In each case, the shape of the mass function is, as expected, a power-law with slope dN / dM ∝ M -( 1 . 6 -1 . 8 ) (very similar to observed; see Blitz & Rosolowsky 2005; Rice et al. 2016), and turnover at a maximum mass about the Toomre mass M Toomre defined above (also as expected from observations and analytic theory; Murray 2011; Hopkins 2012a). The lower limit is purely a resolution effect: we only keep structures with ≥ 3 particles. We also compare the internal properties of the clouds, specifically the linewidth-size relation (the onedimensional velocity dispersion of gas within clouds, versus their projected mass-weighted rms radius), to observations of nearby galaxies ((Bolatto et al. 2008; Fukui et al. 2008; Heyer et al. 2009; Muraoka et al. 2009; Roman-Duval et al. 2010; Colombo et al. 2014; Heyer & Dame 2015; Tosaki et al. 2017); note our definition of R cloud is equivalent to their σ r ). The generic power-law scaling  \n15 For simplicity, we use a friends-of-friends group finder with a linking length = 0 . 2 times the mean inter-particle separation within the galaxy, to identify substructures in the gas at temperatures below ≤ 8000K and densities n > 5cm -3 , within 0 . 1 virial radii of the center of the main galaxy, at ∼ 30 uniformly-spaced snapshots in time between z = 0 . 06 -0, and plot the time-averaged mass function of gas structures, MdN / d log 10 M for each simulation.  \nhere is also predicted in turbulent fragmentation models (references above) and similar to observations.  \nRemarkably, the mass function and linewidth-size relation appear independent of resolution down to clouds with just a few gas particles - we simply sample further and further down the mass function as we increase the mass resolution. Since most of the mass is in Toomre-mass structures, the total mass in clouds is also identical to within ∼ 30% in all the runs plotted (likewise for most of the turbulent power/kinetic energy). The only case which may be biased is the lowest-resolution example ( mi 1000 = 470), where the most massive clouds appear slightly more massive: this is because the clouds with ∼ 10 6 M glyph[circledot] (not quite the peak of the MF, but just below it) contribute a significant total mass but cannot be resolved (this is ∼ 2 particles), so that mass is 'shifted' into more massive structures numerically. But even in this marginal case, the gas-mass weighted mean M GMC is only over-estimated by a factor ∼ 2. In our best-resolved case, > 75% of the GMC population gas mass (and > 75% of the SFR of the entire galaxy) is contained in clouds with > 10 6 M glyph[circledot] , similar to what is observed in the Milky Way (Williams &McKee 1997).  \nNote that Hopkins et al. (2012c), using the same feedback physics (but in non-cosmological simulations using a different code), consider a much more detailed study of simulated GMC mass functions as well as size-mass relations, linewidth-size relations, virial parameters, internal column density distribution functions, lifetimes, and star formation efficiencies. Their resolution studies (and agreement with observations) are all consistent with our simple comparison here.  \nThis further demonstrates, as we argued above, that resolving the thermal Jeans mass in cold gas is not important to cloudor-larger scale dynamics (it only should matter if we are trying to resolve individual proto-stars). We will further demonstrate this below, when we consider adding artificial pressure floors to the ISM, or remove the low-temperature ( glyph[lessmuch] 10 4 K) cooling physics: because of the strong dependence of Jeans mass on pressure, these changes in turn change the thermal Jeans mass in the cold phase gas by several orders of magnitude (factors ∼ 10 3 -10 5 ), yet they have no appreciable systematic effect on any results we measure.  \nFinally, Fig. 14 compares the mass function, internal structure, and spatial distribution of subhalos within our z = 0 m12i halo, as a function of mass resolution. Down to a subhalo mass (and corresponding maximum rotation velocity) corresponding to ∼ 10 -30 DM particles, the predictions are independent of resolution, consistent with other studies. 16 Given the strong dependence of galaxy baryonic mass on halo mass, and much stricter resolution criteria for galaxy baryonic properties, the DM substructure mass functions and orbital distributions are always extremely well-resolved when we consider a galaxy resolved.", '4.2 Spatial Resolution': 'Because our simulations are Lagrangian, there is no single meaningful definition of the spatial and/or force resolution. Here we discuss different criteria, and study their importance for our conclusions.  \n16 Compare, for example, Fig. 14 here to Fig. 9 in Springel et al. (2008) (our runs with mi , 1000 = 450 , 56 , 7 correspond approximately to the mass resolution in their runs Aq-A-5, Aq-A-4, Aq-A-3). The results are similar, although the low-resolution mass functions accurately track the highresolution solution down to somewhat smaller masses here (likely owing to a combination of somewhat different force softening, and a more restrictive but accurate timestep and force-tree node opening criterion).', '4.2.1 Hydrodynamic Resolution': "The meaning of hydrodynamic resolution in our mesh-free Godunov methods is discussed extensively with dozens of numerical examples in the methods papers (Hopkins 2015, 2016, 2017; Hopkins & Raives 2016). We briefly review it here. For hydrodynamics the resolution is fully adaptive, set by the inter-particle spacing hi - it can, in principle, become arbitrarily small. In our Lagrangian method, the mean inter-particle spacing is directly related to the gas density n H, via  \nh gas i = 16pc m 1 / 3 i , 1000 ( n H 10cm -3 ) -1 / 3 (8)  \nIn Hopkins (2015); Hopkins & Raives (2016) we show that the 'effective' resolution for sound waves is identical to secondorder grid methods such as ATHENA, meaning one element ('particle') is equivalent to one cell, and the appropriate 'resolution' should be taken to be the inter-particle spacing, as opposed to the extended 'kernel search radius' (maximum distance to any interacting neighbor cell), which is more akin to the gradient stencil in grid-based codes (and depends on the detailed kernel shape). This is different from SPH methods, where properties are (by definition) 'smoothed' over a kernel. In fact, because of the higherorder reconstruction, sound waves can be accurately re-constructed to scales about ∼ 1 / 2 the inter-particle separation (Hopkins 2015 Fig. 2, Hopkins & Raives 2016 Fig. 1). Contact discontinuities can be captured and conserved across ∼ 2 -3 elements (two inter-element spacings), much sharper than non-Lagrangian codes if the discontinuity is moving (see Hopkins 2015 Figs. 16-23). Strong shocks are resolved over a similar ∼ 3 elements, close to ATHENA (Hopkins 2015 Figs. 3, 10-15, 29-30). Intrinsically multidimensional problems, such as conserving vorticity of a rotating, pressure-equilibrium disk for a few orbits are more demanding; this requires a few hundred total elements (about ∼ 10 elements 'per radian' around the structure; Hopkins 2015, Fig. 4-5). However that is easily satisfied for any disk with a resolved vertical scale-height H glyph[greaterorsimilar] 1 -2 hi . In multi-dimensional super-sonic turbulence, numerical dissipation and noise truncates the inertial range at wavelengths approximately ∼ 5 times the inter-particle spacing (comparable to AREPO and ATHENA; Hopkins 2015 Fig. 26-28).", '4.2.2 Force Softening: Definitions & Optimal Choices': "The meaning of 'gravitational spatial resolution' is somewhat ambiguous. Commonly, however, this is used to refer to the force softening glyph[epsilon1] . As discussed in § 2.2, we solve gravitational forces for the same gas-mass distribution as hydrodynamic forces: this means setting the gravitational force softening for gas adaptively,  \nglyph[epsilon1] gas i ≡ h gas i = 16pc m 1 / 3 i , 1000 ( n H 10cm -3 ) -1 / 3 (9)  \nFor clarity, we use the same definition of force softening and spatial resolution, corresponding to the inter-particle separation hi =∆ x (defined such that ρ i = mi / h 3 i ). For our default kernel, this means that the commonly quoted 'Plummer equivalent' softening is glyph[epsilon1] plummer ≈ ( 2 / 3 ) glyph[epsilon1] =( 2 / 3 ) hi . For the reasons in § 2.2 and § 4.2.2 below, we adopt constant force-softening parameters for DM ( glyph[epsilon1] DM) and stars ( glyph[epsilon1] ∗ ), with their values set as described below (for specific values, see Table 1). We define these for clarity the same way as for gas, so that they are softened according to the same kernel shape (in other words, a gas, star or DM particle with the same mass and glyph[epsilon1] will exert identical gravitational forces at all radii).  \nWe will explore variations below, but first summarize our best estimate of 'optimal' parameters here:  \n1  \nFigure 15. Spatial resolution study, as Fig. 8. We take the m10q simulation with fixed mass resolution mi , 1000 = 2, and re-run with different gravitational force softening. Left (and middle-left): Varying force softening for collisionless particles (DM glyph[epsilon1] DM and stars glyph[epsilon1] ∗ ). We multiply the default values from Fig. 8 at this mass resolution (constant glyph[epsilon1] DM = 40pc , glyph[epsilon1] ∗ = 10pc) by constant values from ∼ 0 . 5 -7 . 5; we also compare softening fixed in co-moving units (default values × 5 / ( 1 + z ) ); and we consider fully adaptive glyph[epsilon1] DM and glyph[epsilon1] ∗ . These changes have no systematic effect (small differences in mass are consistent with stochastic variations, given the bursty SFH). We have verified the same in more limited surveys of glyph[epsilon1] DM and glyph[epsilon1] ∗ in simulations m10y , m11q , m11v , m12i . Right (and middle-right): Varying force softening for gas. We compare our default, fully-adaptive softening ( glyph[epsilon1] gas = hi , the inter-particle spacing), to simulations with fixed gas softening set to 1 -200pc. Our default (adaptive) run and those with glyph[epsilon1] gas = 1 , 4pc enforce a minimum gas density for star formation n crit = 1000cm -3 (as labeled). With 1pc softening, this produces identical results to our adaptive-softening run; but with 4pc softening, the maximum gravitationally-resolved density n max =( mi /glyph[epsilon1] min gas ) / mp ≈ 1000cm -3 ( mi , 1000 / 2 )( glyph[epsilon1] min gas / 4pc ) -3 barely reaches n crit, so our self-gravity criterion cannot properly resolve if gas is selfgravitating at n > n crit and SF is artificially suppressed (slightly lower M ∗ and [ Z ∗ / H ]). Taking glyph[epsilon1] gas = 20 pc ( n max ≈ 10cm -3 ) with n crit = 100 shows this problem more severely. With fixed (non-adaptive) softening, it is necessary to choose n crit < n max ; lowering n crit = 100 for glyph[epsilon1] gas = 4pc shows excellent agreement with our default (adaptive) run. With glyph[epsilon1] gas = 200pc, n max = 0 . 01cm -3 is much less than n crit and the mean galaxy gas density - this means we cannot resolve Toomre-scale structures (large GMCs), so star formation is completely suppressed ( M ∗ = 0, exactly, here), and a super-dense (but non-fragmenting) gas disk forms.  \n<!-- image -->  \n· Gas: For gas, the optimal force softening is un-ambiguous: it should be set adaptively to match the hydrodynamic resolution. This is optimal for several reasons. (1) It is the physically correct set of equations for a collisional fluid. The hydrodynamic solver describes a mass distribution (not just at particle locations, but everywhere in the domain, according to the reconstruction of the method). We should therefore solve the Poisson equation for the same mass distribution - anything else is fundamentally ill-posed and can produce unphysical outcomes (Bate & Burkert 1997). 17 (2) Many authors have shown that this provides an optimal nu-  \nrical softening, in the sense that it automatically provides the most accurate solution (converges most rapidly) while minimizing N -body integration errors (Merritt 1996; Bate & Burkert 1997; Romeo 1998; Athanassoula et al. 2000; Dehnen 2001; Rodionov &Sotnikova 2005; Price & Monaghan 2007; Barnes 2012; Hubber et al. 2013). (3) It guarantees that numerical hard gas-gas scattering can never dominate over physical self-gravitating motions. (4) It removes ambiguity about the meaning of force/spatial/hydrodynamic resolution (it makes the mass resolution the un-ambiguous resolution scale). (5) In Lagrangian codes, spatially, self-gravitating objects will be able to collapse (correctly) to arbitrarily high densities, where some star formation criterion should identify them. This makes it essentially impossible to artificially suppress star formation by 'tuning' some star formation threshold. (6) It automatically behaves correctly in all density regimes, and naturally removes any  \nFigure 16. Spatial resolution study in re-starts of a MW-mass galaxy ( m12i ) at late times, as Fig. 10, keeping fixed mass resolution mi , 1000 = 56 but varying the spatial resolution, as in Fig. 15. All runs unless otherwise labeled take n crit = 1000cm -3 . Our default run uses adaptive softening for gas, fixed softening for DM and stars ( glyph[epsilon1] DM , glyph[epsilon1] ∗ ); the run labeled 10 × ( glyph[epsilon1] DM , glyph[epsilon1] ∗ ) increases these by 10 × , but keeps the adaptive gas softening - this has no effect on our prediction. Runs labeled glyph[epsilon1] gas fix the gas softening; for glyph[epsilon1] gas = 1 . 4pc and 14pc ( n max = 10 6 , 10 3 cm -3 , respectively), this has little or no effect. For glyph[epsilon1] gas = 140pc ( n max ≈ 1cm -3 ), we see no star formation unless we lower n crit and allow the hydrodynamic resolution hi to decrease to values much smaller than glyph[epsilon1] gas ( < 0 . 1 glyph[epsilon1] gas, here) - then shocks driven by SNe produce gas at high densities (i.e. hydrodynamic, not gravitational, effects allow the gas to reach high densities), but the most massive GMCs are still only marginally resolved. For glyph[epsilon1] gas ∼ kpc, the galaxy never forms locally self-gravitating structures, even lowering n crit = 1 and allowing hi = 0 . 01 glyph[epsilon1] gas.  \n<!-- image -->  \nambiguity about co-moving or physical softenings in cosmological integrations. 18  \nIf constant gravitational softenings for gas must be used, they must be chosen sufficiently small to resolve the vertical scaleheights of the cold gas disk, and Toomre lengths of the most massive Toomre-mass objects discussed in our mass-resolution criteria ( glyph[epsilon1] gas glyph[lessmuch] 100pc).  \n· Dark Matter: For dark matter, we will show that the force softening makes little or no difference to any conclusions, so long as it is not extremely small ( ∼ 100 × smaller than our default, which would trigger hard-scattering effects), or extremely large (which would over-soften the central DM cusp). One option is to use fully adaptive softening. However, (1) this requires stronger timestep constraints (§ 4.3.3), which considerably increase computational expense for no apparent improvement in accuracy; and (2) for a collisionless fluid, the 'correct' adaptive scaling is physically ambiguous. It is not simply the case that the spatial domain of the DM represented by particle ' a ' must shrink, if the DM particles b around a move inwards (i.e. if ∇· v a DM < 0), because those particles can move through particle a without compressing it (formally, there is no unique relation between the real-space DM N -body par-  \n18 Adaptive gas softening is the standard in most grid-based codes (Kravtsov et al. 1997; Teyssier 2002; Bryan et al. 2014) and moving-mesh codes (Springel 2010) and in particle-based codes as well in e.g. the fields of star and planet formation (see e.g. Bate & Burkert 1997). Particle-based codes in galaxy formation have proved a surprising historical exception.  \nFigure 17. Effects of force softening on the DM profile at z = 0, at fixed mass resolution. Top: DM-only simulation of m12i , with poor mass resolution ( mi , 1000 = 450). Middle: DM-only m12i , with better mass resolution ( mi , 1000 = 7). Bottom: Full-physics simulation of m10q , with mass resolution mi , 1000 = 2. Vertical dashed (dotted) lines show 'convergence radii' r 0 . 06 ( r 0 . 6), as Fig. 11; thin green line is the best-fit NFW profile in the highest-resolution simulation we have run. DM softenings glyph[epsilon1] DM ∼ 20 -1000pc ( ∼ 0 . 05 -5 times the mean inter-particle spacing inside ∼ 1kpc) have almost no effect at > r 0 . 6 and only small ( < 0 . 2dex) effects within r 0 . 06 (suppressing the densities when glyph[epsilon1] DM is too large, glyph[greaterorsimilar] 0 . 5 r 0 . 06). Mass resolution is clearly much more important, compared to force softening (compare top vs. middle ). Only the smallest glyph[epsilon1] DM = 10pc at very poor mass resolution mi , 1000 = 450 shows hard-scattering effects (the spurious cusp) - this occurs when glyph[epsilon1] DM glyph[lessorsimilar] 0 . 01 r 0 . 06.  \n<!-- image -->  \nFigure 18. Close-up of the two higher-resolution simulations from Fig. 17, showing the DM density profile in the halo center in more detail. We can now clearly see that while the profiles agree very well at large radii ( glyph[greaterorsimilar] kpc), at small (sub-kpc) radii, the runs with excessively large force softening glyph[epsilon1] DM suppress the density profile relative to the higher-resolution solutions, by a modest factor ∼ 1 . 2 ( ∼ 1 . 5) at ∼ 500pc ( ∼ 100pc). The suppression can even extend to radii of order r 0 . 6 ( dotted vertical line ). The excess suppression appears in all cases we have tested when glyph[epsilon1] DM glyph[greaterorsimilar] 0 . 5 r 0 . 06 (our 'convergence radius' r 0 . 06 is shown as the vertical dashed line) - approximately glyph[epsilon1] DM glyph[greaterorsimilar] 100pc in both cases shown here. Once glyph[epsilon1] DM is less than ∼ 0 . 25 r 0 . 06, we see no evidence for improved accuracy in the simulations (differences between these runs and adaptive-softening runs are consistent with shot noise).  \n<!-- image -->  \n8  \n6  \nticle distribution and their phase-space distribution; for discussion see Abel et al. 2012).  \nTherefore, fixed softenings appear preferable. We adopt softenings fixed in physical units at z < 10 (co-moving above this), since the halo centers do not change significantly in density over this redshift range, so this maintains an approximately fixed ratio of glyph[epsilon1] DM to the inter-particle spacing h DM i ∼ 17pc m 1 / 3 i , 1000 ( ρ DM / 10 9 M glyph[circledot] kpc -3 ) -1 / 3 .  \nFor fixed DM softenings, avoiding over-softening in the central DM profile - i.e. converging as accurately as possible to the solution of higher-resolution simulations - requires a force softening smaller than the 200-particle Power-like radius (where t relax ≈ 0 . 06 t circ ( R 200 ) ) by a factor of at least a couple, i.e. glyph[epsilon1] DM glyph[lessorsimilar] 0 . 5 r 0 . 06. This is approximately ensured if  \n3  \n3  \nFigure 19. Mass function ( top ), circular velocity distribution ( middle ), and spatial distribution ( bottom ) of satellites (DM subhalos) within the primary halo in DM-only runs of our m12i simulation (at z = 0), as Fig. 14, at fixed mass resolution ( mi , 1000 = 7) but varied force softening glyph[epsilon1] DM. Top: Number of subhalos versus bound subhalo mass. Objects are plotted down to < 30 DM particles. Middle: Number of subhalos versus maximum circular velocity. Bottom: Distribution of subhalos in radial distance from the center of the primary halo. All the distributions agree well, independent of force resolution (the small deviation in the mass function of glyph[epsilon1] DM = 60pc owes to a single subhalo which falls just outside, instead of inside, the radius cut used to identify substructures at z = 0).  \n<!-- image -->  \nglyph[epsilon1] DM < 30pc m 1 / 2 i , 1000 ( M vir / 10 12 M glyph[circledot] ) -0 . 2 (we estimate this by taking M enc ( < r 0 . 06 ) ≈ 220 mi , and assuming NFW profiles with concentration c ≈ 10 ( M vir / 10 12 M glyph[circledot] ) -0 . 15 ). The minimum softening is given by the value where we begin to see hard scattering effects cause excessive N -body relaxation in the profile and/or velocity distribution function, which is roughly glyph[epsilon1] DM > 0 . 03pc mi , 1000 ( M vir / 10 12 M glyph[circledot] ) -2 / 3 (obtained by comparing the hard-scatter ∆ v to Vc ; § 4.2.6). Note the large dynamic range in between (e.g. for our default runs, this gives 1pc glyph[lessorsimilar] glyph[epsilon1] DM glyph[lessorsimilar] 75pc). We  \nwill show we can vary glyph[epsilon1] DM by multiple orders of magnitude without changing out results. 19 To be conservative we choose softenings about ∼ 2 -3 times smaller than the upper limit above (giving glyph[epsilon1] DM ∼ h DM i in halo centers). Various tests have shown this is optimal to reduce both noise and over-softening errors (see, e.g. Merritt 1996; Athanassoula et al. 2000; Dehnen 2001; Barnes 2012; Hubber et al. 2013). 20  \n- · Stars: Stars are the most ambiguous softening case. Fortunately, like dark matter, our results appear almost completely independent of how gravity from stars is softened (§ 4.2.4). As with DM, the same physical ambiguities apply to using adaptive softenings for a collisionless fluid. Moreover, if the adaptive softening is based on the inter-star distance, a star born in a pure gas cloud would instantly 'jump' to a huge force softening; if based on inter-gas distance, a pure-stellar bulge would be over-softened out to the gaseous halo. In runs with adaptive DM+stellar softening, we therefore set glyph[epsilon1] ∗ to scale with the inter-baryon distance (gas+stars), which at least has the advantage that it handles both extremes above correctly. Even in this case, however, a 'real' (selfgravitating, resolved) star cluster can artificially expand once gas is blown away, because the neighbor search expands, even though self-gravity should fix the size (hence softenings). With constant softenings, the ambiguity is that stars form with very different densities. We have also tested glyph[epsilon1] ∗ fixed in time but variable particle-toparticle, to equal the softening of the gas from which the particle formed; while this does not change our results, it leads to wildly variable glyph[epsilon1] ∗ for stars in the same location at late times (which formed from gas at different densities at different times), which is numerically problematic. A fixed glyph[epsilon1] ∗ chosen to be large ensures smoothness of the potential but this is not actually physical for stars (because star formation is clustered, the stellar potential is 'lumpy' on small scales), and we see this again causes sudden, massive expansion of the softening after a star particle forms. A reasonable compromise is to set glyph[epsilon1] ∗ similar to the gas softening at the mean density of star formation, glyph[epsilon1] ∗ ∼ 3 . 4pc m 1 / 3 i , 1000 ( 〈 n 〉 SF / 10 3 cm -3 ) -1 / 3 , and our default choice follows this criterion.", "4.2.3 Other 'Spatial Resolution' Definitions": "As noted above, there is no single 'spatial resolution' in our simulations. In Table 3, we therefore provide a number of other 'effec-  \n19 Note that it becomes impossible to satisfy both criteria for particle masses mi corresponding to glyph[lessorsimilar] 40 DM particles in the halo; of course this is because the halos then cannot be internally spatially resolved.  \n20 Another common criteria for DM softening is that the maximal 2-body acceleration ( ai ∼ Gmi /glyph[epsilon1] 2 DM ) not exceed the bulk acceleration at some radius ( a bulk ∼ GM enc ( < r ) / r 2 ). If we assume an NFW profile, and require ai < 〈 a bulk 〉 where 〈 a bulk 〉 is the massaveraged mean acceleration of DM within the halo, we obtain glyph[epsilon1] DM glyph[greaterorsimilar] 6pc ( c / 10 ) -0 . 45 m 1 / 2 i , 1000 ( M vir / 10 12 M glyph[circledot] ) -1 / 2 (we approximate the exact scaling with the halo concentration c by a power-law, good to ∼ 10% over the entire range of interest). Our fixedglyph[epsilon1] DM choices agree reasonably well with this scaling. However, we stress that we see no measurable errors or deviations in the mass profile, rate-of-growth of structure, or velocity distribution function, even at/outside R vir, using order-of-magnitude smaller glyph[epsilon1] DM. This is because this acceleration criterion is not meaningful when glyph[epsilon1] DM is below the inter-particle separation, because only very rarely will particles approach within separations glyph[lessmuch] glyph[epsilon1] DM, and when they do, the net velocity imparted by the encounter will be given by the scaling in § 4.2.6: (∆ v /σ v , DM ) ∼ 10 -4 mi , 1000 ( 20pc /glyph[epsilon1] DM )( σ v , DM / 100kms -1 ) -2 glyph[lessmuch] 1. Equivalently, the timescale required for N -body heating at ∼ R vir to perturb a particle orbit by ∼ 10% is much greater than t Hubble for all glyph[epsilon1] DM glyph[greaterorsimilar] 1pc.  \ntive spatial resolution' values for one of our resolution studies (our m12i study) from Fig. 8, at each resolution level we have considered. Specifically we quote the 'spatial resolution' (inter-particle spacing hi ) at:  \n- (i) MaximumDensity/Minimum Softening: This is simply the minimum glyph[epsilon1] gas i ≡ h gas i reached in the simulation, at any time. This (by definition) represents the extreme in the simulation, not 'typical' values - in our highest-resolution simulations, this corresponds to densities n max ∼ 5 × 10 6 cm -3 , which are reached by the simulation only briefly in the galactic nucleus in an intense, high-redshift starburst.\n- (ii) Star-Forming Densities: We output from our simulations the gas density at which every individual star particle forms, which maps one-to-one to the equivalent hi . We show the mean density, 〈 n SF 〉 , weighted by total star formation, integrated over the entire galaxy history to z = 0, and corresponding h 〈 SF 〉 i . We emphasize that because star formation is only allowed in self-gravitating gas which is also self-shielding and molecular, the mean 〈 n SF 〉 is always significantly larger than the minimum density at which star formation is allowed ( n SF , min), even when the latter is set to n SF , min = 1000cm -3 .\n- (iii) Thermal Jeans Mass (Warm/Hot Gas): As discussed in § 4.1.3, the smallest resolveable Jeans mass simply tells us where the 'fragmentation cascade' in self-gravitating objects will be truncated (and our sub-grid star formation model will take over). We can calculate this minimum scale, but the correct 'effective Jeans mass' (the actual characteristic mass of structures) depends on turbulence, not on the sound speed, in a super-sonically turbulent medium (see e.g. Padoan & Nordlund 2002; Hennebelle & Chabrier 2008; Hopkins 2012b,a, 2013d,c,b; Hopkins & Christiansen 2013; Guszejnov & Hopkins 2015, 2016). To get some insight, however, for the warm ( glyph[greaterorsimilar] 10 4 K) gas, the turbulence is expected to be trans-sonic, so we can reasonably simplify by considering the thermal Jeans properties. We calculate the smallest possible Jeans mass resolvable with at least ∼ 10 elements at these temperatures, and the corresponding resolved Jeans length (with ∼ 3 hi ) and gas density. 21 Any lower density or larger Jeans lengths have resolved thermal Jeans fragmentation - important for our purposes, these values need to be sufficient to capture the largest GMCs (sizes ∼ 100 -200pc and densities ∼ 10cm -3 ), which form out of the warm gas. Note that if the gas is hotter, it becomes easier to resolve Jeans-scale structures.\n- (iv) Turbulent Jeans Mass (Cold Gas): Once the gas cools below temperatures ∼ 10 4 K (in e.g. GMCs or molecular disks), the turbulence is highly super-sonic - by definition, thermal pressure does not control the dynamics, and the appropriate Jeans scale for fragmentation is the turbulent Jeans scale. This replaces λ J = cs / √ G ρ by λ turb ≈〈 v 2 turb ( λ turb ) 〉 1 / 2 / √ G ρ , where v turb ( λ turb ) is the rms turbulent velocity measured within the same region. A  \n21 Define the Jeans length λ J and mass mJ = ( 4 π/ 3 ) ρ ( λ J / 2 ) 3 . Following § 4.1.1, a 'resolved' structure has some number N ∼ 10 N 10 elements, so mass ∼ Nmi . By definition of the inter-particle spacing, we also have ρ = mi / h 3 i . This gives λ J ≈ 3 hi N 1 / 3 10 , always - in other words, our adaptive softening ensures the Jeans length λ J is always spatiallyresolved provided the Jeans mass mJ is mass-resolved. Now use this and the definition of Jeans length, λ J = cs / √ G ρ , to solve for the Jeans radius λ J / 2 ≈ 9 GmiN 10 / c 2 s ; for T = 10 4 K (and assuming fully-ionized gas), this becomes ≈ 0 . 3pc mi , 1000 N 10. Taking a minimum resolved size N 10 = 1, this gives the minimum resolved λ J ; the corresponding hi gives, in turn, the corresponding maximum density n J max .  \nrigorous definition and derivation of the corresponding turbulent Jeans length, mass, and fragmentation cascade is given in Hopkins (2013b). If we assume a linewidth-size relation seen in the ISM and expected for super-sonic turbulence ( v turb ∝ λ 1 / 2 ), and that the 'parent' clouds have virial parameter of unity and (according to Larson's Laws) a universal surface density Σ cloud = Σ 300 300 M glyph[circledot] pc -2 , and take a minimum resolved-object mass ∼ Nmi , then following Hopkins (2013b) gives a minimum-resolved length λ turb ≈ 4 . 0pc ( mi , 1000 N 10 / Σ 300 ) 1 / 2 and maximum-resolved density n turb max ≈ 1 . 2 × 10 4 Σ 300 ( mi , 1000 N 10 / Σ 300 ) -1 / 2 .  \n(v) Dark Matter Power-Type Radius: None of these criteria apply to dark matter. We therefore quote the DM inter-particle separation: both its minimum value h DM , min i , as well as rms value within twice the effective radius of the baryonic galaxy at z = 0, h DM , core i . However for the DM profiles, § 4.1.4 above shows that the internal structure of a collisionless object evolved for a Hubble time is well-converged inside a radius enclosing ∼ 200 DM particles. We therefore quote this as well, for the DM.  \nQualitatively, the results from our other MW-mass series ( m12f ) are nearly identical to those in Table 3. Our dwarfs ( m10v and m10q ) have superior resolution owing to their much smaller particle masses.", '4.2.4 Effects of Collisionless Force Softening: Resolution Studies': "Figs. 15-16 consider the effects of varied force softening for both collisional (gas) and collisionless (DM, stars) particles in our cosmological simulations, for both a dwarf and MW-mass galaxy. First, consider the effects of the collisionless softening. We vary both DM and stellar softening simultaneously by multiplying both by a constant ( ∼ 0 . 5 -10) relative to their default values, setting them to be fixed in comoving (instead of physical) units at all redshifts, or setting them purely adaptively (with appropriately careful time-stepping, see § 4.3.3). 22  \nTo first approximation, we see no effect from these changes. There are some variations in the growth history in Fig. 15, but given the bursty nature of dwarf star formation histories, this appears to be primarily stochastic. Still, there is a (weak) tendency, on average, for the runs with larger softenings to produce slightly higher stellar masses (and corresponding metallicities). Although not apparent by-eye in Fig. 15, these larger-softening runs appear to have slightly less-bursty SFHs (if we quantify this by measuring the logarithmic variance in the SFR measured in 10Myr bins relative to a rolling Gyr-average value) - and in previous studies, we have shown that more-bursty SF produces more efficient outflows (because the same feedback is more concentrated; Muratov et al. 2015). We have also verified this by re-running our m11q and m11v simulations (the runs where 'burstiness' has the most dramatic effect on the DM halo structure) with 5 × larger softening ( glyph[epsilon1] DM , glyph[epsilon1] ∗ ) = ( 200 , 20 ) pc; the 'burstiness' is still obviously present but is quantitatively suppressed, and the stellar masses of both increase by a factor ∼ 1 . 4. It makes sense that such large softenings  \n22 For adaptive DM softening, we determine glyph[epsilon1] DM using the same methods as for gas, but only counting other DM as 'neighbors.' In other words, we set glyph[epsilon1] DM i = h DM i , where h DM i is determined from the dark matter particle neighbor distribution in the identical manner to h gas i , ensuring h DM i is the mean (kernel-averaged) inter-particle spacing of dark matter particles within the kernel. This gives glyph[epsilon1] DM i ≡ h DM i = 17pc m 1 / 3 i , 1000 ( ρ DM / 10 9 M glyph[circledot] kpc -3 ) -1 / 3 . For adaptive softening for stars, we follow the same exercise but include all baryonic (gas+star) particles as 'neighbors.'  \nsuppress burstiness to some extent: they smear out the DM mass profile, meaning the potential is weaker (so SF needs to be less concentrated/dramatic before it can drive outflows) and smear out any star clusters or small merging galaxy stellar components that would otherwise provide more concentrated feedback. But these are unphysical effects; furthermore, the 'level of burstiness' in the adaptive softening runs agrees well with our default simulations together this gives us confidence in our default choices. But in any case, the magnitude of this effect is much smaller than systematic uncertainties in the stellar mass predictions. Not surprisingly, for the MW-mass system, which has a more smooth SFH, the effect is minimal ( < 10%).  \nFig. 17 confirms the result of previous studies (e.g. Power et al. 2003) that at fixed mass resolution, changing the DM force softening (either in physical or comoving units or using fully-adaptive softenings) has very little effect. We see that outside the radius containing N ∼ 200 particles, these choices (with fixed glyph[epsilon1] DM varied from ∼ 10 -1000 pc) have almost no effect on the mass profile. Of course, if glyph[epsilon1] DM is too large, it will eventually suppress any smaller-scale structure (e.g. artificially flattening the DM cusp) this occurs when glyph[epsilon1] DM is larger than the convergence radius r 0 . 06. Power et al. (2003) show that convergence to the correct solution in the central structure of N -body dark matter halos also requires glyph[epsilon1] DM < c -2 [ln ( 1 + c ) -c / ( 1 + c ) ] R 200 glyph[lessorsimilar] 0 . 01 R 200; this is easily satisfied by any of our simulations that satisfy glyph[epsilon1] DM glyph[lessorsimilar] r 0 . 06. Note that most previous studies have considered only pure-DM simulations; we show here the same conclusions apply even in our 'full physics' runs.  \nFig. 17 does raise one important caveat: if glyph[epsilon1] DM glyph[lessmuch] hi (the inter-particle spacing), and the 'hard scattering' velocity deflection between two particle is comparable to their velocity dispersion, then runaway N -body effects can produce a gravito-thermal catastrophe over a Hubble time. We see this in our test with m DM = 2 . 8 × 10 6 M glyph[circledot] and glyph[epsilon1] DM = 10pc; much smaller force softening for dark matter than used in our production FIRE-2 simulations at comparable mass resolution. Considering the hard-scattering velocity from our N -body heating rate calculation below, we estimate that avoiding this requires glyph[epsilon1] DM glyph[greaterorsimilar] 0 . 02pc mi , 1000, easily satisfied in all our production simulations.  \nFig. 18 shows the central portion of the profiles from Fig. 17 in closer detail. Here we can see that there is some systematic suppression of the predicted mass profile, relative to the high-resolution converged solution, in runs with large softenings glyph[epsilon1] DM glyph[greaterorsimilar] r 0 . 06. Although the effect is small, this is expected: if the profiles are converged down to ∼ r 0 . 06 with appropriate softening, then oversoftening by setting glyph[epsilon1] DM glyph[greaterorsimilar] r 0 . 06 will necessarily smear the profile out at these radii.  \nFig. 19 shows the subhalo mass function, V max function, and spatial distribution (as Fig. 14), at fixed mass resolution but varying again the DM force softening. These are almost entirely insensitive to the DM force softening for reasonable choices. We have also examined halo formation times and internal kinematics, and find the same result.", '4.2.5 Effects of Gas Force Softening: Resolution Studies': "As noted previously, Hopkins (2015) present a large number of tests demonstrating the accuracy and near-ideal convergence rate of the implementation of adaptive gas softening in GIZMO on test problems with known solutions, including self-gravitating polytropic collapse (Evrard 1988), cosmological collapse of Zeldovich pancakes (with baryons and with/without DM), and steady-state orbit integration of stable (Toomre Q > 1) Keplerian disks, as well as  \ngood agreement with other state-of-the-art codes such as AREPO on popular code-comparison tests such as the adiabatic 'Santa Barbara Cluster' (Frenk et al. 1999). Of course, while necessary, idealized tests do not ensure ideal results in complicated multi-physics simulations like those here, so we explore changes to the gas softening here.  \nIn Figs. 15-16, we also considered the effects of changing the gas force softening, replacing our default self-consistent adaptive softening ( glyph[epsilon1] gas = hi ) with a fixed, constant physical softening ( glyph[epsilon1] gas). In both our dwarf and MW-mass simulations, we see, reassuringly, that for sufficiently small fixed glyph[epsilon1] gas, the differences are essentially negligible (entirely consistent with stochastic fluctuations).  \nHowever, recall that our star formation model is based on identifying self-gravitating gas above some density threshold. At fixed mass resolution, a fixed glyph[epsilon1] gas sets a minimum inter-particle separation, hence maximum gas density, at which self-gravity will be correctly calculated: this density is n max = m -1 p mi ( glyph[epsilon1] min gas ) -3 ≈ 1000cm -3 mi , 1000 ( glyph[epsilon1] min gas / 3 . 5pc ) -3 . So not surprisingly, when we make glyph[epsilon1] gas large enough that n max glyph[lessorsimilar] n crit (the minimum density for star formation), we can artificially suppress the SFR. Obviously, one should not therefore, with fixed glyph[epsilon1] gas, choose n crit glyph[greaterorsimilar] n max . If we lower n crit sufficiently so that n crit glyph[lessmuch] n max , then we recover nearlyidentical behavior to our default simulations with adaptive glyph[epsilon1] gas and larger n crit.  \nThis is important: some simulations allow n crit glyph[lessmuch] n max (e.g. Guedes et al. 2011; Shen et al. 2014 where glyph[epsilon1] gas ∼ 100 -200pc, n max ∼ 0 . 1cm -3 , and n crit ∼ 5 -100cm -3 ) - external forces (e.g. shocks) can still produce n > n crit (if hi is allowed to be smaller than glyph[epsilon1] gas), but this means star formation has nothing to do with self-gravity locally. We therefore consider two similar experiments ( m10q with glyph[epsilon1] gas = 20pc, n max = 10cm -3 , and m12i with glyph[epsilon1] gas = 140pc, n max = 1cm -3 ; both with n crit = 100cm -3 ). Surprisingly, these still behave reasonably. Both galaxies have large turbulent motions an external perturbations (e.g. mergers) that produce a broad distribution of density fluctuations; 23 the overdensities cannot 'detach' from the turbulent flow and collapse to still higher densities (like real GMCs) if glyph[epsilon1] gas is too large (which means the properties of the cold, dense gas will be incorrect), but they can still reach n crit, at which point the galaxy-averaged SFR is self-regulated by a balance of inflow and feedback. In the m10q ( glyph[epsilon1] gas = 20pc) run, the too-large glyph[epsilon1] gas still suppresses the SFR and stellar mass by a factor ∼ 2 . 5, and makes the SFR artificially bursty (more like a much lower-resolution simulation), while the effects are less dramatic in our m12i experiment. This is expected as dwarfs have much lower turbulent Mach numbers, so the mechanism above requires some large-scale perturbation (it appears that gas accumulates until the disk goes globally gravitationally unstable, causing excessively large bursts in the galactic center instead of local fragmentation into clouds). Of course, if we make glyph[epsilon1] gas still larger, we eventually suppress any SF or ISM substructure/GMCs, giving an unphysical warm gas-pressure supported, non-self-gravitating disk.", '4.2.6 N-Body Heating Rates': "One natural concern is that N -body heating of gas/star particles by the dark matter (or stars, or other gas particles) could thicken the stellar disk, or inject spurious thermal energy/turbulence into  \nthe gas. Consider the worst-case hard-scattering scenario: two particles approach one another in a vacuum (no other particles nearby) and scatter near the particle centers, with encounter velocity v encounter. Since the duration of the encounter is short, the impulsive velocity change is ∆ v = | ∆ v | ∼ Gmi / ( hi v encounter ) (a proper integration over the exact kernel shape used in the simulations gives a maximal deflection ≈ 1 . 3 Gmi / ( hi v encounter), which occurs for encounters with impact parameter ≈ hi ). For gas, this is ∆ v ∼ 0 . 01kms -1 m 2 / 3 i , 1000 ( n / cm -3 ) 1 / 3 ( v encounter / 10kms -1 ) -1 . For dark matter, the particle masses are larger, but so are the softenings and typical encounter velocities since the dark matter has approximately isotropic dispersion at the virial velocity, ∼ 100kms -1 . So ∆ v ∼ 1 . 3 Gm DM / ( h DM v encounter ) ∼ 0 . 017kms -1 mi , 1000 ( h DM / 20pc ) -1 ( v encounter / 100kms -1 ) -1 , only slightly larger even when we adopt the smallest DM softening seen in our simulations. 24  \nWe can translate this to an 'N-body heating rate,' assuming each encounter adds to the velocity incoherently, with particle encounter rate ∼ v encounter / hi . Accounting for the number of baryons per gas particle ( ∼ mi / mp ), this gives a heating rate (rate-of-change of kinetic energy) of:  \nQ gas -gas heat 3 -1 ∼ 2 × 10 -32 mi , 1000 ( 10kms -1 )  \nQ gas -DM heat ergcm 3 s -1 ∼ 6 × 10 -32 mi , 1000 ( 100kms -1 v encounter )( ρ DM ρ gas )  \nergcm s v encounter (10) (11)  \nThese should be compared to the physical heating/cooling rates of the gas, Λ ∼ 10 -23 -10 -22 ergcm 3 s -1 for typical ISM conditions; they are far smaller than almost any other source of error in the baryonic physics, even for the lowest-resolution simulations in our tests. 25  \n24 In fact, the ∆ v in the text above is actually a significant over-estimate when we use adaptive gravitational softenings, as particle-particle encounters are never 'in a vacuum.' While a particle a traverses a domain h within the kernel of particle b , it of course feels some acceleration. But provided there is actually mass being represented by b , this is completely physical. The numerical noise/error depends not on the absolute magnitude of the acceleration from b onto a , but on the deviation of the potential from b , owing to finite sampling, from the smooth potential that would be represented if we had infinite resolution. This is straightforward to estimate. Following Dehnen & Aly (2012), the acceleration from each particle is constructed assuming its mass distribution follows the kernel function; if we assume the correct background is a uniform density field, then discretize this into particles, it is straightforward to compute the fractional deviation from the correct (infinite-resolution) solution for different particle configurations within the kernel. For our standard cubic-spline kernel, this is essentially what is shown in Fig. 3 of Dehnen & Aly (2012); the typical deviation considering various different particle configurations is about | ∆ ρ | / | ρ | ∼ 0 . 005. Thus the actual hard-scattering amplitudes are suppressed by a factor of ∼ 100 from their already-small values.  \n25 We have directly verified this in a series of numerical tests: setting up a periodic box with dark matter particles (matched to the particle mass, velocity dispersion, and space density of the cosmological simulations at the radii within the halo where the N -body heating in Eq. 11 is maximized), and an equilibrium, isothermal, non-radiative gas disk. In tests where the gas feels only DM, we measure Q gas -DM heat a factor of a few smaller than Eq. 11; in tests where the gas feels self-gravity and is 'stirred' (so some relative gas-gas motion exists), the heating contribution from Q gas -gas heat is unmeasurably small.  \nρ  \n∆  \nFigure 20. Dangers of using adaptive force softening for collisionless (nongas) particles, without careful timestepping. We show the z = 0 mass profile in DM-only simulations (here our m10v run, but we repeated a subset of these tests in m09 , m10q , m11q , m11v , and m12i ), at a series of resolution levels (labeled). For our default (fixedglyph[epsilon1] DM) runs, we show the conservative ∼ 2000-particle Power et al. (2003) radius r 0 . 6, but plot the full profile down to our 'convergence radius' r 0 . 06 - the profiles in all previous tests agree well outside this radius, independent of glyph[epsilon1] DM, and we see the same in the fixedglyph[epsilon1] DM runs plotted here. We compare our adaptivesoftening implementation with three different timestep choices. If we use only the same timestep criteria employed for fixed-softening runs ('fixedglyph[epsilon1] ∆ t limiter only'), the simulations clearly exhibit a spurious, non-converged 'cusp' at the center indicative of integration errors that arise because particles 'move through' one another in a timestep, preventing the adaptive terms from properly adapting (this cannot occur with gas, owing to the Courant condition). Lowering the timestep by a uniform factor of 30 resolves the problem, and produces excellent agreement with the fixedglyph[epsilon1] runs, but at great computational expense. Our default implementation uses the 'new ∆ t limiter' described in § 4.3.3 (Eq. 13) to control the timesteps with adaptive softening for collisionless particles: this also produces excellent agreement with fixedglyph[epsilon1] runs.  \n<!-- image -->", '4.3 Time Resolution': 'The time resolution in the simulations is set by the timestep ∆ t , which is always set according to the minimum of the various criteria in § 2.7. Together, these criteria ensure that gravitational and fluid dynamics, as well as stellar evolution, are always explicitly time-resolved. Usually, the acceleration-based criterion ∆ t < 0 . 2 ( hi / | a i | ) 1 / 2 from Power et al. (2003) is the most demanding, for both gas and N -body particles (because the dense gas is highly super-sonic). Occasionally, however, in gas shock-heated by many SNe, the Courant criterion is most important.  \nIn either case, the minimum timesteps reached in our simulations are ≈ 100yr (see Table 3). This is reached regularly by some elements from redshifts z ∼ 0 -2, although always by a very small fraction of the total population at any given instant.', '4.3.1 Standard Criteria & Results of Variations': 'Wehave tested variations in our standard time-stepping criterion, in a limited sub-set of simulations with both dark matter and baryons: the DM-only runs are run to z = 0, but baryonic runs are only integrated to z = 4, because inappropriate timesteps usually lead to catastrophic numerical instability which is evident quickly. Our conclusions are identical to canonical studies in the field. Like  \nPower et al. (2003), we find that a coefficient smaller than α = 0 . 2 in the timestep limiter ∆ t < α ( hi / | a i | ) 1 / 2 produces no appreciable gain in accuracy, but much larger values ( α glyph[greaterorsimilar] 0 . 5) can seriously degrade mass profiles and orbit integration (e.g. angular momentum conservation); our gravity solvers are very similar, following Springel (2005), so this should not be surprising, especially for DM-only tests. Even though this acceleration-based criterion usually dominates for gas (recall a i includes all accelerations, i.e. hydrodynamic and gravitational), a Courant-type condition is still necessary, with coefficients α > 0 . 8 in ∆ t <α hi / v max sig , i giving rise to numerical instability, while values α < 0 . 4 (our preferred value) do not produce any obvious improvement. This is consistent with standard hydrodynamic tests and analytic numerical stability analysis (see Hopkins 2015 for discussion).  \n4.3.2 Stellar-Evolution Timestep Limits  \nRecall, we do not allow the timestep for star particles to exceed  \n∆ t ∗ < MAX ( 10 4 yr , t ∗ 300 ) (12)  \nwhere t ∗ is the age of the star, to ensure stellar evolution is timeresolved and the expected number of SNe (per star particle per timestep) is always < 1 at our production resolution. We have experimented with weakening this limiter, increasing it by a factor of ≈ 3 (which produced no measureable effect), or most radically using ∆ t ∗ < 20Myr m -1 i , 1000 - for mi , 1000 glyph[greaterorsimilar] 1, this is equivalent to allowing as many as ∼ 10SNe per timestep in an extremely young star particle. We did not see any significant difference running tests of m10q and m12i at low resolution (the difference in timestep is maximized at low resolution) to z = 2; however, we found that usually the number of SNe exploding at once was much smaller ( ∼ 2 -3), because the youngest stars are in dense regions where other timestep limits (e.g. the acceleration criterion above) impose much stricter limits than ∼ Myr. However, if we removed the limiter entirely ( ∆ t ∗ →∞ ), we saw clear (albeit rare) pathological activity: for example lone star particles formed in poorly-resolved dwarfs in the outskirts of the high-resolution region or ejected via tidal interactions might be assigned extremely long dynamical timesteps and have a huge single-timestep injection of SNe, continuous stellar mass loss, and radiation (giving for example a small number of gas elements with artificially high metallicity).', '4.3.3 Adaptive Force Softening for Collisionless Particles: The Need for Additional Timestep Criteria': "As described in § 2.2, for a subset of our tests presented in § 4.2.4, we use adaptive softening for collisionless (DM and stellar) particles. Variations of these methods have been explored in a number of studies on collisionless, self-gravitating systems (Athanassoula et al. 2000; Price & Monaghan 2007; Bagla & Khandai 2009; Iannuzzi & Dolag 2011; Barnes 2012; Iannuzzi & Athanassoula 2013). We have validated our numerical implementation with all tests presented in Price & Monaghan (2007) and Iannuzzi & Dolag (2011).  \nHowever, introducing these adaptive softenings requires stronger timestep criteria for collisionless particles, especially in dense halo centers where DM particles on highly-radial orbits may 'plunge' and interact with particles of widely-differing glyph[epsilon1] over the course of their orbits.  \nFig. 20 shows the results of using adaptive softening in DMonly simulations ( m10v ) of varying resolution, without imposing any additional timestep limiter beyond what is described above for collisionless particles with fixed softening glyph[epsilon1] . Clearly, the central sub-kpc regions of the DM profile exhibit a spurious density 'cusp' well in excess of the converged (much higher-resolution)  \nresults from constantglyph[epsilon1] DM runs. We confirm that this feature is artificial by simply re-running the simulations enforcing a factor of 30 smaller timestep; this agrees well with the fixedglyph[epsilon1] runs. We have also confirmed that both sets of runs maintain global conservation (as expected) - the error appears to be associated with local integration errors between particle neighbors when the timesteps are too large.  \nWith fixed softening, there is no Courant-like condition required for collisionless particles: since the equations of motion depend only on the collective long-range forces (the gravitational potential), if the potential is sufficiently smooth ( | a | is small), particles can safely 'move through' one another (Monte Carlo-sampling the phase-space distribution function). But this is not true with adaptive softening (again, for collisionless particles), because the softening length (hence self-gravity) of a particle depends on the local neighbor configuration, introducing local correction terms that must be integrated smoothly as particles move through one another - otherwise these terms are under-sampled and effectively scatter the orbits.  \nFor gas, this error does not occur because of the Courant condition. The solution for collisionless particles with adaptive softening is therefore straightforward: we should implement a similar timestep criterion. We require:  \n∆ t a AGS < 0 . 25MIN { 1 |〈 ˜ ∇· v 〉 a | , glyph[epsilon1] a v a sig , AGS } v a sig , AGS ≡ MAX ( b : | x | ba < Ha , H b ) {| ( v b -v a ) · ˆ x ba |} ˜ ∇· v 〉 a ≡ ∑ b : | x | ba < Ha ( v b -v a ) · ∇ W ( x ba , Ha ) Ω a ∑ b W ( x ba , Ha )  \n〈  \nHere Ha ≈ 2 glyph[epsilon1] a (with glyph[epsilon1] a = ha ) is the domain of the nearest-neighbor search around particle a , and the sums over b represent sums over all interacting neighbors of the relevant particle type such that they contribute to defining the inter-particle spacing and softening length glyph[epsilon1] of particle a . By analogy to the Courant condition, we define v a sig , AGS as the maximum approach/recession velocity of any neighbor within this interaction kernel - this requirement is simply that two particles cannot 'cross' more than ∼ 1 / 2 their relative softening lengths in a single timestep. The particle-divergence 〈 ˜ ∇· v 〉 a is defined exactly the same as the traditional SPH velocity divergence, 26 in such a manner that this requirement prevents glyph[epsilon1] a from changing by more than ∼ 10% within a single timestep. Because we use adaptive timesteps, we also enforce a 'wakeup' condition identical to that used for the hydrodynamics (see Saitoh & Makino 2009; Durier & Dalla Vecchia 2012); specifically, if an 'active' particle in a sub-step interacts with an 'inactive' particle in a timestep > 4 times larger ( v sig , AGS from its previous active step > 4 times larger), the inactive particle is stopped from taking the larger timestep and moved to the smallest active timebin in the active hierarchy. This prevents particles moving very rapidly from artificially moving 'through' a particle with a long timestep.  \n26 With adaptive softening for DM, in 3D, we define glyph[epsilon1] a = ha = ( 3 N eff / 4 π ) 1 / 3 Ha such that, in 3D, h 3 a ¯ na = 1 and 4 π/ 3 H 3 a ¯ na = N eff (with N eff = 32 for our standard cubic spline kernel), where ¯ na ≡ ∑ W ( x ba , Ha ) is a kernel-averaged particle neighbor number density. With some straightforward algebra, this gives an exact discrete equation for the Lagrangian derivative of ha , Dha / Dt = -( ha / 3 Ω a ¯ na ) ∇ ¯ na · ∂ ( x b -x a ) /∂ t = -( ha / 3 ) 〈 ˜ ∇· v 〉 a , where Ω a ≡ 1 +( ha / 3 ¯ na ) ∂ ¯ na /∂ ha (for derivations, see Price 2012; Hopkins 2013a). So if we do not want the softening ha to change by more than a factor | (∆ ha ) / ha | < α in one timestep, we require a timestep ∆ t < 3 α/ |〈 ˜ ∇· v 〉 a | .  \n(13)  \n1  \nFigure 21. Effects of the method for solving the hydrodynamic equations, in cosmological simulations of both a dwarf ( m10q ; left ) and MW-mass galaxy ( m12i ; right ), as Fig. 4. We compare our default method in GIZMO - MFM (solid lines), a higher-order accurate, mesh-free finite-volume Godunov method - to the GIZMO implementation of smoothed-particle hydrodynamics (SPH) which was used for FIRE-1 (dotted lines). GIZMO is a multi-method code so we can change the hydro solver while keeping all other physics and numerics identical here. We repeat each comparison at multiple resolution levels. For dwarfs, we see excellent agreement between MFM and SPH, at every resolution level. For MW-mass systems, on the other hand, FIRE-1 SPH predicts somewhat lower SFRs, stellar masses, and central Vc . The critical difference between dwarfs and massive galaxies is likely to be the presence of the 'hot gaseous halo' around massive systems (absent in dwarfs), which determines the cooling rate onto the galaxy and into which galactic winds propagate. Known issues in SPH can suppress fluid mixing in the halo, even with 'state of the art' SPH formulations, which in turn leads to easier escape of winds, less efficient cooling, and (in turn) lower masses.  \n<!-- image -->  \nFig. 20 demonstrates that this timestep limiter cures the errors seen before. The added timestep criterion does add a significant cost to the DM-only run, though far less costly than the uniform factor ∼ 30 smaller-timestep case, since only a small number of particles are affected at each time. In any case the agreement between adaptive softening with appropriate timesteps and fixed-softening runs is excellent down to r 0 . 06. We have also recently learned that other authors who have implemented adaptive DM softening following Price & Monaghan (2007) have reached the same conclusions and found it necessary to include similar timestep criteria to maintain numerical stability (V. Springel, private communication).  \nFigure 22. Effects of the hydrodynamic method in our cosmological m12i simulation at two resolution levels ( left and right ), as Fig. 21. Here we compare different 'flavors' of SPH (see § 5.1.1): (a) the FIRE-1 implementation from Fig. 21, which uses the pressure-energy formulation and a larger smoothing kernel, designed to reduce fluid-mixing errors; (b) the FIRE-1 SPH model with 'stronger mixing' - an explicitly increased thermal energy/entropy mixing term (larger 'artificial conductivity'); (c) D-SPH: a simpler SPH implementation which uses the 'density-energy' formulation of the equations and a smaller kernel (but still uses the higher-order artificial viscosity and conductivity switches of FIRE-1 SPH), which strongly suppresses the ability of the method to capture fluid mixing instabilities. The 'FIRE-1 SPH + Stronger Mixing' run agrees well with MFM; the 'D-SPH' run much more strongly suppresses the SFR, stellar mass, and central Vc , in a manner which appears to diverge with resolution (consistent with the fact that the SPH errors are zeroth-order). This demonstrates that the effects of SPH on fluid mixing physics dominate the differences between runs in Fig. 21.  \n<!-- image -->  \n2  \n2", '5.1 Finite-Volume Godunov Methods (MFM) vs. SPH': "As noted in § 2.1, a range of idealized test problems demonstrate that our default, finite-volume MFM method for the hydrodynamics is significantly more accurate than SPH, at fixed mass resolution. However, given that gravity and feedback overwhelmingly dominate over pure hydrodynamic forces in the simulations here, does the hydrodynamic accuracy matter? There have been many studies arguing that the treatment of feedback is more important than details of the hydrodynamic solver in this context (Scannapieco et al. 2012; Power et al. 2014; Hu et al. 2014; Few et al. 2016; Davé et al. 2016; Zhu & Li 2016; Stewart et al. 2017; Kim et al. 2016), although in some of these cases it is difficult to isolate the effects of the hydrodynamics within the context of a single feedback model.  \nFigure 23. Effects of the hydrodynamic method on the SFR at low-redshift of a re-started MW-mass ( m12i ; mi , 1000 = 56) simulation as Fig. 10. We compare MFM and various SPH 'flavors' as Fig. 21. For fixed initial conditions (early times) the SFR is feedback-regulated and identical. The two slowly diverge by ∼ 30% over ∼ Gyr as the SPH flavors with less accurate fluid-mixing treatments are able to more easily eject gas from the galaxy. However the effect is small compared to the cosmologically time-integrated effects seen in Fig. 21. This is consistent with the idea that the mixing and re-cycling of these outflows, and subsequent hot gas cooling (which occurs on timescales of order the Hubble time), not the generation of outflows or self-gravitating fragmentation and star formation in the dense galactic disk, is the dominant reason for the difference between certain SPH flavors and MFMin massive galaxies in Fig. 21.  \n<!-- image -->  \n<!-- image -->  \nFigure 24. Mock images, as in Fig. 3, showing the effects of the hydrodynamic solver on the z = 0 morphology of the simulated galaxies, for the highest-resolution ( mi , 1000 = 56) MW-mass ( m12i ) systems in Fig. 21. Dwarfs are not shown, as they have irregular morphologies independent of the hydrodynamics or resolution (as expected). The large-scale visual morphology is similar, although the MFM simulation exhibits a slightly more compact, thinner disk and more well-ordered spiral structure. The major differences in stellar mass and central circular velocity are not obvious in the stellar morphology. The same is true for all SPH flavors in Fig. 22.  \n<!-- image -->  \nFigure 25. Gas morphology around the galaxy in the m12i simulations (with mi , 1000 = 56) from Fig. 21, at z = 0 . 9 in 200 kpc boxes, with different hydrodynamic methods (as Fig. 21; labeled). Images are logarithmically-weighted surface-density projections, with red/green/magenta showing hot ( > 10 6 K), warm ionized ( ∼ 10 4 -10 5 K), and cool neutral ( < 8000K) gas. The most dramatic differences appear in the CGM; we choose a time where a merger has triggered violent outflows of cool gas to maximize these differences. In some less-accurate SPH formulations, such as 'D-SPH' here ( right ; density-SPH with smaller smoothing kernel), the outflow has broken into the well-known artificial 'SPH blobs' that result from errors treating fluid mixing interfaces. In our non-SPH, finite-volume MFM method ( left ), these are absent. The FIRE-1 implementation of SPH ('P-SPH' with larger kernel) was specifically formulated to reduce the fluid-mixing errors in SPH; it dramatically reduces (but does not completely eliminate) the 'blobs' (FIRE-1 SPH+Stronger Mixing closely resembles this, but with slightly fewer 'blobs'). Also, in MFM the hot halo gas is more compact/dense, with a sharper shock front, while in SPH it is lower-density with an extended boundary, owing to difficulties in shock-capturing and numerical dissipation/mixing (Bauer & Springel 2012). These effects in SPH make it 'easier' for cold winds to be ejected and avoid mixing in the halo, lowering the CGM gas cooling rate, explaining the suppression in SF in Figs. 21-23.  \n<!-- image -->  \nFigure 26. Comparison of hydrodynamic methods (MFM and SPH 'flavors' from Fig. 22) in the temperature-density phase diagram (at z = 0) in the MW-mass ( m12i , with mi , 1000 = 56) systems in Figs. 21-25. We compare all gas inside the virial radius ( top ) and gas within the CGM (excluding the galaxy; bottom ). Colors are a mass-weighted heat map with the density of gas mass in the space shown increasing logarithmically from black-blue-green-yellow-red-white. The 'D-SPH' (simpler, less accurate fluid-mixing) flavor produces no hot halo gas with densities n glyph[greatermuch] 10 -4 cm -3 , owing to a combination of poor shockcapturing, numerical dissipation, and artificial 'ease' with which cold 'blobs' in shredded galactic winds escape the halo (instead of shocking); there is also a substantial population of gas with temperatures T ∼ 10 -1000K and densities n ∼ 10 -2 -10 3 cm -3 in the CGM, from the same 'blobs.' The more accurate fluid mixing treatment in the FIRE-1 SPH mostly eliminates the cold CGM blobs and restores hot dense gas - but there is still a significant amount of gas at ∼ 10 4 K even at very low densities which disappears when we add 'stronger mixing,' in agreement with our MFM simulations.  \n<!-- image -->  \nFortunately, GIZMO is an inherently multi-method code, so we can compare simulations with otherwise identical physics and numerics, replacing only the hydrodynamic solver.  \nFigs. 21, 22, 23, & 24 compare GIZMO simulations using both  \nour (1) default (non-SPH) MFM method, 27 (2) various 'flavors' of SPH.  \n27 We have also considered a limited comparison of the 'meshless finite volume' (MFV) method in GIZMO. This differs from MFM only in a second-order mass-flux term, and gives nearly-identical results in test prob-", "5.1.1 'Flavors' of SPH": "There is an extensive literature of various SPH 'flavors' which attempt to reconcile the differences between SPH and finite-volume methods like our MFM here (see e.g. Wadsley et al. 2008; Price 2008; Saitoh & Makino 2013; Hopkins 2013a; Read & Hayfield 2012; Hu et al. 2014; Rosswog 2015, and references therein). In this section we will consider three such flavors:  \n- · FIRE-1 SPH: This is the SPH formulation used in FIRE1, and described in Hopkins et al. (2014). Briefly, we use the 'pressure-energy' formulation of SPH (P-SPH) from Hopkins (2013a) to eliminate the spurious 'surface tension' error present in 'density-energy' (D-SPH) formulations at contact discontinuities; a higher-order kernel (the quintic spline; with ∼ 64 effective neighbors instead of ∼ 32 for the cubic spline in MFM) to reduce zeroth-order SPH errors (see Dehnen & Aly 2012; Zhu et al. 2015); more accurate gradient estimators using the moving-leastsquares approach (as in García-Senz et al. 2012; Rosswog 2015); higher-order switches to minimize the artificial viscosity following Cullen & Dehnen (2010) with further improvements described in Hu et al. (2014) and Hopkins (2015); and added 'artificial conductivity' terms to allow entropy/thermal energy diffusion similar to artificial viscosity (Wadsley et al. 2008; Price 2008). These improvements to SPH are designed to improve behavior in fluidmixing and shock capturing, while reducing noise and artificial numerical diffusivity away from shocks (see references above).  \nReaders interested in more details of the SPH method here should consult the public version of GIZMO: this is the default SPH implementation.  \n- · FIRE-1 SPH + Stronger Mixing: Acommon feature in most of the SPH flavors above is the use of an 'artificial conductivity' term (as in FIRE-1 SPH) which allows for diffusion of entropy between particles: without this, entropy becomes 'particle-locked' leading to artificial resolution-scale discontinuities that suppress fluid mixing instabilities.  \nWith these 'artificial diffusion' terms in SPH, there is considerable freedom to adjust the form and normalization of the diffusivity. While the prescriptions are usually tuned to give some desired balance of accuracy on a mix of different idealized test problems, there is usually no single 'correct' prescription. We therefore consider the effects of adjusting this diffusivity term to give 'stronger mixing' in SPH, in order to eliminate the spurious 'blobs' discussed below. Specifically, we take the artificial conductivity in the GIZMO SPH implementation to be:  \nd E SPH a dt = ∑ b α ab ˜ v sig ab mamb ¯ ρ a +¯ ρ b ( ua -ub ) ˜ Wba (14) ˜ Wba ≡ ∂ W ( | x ba | , ha ) ∂ | x ba | + ∂ W ( | x ba | , hb ) ∂ | x ba | ˜ v sig ab ≡ cs , a + cs , b -3 ( v b -v a ) · ˆ x ba α ab ≡ Θ( ˜ v sig ab ) | Pa -Pb | Pa + Pb    α V a + α V b 8 ( 'Default SPH '' ) 1 ( 'Stronger Mixing '' ) Θ( ˜ v sig ab ) ≡ { 0 ( ˜ v sig ab ≤ 0 ) 1 ( ˜ v sig ab > 0 )  \nlems (Hopkins 2015; Hopkins & Raives 2016), so (unsurprisingly) the results are very similar to MFM. We prefer MFM as our 'default' method because it maintains element masses, minimizing particle-splitting and reducing N -body noise.  \nwhere x ba ≡ x b -x a is the particle separation, m mass, ¯ ρ the SPHestimated density, W the SPH smoothing kernel as a function of x ba and smoothing length ha , cs the sound speed, v the velocity, P the pressure, and α V a is the artificial viscosity coefficient defined in Hopkins 2015, Appendix F2, Eq. F16-F17. This α V varies between 0 . 05 -2 depending on the velocity divergence ∇· v and its time-derivative, reaching large values when particles approach increasingly rapidly ( ∇· v < 0 and d [ ∇· v ] / dt < 0) and decaying rapidly otherwise.  \nThe difference between our 'default' FIRE-1 SPH and 'Stronger Mixing' simulations, therefore, is that 'FIRE-1 SPH' (by design) only allows for mixing/diffusion of entropy and thermal energy when particles are effectively inter-penetrating or 'move through one another.' In the 'Stronger Mixing' case, we make the effective diffusivity/mixing stronger in those cases (by a factor ≈ 4), but more importantly, we allow mixing in shear flows. The latter case has been shown to improve the accuracy of SPH in modeling fluid mixing instabilities such as the Kelvin-Helmholtz instability (Price 2008) and more accurately represents the 'unresolved turbulent mixing' terms advocated in SPH by Wadsley et al. (2008); it is also numerically closer to the implementation advocated by Read &Hayfield (2012); Hobbs et al. (2013), who have argued this gives better convergence modeling cores of cooling-flow halos (Power et al. 2014). On the other hand, the 'Stronger Mixing' implementation tends to 'smear out' shock fronts and Keplerian shear flows much more significantly, further reducing the 'effective resolution' of SPH (Hu et al. 2014; Few et al. 2016) - this is why it was not our 'default' approach in FIRE-1. Note, though, that even in the 'Stronger Mixing' case, there is still a 'switch' α ab which suppresses diffusion in supersonically receding flows and sharp phase discontinuities.  \n- · Density-SPH + Smaller Kernel (D-SPH): Alternatively, we can compare an implementation of SPH where we remove some of the improvements used in our FIRE-1 SPH method. Specifically, in these simulations, we use the 'density-energy' SPH (D-SPH) formulation instead of the 'pressure-energy' SPH formulation, which introduces a non-convergent (sub-zeroth-order) error that has the functional appearance of a surface tension force at phase discontinuities (see Saitoh & Makino 2013; Hopkins 2013a). We also reduce the size of the SPH smoothing kernel, using a cubic spline with ∼ 32 effective neighbors; this increases the zeroth-order SPH errors and reduces the ability of the code to capture fluid-mixing instabilities. This is the 'T-SPH' formulation studied in Hopkins (2015). We retain all other aspects of FIRE-1 SPH (specifically the higher-order gradient, artificial viscosity, and artificial conductivity estimators). We consider this implementation because it allows us to specifically highlight the effects of SPH errors that suppress fluid mixing at contact discontinuities.", '5.2 (Weak) Effects in Dwarfs & the ISM of Massive Galaxies': "Fig. 21 shows that at all resolution levels, our dwarf galaxy simulations are very similar in SPH and MFM methods (at all redshifts). It also clearly shows that SPH and MFM agree well on the properties of the main progenitor of massive (MW-mass) galaxies while the progenitor is still a dwarf with stellar mass M ∗ glyph[lessorsimilar] 10 10 M glyph[circledot] at redshifts z glyph[greaterorsimilar] 1 . 5 -2. This holds for all properties plotted and all others we have examined: e.g. the distribution of mass in different phases, CGM gas morphology, covering fraction of absorbers, metallicity gradients, and galaxy rotation.  \nThis is not surprising: such galaxies are in the 'cold accretion' regime where there is little or no 'hot halo' of virial shock-heated  \ngas surrounding the galaxy (e.g. Kereš et al. 2005). Rather, gas falls onto galaxies on a free-fall time before being expelled by feedback. As such the properties are simply regulated by a combination of feedback and self-gravity; subtleties of fluid mixing (where SPH differs most dramatically from other hydrodynamic methods) are unimportant. This is consistent with previous studies in FIRE-1 using different 'flavors' of SPH, which had no effect on dwarf galaxy properties (Hopkins et al. 2014).  \nFigs. 21-23 show that for MW-mass galaxies, which do have hot halos, the qualitative behavior is similar in MFM and SPH, but some significant quantitative differences appear (see below). However, Fig. 23 shows that if we re-start an identical initial condition in MFMand SPH at low redshift, the SFRs are initially identical, and only slowly drift apart (here by just ∼ 30% over ∼ 1Gyr). Moreover Fig. 24 shows that the visual galaxy morphologies in SPH and MFMare very similar.  \nThese results are consistent with a series of studies of noncosmological simulations of MW-mass galaxies, comparing different SPH flavors (Hopkins et al. 2012b, 2013c; Hu et al. 2014; Kim et al. 2016) and SPH versus moving-mesh codes (Hayward et al. 2014). These studies showed that within galaxies known SPH errors have little effect on predictions, since the turbulence is primarily super-sonic, cooling is fast compared to dynamical times, and phase structure is primarily driven by gravitational collapse and feedback - all limits where SPH performs well (see Kitsionas et al. 2009; Price & Federrath 2010; Hopkins 2015).", "5.3 (Significant) Effects in the CGM of Massive Galaxies' 'Hot Halos'": "However, despite the similarity in the SFRs in re-started massive galaxies and non-cosmological simulations, Figs. 21-22 clearly show that the cosmological SFHs of MW-mass systems diverge between some SPH flavors and MFM, as they reach stellar masses M ∗ glyph[greaterorsimilar] 10 10 M glyph[circledot] , corresponding to halo masses glyph[greaterorsimilar] 10 11 . 5 M glyph[circledot] . In the D-SPH runs (which maximize the difference), the final stellar mass is suppressed relative to the MFM run by a factor of ∼ 1 . 5 (low resolution) or ∼ 3 (intermediate resolution), giving a correspondingly smaller peak circular velocity.  \nSimilarly, Davé et al. (2016) compare low-resolution cosmological large-volume simulations using a simplified, non-FIRE subgrid ISM and feedback model and MFM and 'P-SPH' (FIRE-1 SPH) methods in GIZMO. Over the mass range simulated ( M ∗ ∼ 10 9 -10 12 M glyph[circledot] ), the stellar mass functions agree very well at z > 1, but then begin to differ, with SPH tending towards smaller masses by ∼ 0 . 2dex by z = 0. Yet another study using still different ISM and feedback models in GIZMO presented in Zhu & Li (2016), who see a similar 'divergence' between MFM and some SPH flavors.  \nNot coincidentally, this epoch where SPH and MFM diverge corresponds precisely to the formation of the 'hot halo' (where virial-shocked gas has a cooling time longer than the dynamical time and establishes a steady-state atmosphere within the halo; see Kereš et al. 2005). In fact, many studies have shown that much, if not most, of the fuel supply for massive galaxies at later times owes to recycled wind material trapped in the CGM in these hot halos (see Oppenheimer et al. 2010; Davé et al. 2011b; Faucher-Giguère et al. 2011; Christensen et al. 2016; Muratov et al. 2017; AnglésAlcázar et al. 2017a). But this requires following the interaction of multi-phase winds moving sub-sonically or trans-sonically through a pressure-supported medium - precisely the regime where the known SPH errors are most problematic (O'Shea et al. 2005; Agertz et al. 2007; Read et al. 2010; Bauer & Springel 2012; Sijacki et al. 2012; Torrey et al. 2012).  \nFig. 25 shows a visual comparison of the gas morphology and phase structure in the CGM at redshifts z ∼ 1, where the differences between methods are most apparent. In D-SPH - which, we emphasize, is a formulation of SPH known to introduce larger errors compared to our 'FIRE-1 SPH' - cool/cold gas in the CGM is primarily locked into numerically spurious 'blobs.' These are a well-known result of numerical errors in SPH, specifically its inability to capture multi-phase fluid mixing interfaces near the resolution scale (see Agertz et al. 2007; Kereš et al. 2012; Saitoh & Makino 2013; Power et al. 2014; Hu et al. 2014; Few et al. 2016). This causes inflows and outflows to 'shred' into blobs, rather than properly mixing. The FIRE-1 'P-SPH' formulation is specifically designed to remove the specific 'surface tension' error in SPH that allows the blobs to be self-insulating and long-lived. P-SPH therefore reduces, but does not completely eliminate, the spurious 'blobs.' We note that all SPH variants shown include the 'artificial conductivity' (entropy diffusion term) above, which has been argued to eliminate such spurious structures (see Price 2008; Wadsley et al. 2008; Read & Hayfield 2012); we will discuss this further below.  \nFig. 25 also shows that in SPH, the halo is more extended with more diffuse boundaries. This is also related to well-known SPH issues, specifically grid-scale heating/noise and difficulty shockcapturing, as well as excessive numerical dissipation of sub-sonic and trans-sonic turbulence (see e.g. Bauer & Springel 2012; Sijacki et al. 2012; Kereš et al. 2012). Both Springel (2010) and Hopkins (2015) show that SPH (even with state-of-the-art artificial viscosity switches from Cullen & Dehnen 2010) produces larger velocity noise around cosmological shocks (compared to MFM or moving-mesh methods), damps turbulence strongly below Mach numbers M∼ 1, and requires ∼ 8 -10 inter-particle spacings ( ∼ 30 -50kpc at the virial radius for the low resolution in Fig. 25) to fully capture shock jumps.  \nThe differences are more striking when we examine the z = 0 temperature-density phase diagram in Fig. 26. While both SPH and MFM produce a significant amount of warm ( ∼ 10 5 K) and cool ( ∼ 10 4 K) CGM gas, we see clearly that in D-SPH there is a large amount of cold ( T glyph[lessmuch] 10 4 K) gas which has survived in the CGM being 'protected' in cold, dense lumps, despite having very short physical mixing times (without magnetic fields - not included here - to 'protect' the clouds). Also, in D-SPH, there is almost no hot, intermediate-density gas ( T glyph[greaterorsimilar] 10 6 Kwith nH glyph[greaterorsimilar] 10 -4 cm -3 ), owing to the difficulties of shock-capturing.  \nTogether, these issues suppress cooling in some SPH flavors, especially from hot gas in more massive halos (again consistent with previous studies; see Kereš et al. 2012; Torrey et al. 2012; Hobbs et al. 2013; Zhu & Li 2016), and they make it relatively 'easier' for cold outflows to escape from the galaxy. Perhaps most troubling, however, in the simplest SPH formulations (e.g. D-SPH here), the 'blobs' and associated under-mixing errors do not converge away at higher resolution, but persist at the resolution scale and actually reach larger density contrasts and contain more mass at higher resolution (this is because they owe to sub-zeroth-order errors; see Agertz et al. 2007; Read et al. 2010). This is why our D-SPH runs at MW-mass deviate more dramatically from MFM at high-resolution, and do not appear to be converging. Of course in MFM (as with all Godunov methods) there is almost certainly some numerical over -mixing at low resolution; a key difference is that this converges away with increasing resolution in a numerically well-defined manner (Harten et al. 1983).  \nWestress, however, that the FIRE-1 SPH (used in all published FIRE-1 simulations) dramatically reduces these discrepancies.", "5.4 Comparing SPH 'Flavors': Adding Diffusion Explains Differences in the CGM": "Fig. 22 also shows the results of using the 'FIRE-1 + Stronger Mixing' SPH flavor instead of our default 'FIRE-1 SPH' in the cosmological m12i simulation. With this enhanced mixing, the results agree well with our MFM simulations. In contrast, as noted above, the D-SPH simulations using the 'density-energy' formulation and a smaller kernel (both of which increase the SPH errors on fluid-mixing problems) show the largest disagreement with MFM.  \nThis demonstrates that the MFM-SPH differences, where present, owe to the treatment of fluid mixing, particularly how winds in 'hot halos' do or do not mix and recycle onto the galaxy. This is a challenging numerical problem, and is almost certainly incompletely resolved in any numerical galaxy formation simulation. Therefore it is difficult to say which method is 'more correct' at low resolution - rather we simply urge caution in interpreting these results at present. The major caveat which must be borne in mind for SPH is that, if the answer depends on an arbitrarily adjustable parameter (e.g. the artificial conductivity), and is non-convergent (as we find here), then SPH has unfortunately limited predictive power - we are forced to calibrate the SPH method at each resolution level to calculations with other codes, rather than simply increase the SPH resolution directly and trust that the errors should converge away (for more discussion, see Hobbs et al. 2013; Zhu et al. 2015; Zhu & Li 2016).  \nThis also suggests that the large central Vc 'spikes' we see in some of our low-resolution MFM runs, which are also sensitive to resolution, may be related to the same wind-mixing physics in hot halos.  \nFor additional systematic comparison of 'improved' SPH and MFM methods in recent cosmological hydrodynamic simulations of galaxy formation, we refer readers to Zhu & Li (2016). Although the simulations there use a completely different treatment of feedback, cooling, and star formation, many of the conclusions - most importantly regarding the effects of fluid mixing in different methods - are identical.", "6 EFFECTS OF 'ARTIFICIAL PRESSURE' TERMS": "In order to follow ISM structure and star formation selfconsistently, we specifically avoid the use of artificial or effective equations of state in our simulations - i.e. in FIRE-2 we do not adopt, as some previous studies (including FIRE-1) have, an 'artificial pressure floor,' numerically forcing the Truelove et al. (1997) criterion to have some value. In those approaches, an additional pressure of the form P artificial ∼ G ( Nh ρ ) 2 /γ is added to the equations of motion, where h is the resolution (local gravitational softening), ρ the density, and N is approximately the ratio of the Jeans length to h (after the artificial pressure is added) - this artificially moves the thermal Jeans mass to a desired scale (and in Lagrangian codes, forms an adiabat with P artificial ∝ ρ 4 / 3 ).  \nA common misconception is that this 'prevents artificial fragmentation,' but this is not accurate (see Teyssier 2015); in cases where P artificial dominates, the medium is by definition Jeans unstable, so fragmentation (usually down to unresolved scales) is physically correct. Rather, what the artificial pressure does is move the fragmentation scale up to a larger scale, set by ∼ Nh , and suppress fragmentation on smaller scales. This is numerically convenient in some circumstances. However, given our star formation prescription, it is not appropriate to suppress fragmentation. Indeed, a cold (weakly pressure-supported), Jeans-unstable, self-gravitating region, with a Jeans scale that is small (stellar-mass scale), is pre-  \nFigure 27. Effects of cooling physics details (§ 7) and artificial 'pressure floors' (§ 6) on a restart of our MW-mass m12i simulation near z = 0, as Fig. 10. Our default model is full-physics, with no 'artificial pressure.' (1) Toy Model Chemistry/Cooling: We replace our default, detailed cooling physics with the toy model in § 7 which puts all gas on a single, solarmetallicity cooling curve, artificially turns off self-shielding (preventing cooling to T glyph[lessmuch] 10 4 K), and removes the requirement that SF be in molecular gas. Because the cooling time is much faster than the dynamical time in all cases, details of the cooling functions have almost no effect on SF. (2) Artificial Pressure Floor: We add an artificial numerical 'pressure floor' as described in § 6, with 'modest' values designed to artificially inflate the thermal Jeans length to always be equal to four softening lengths. For these values this produces some spurious artifacts in e.g. the GMC and star cluster mass function, but does not change any galaxy-scale results, since supersonic turbulence and resolved collapse from much larger scales dominate SF. (3) Extreme Artificial Pressure: We now inflate the artificial pressure so the thermal Jeans length is always > 20 softening lengths. This is far larger than reasonable and makes the entire disk thermally-supported with Q glyph[greatermuch] 1, shutting down substructure. Star formation only occurs when the disk becomes globally unstable leading to 'mini-bursts' - the dynamics are clearly unphysical.  \n<!-- image -->  \nsely that which should be identified as star forming by our algorithm. The un-resolved fragmentation to small scales is already handled, in our method, by assigning the mass to a 'star particle' which represents an aggregate of stars (formed via the un-resolved fragmentation to solar-mass scales). 28  \nFurthermore, 'artificial pressure' prescriptions require careful treatment, as they: (a) can violate energy conservation (providing an infinite source of ' PdV work'); (b) are designed for thermallysupported disks, and it is not clear how to generalize them when turbulent or magnetic pressure dominates (see Myers et al. 2013; Lee et al. 2014); and (c) are not necessary given our numerical methods.  \n28 We note that Truelove et al. (1997) (or Bate & Burkert 1997) did not propose an 'artificial pressure.' Rather, their conclusion was that in order to explicitly resolve fragmentation down to the thermal Jeans mass scale (the brown dwarf scale, ∼ 0 . 1 M glyph[circledot] , in the cold ISM), one needed to maintain glyph[greaterorsimilar] 4 resolution elements per Jeans mass. They suggested this as a local refinement criterion, for AMR-type simulations of proto-star formation. Later, Federrath et al. (2010) and many subsequent authors (e.g. Padoan & Nordlund 2011; Myers et al. 2013; Gong & Ostriker 2013; Federrath 2015; Skinner & Ostriker 2015) pointed out that one can instead use this criterion to determine when local fragmentation is unresolved, and assign the mass to sink particles instead of following it explicitly (standard practice in GMC-scale simulations) - this is, of course, our star formation prescription.  \nSpecifically Kratter & Lodato (2016); Nelson (2006), and Chiaki & Yoshida (2015) show (with very different approaches) that in Lagrangian methods, as long as (1) the Toomre mass is resolved (our own criterion in § 4.1.3), (2) the gas disk scale-height is resolved (identical to the Toomre criterion if Q ∼ 1), and (3) fully-adaptive gravitational softenings are used for gas ( glyph[epsilon1] gas = hi ), the fragmentation cascade converges accurately and is numerically stable (resolution will truncate the lower-limit of the cascade, but all the power is on large scales, just as in turbulent motions). This appears to be confirmed by our tests in Figs. 10 & 13.  \nThat said, because our FIRE-1 simulations used such a prescription, we have checked whether adding an artificial pressure term in the simulations makes a qualitative difference to our results; in Fig. 27, we implement the standard prescription above with N = 4. For N ∼ 1 -5, the difference in bulk galaxy properties (SFR, sizes, stellar masses) are small. This is because, for N = 4, the artificial pressure will stabilize an otherwise unstable spherical overdensity with radius R glyph[lessorsimilar] 30pc m 1 / 3 i , 1000 ( n / 10cm -3 ) -1 / 3 ; so for the resolution and typical gas densities inside the central ∼ 4 kpc (the half-light radius) in Fig. 27, the largest GMCs can still collapse. Manuel et al. (2016) reach similar conclusions (the added terms have little effect) in GMC-scale simulations. However, we caution that there are some potentially un-physical small-scale artifacts introduced by an artificial pressure term, including: noise in the lowtemperature gas thermal properties, a 'bump' in the mass function of GMCs/self-gravitating clouds (owing to 'pileup' of clouds at the scale where the fragmentation cascade is truncated), and spurious stellar clusters (see Nelson 2006; Lukat & Banerjee 2016). If we adopt N glyph[greaterorsimilar] 10 (at the resolution shown here), then we are, by construction, forcing the gas to be smooth on scales larger than the gas disk scale height - i.e. we suppress any substructure in the galactic disks, which in turn artificially suppresses star formation everywhere except at galaxy centers.", '7.1 Cooling, Heating, and Self-Shielding Physics': "We now consider how the details of cooling alter our predictions. Of course, it goes without saying that cooling is critical for galaxy and star formation - 'turning off cooling' would produce no galaxies! Moreover, detailed chemistry and cooling physics are obviously critical for some observables such as line diagnostics (e.g. fine-structure or CO excitation; see Richings & Schaye 2016), but it is not clear whether this in turn produces any dynamical effects on bulk galaxy properties studied here. Our interest here is therefore in exploring whether the more subtle details of cooling, which are often uncertain at the factor ∼ 2 level and treated differently or incompletely in different codes (for example, some include chemical networks for low-temperature and/or non-equilibrium cooling, some include only pre-tabulated cooling tables, some include additional sub-dominant cooling processes), ultimately have significant effects on our galaxy-scale predictions (masses, SFHs, morphologies, sizes, etc.).  \nFigs. 27-28 show the results of turning on or off different portions of the cooling physics in our simulations. 29 Specifically, we  \n1  \nFigure 28. Effects of the details of cooling & chemistry (as Fig. 27), on cosmological simulations of a dwarf ( m10q ) and MW-mass ( m12i ) galaxy, as Fig. 4. We compare our default, full-physics & chemistry treatment, to the extremely simplified 'toy' cooling model from Fig. 27 & § 7 (single, solar metallicity Z = Z glyph[circledot] cooling curve, with no self-shielding so no lowtemperature cooling or chemistry). We also compare a version of the toy model assuming Z = 0 (primordial abundances only, no metal-line cooling at low or high temperatures). The simplified cooling model produces remarkably similar galaxies to the full-physics model, provided the metallicity is similar ( Z = 0 in m10q and Z = Z glyph[circledot] in m12i ). Because cooling times within galaxies are often much shorter than dynamical times so as long as some cooling channel exists, the details of cooling have weak effects. The stellar masses are higher with the Z = Z glyph[circledot] toy model, owing to solar metallicity being assumed in all gas - this is higher than the metallicity of large-scale super-bubbles in m10q or outer parts of the extended 'hot halo' of m12i , artificially enhancing cooling at T ∼ 10 5 -10 7 K. Conversely, with Z = 0 cooling only, SNe bubbles and the hot halo cool less efficiently; notably in m12i this reduces the SFR around z ∼ 1, where the dense central bulge forms, flattening the rotation curve substantially. Note that taking our 'full cooling' model and only disabling high-temperature ( T > 10 4 K) metal-line cooling produces nearly identical results here to the 'toy' Z = 0 model.  \n<!-- image -->  \ncompare our default FIRE-2 cooling physics described in § 2.3 to simulations adopting an extremely simplified toy model. In the toy model, we place all gas on a single, invariant cooling curve Λ( T , Z , I ν , n ) → Λ( T ) for fully ionized, non-self-shielded gas (so there is essentially no cooling to T glyph[lessmuch] 10 4 K), with a single metallicity Z = Z glyph[circledot] or Z = 0 (primordial cooling only). We also remove the requirement that stars form in molecular/self-shielded gas, since this is not self-consistently computed in the toy model. This toy model is not intended to be physically correct. Rather, it is intended only to illustrate how even radical changes to cool-  \nFigure 29. Temperature-density phase diagram of all gas in the simulation box at z = 0, for the m12i simulations in Fig. 28 with our default, full cooling physics implementation and the 'toy model' which puts all gas on a single cooling curve, with no self-shielding, at a uniform metallicity Z = Z glyph[circledot] or Z = 0. As expected, without self-shielding, the toy models produce no cold gas ( T glyph[lessmuch] 10 4 K). However, this is still 'cold' as far as providing dynamical support to the disk (thermal Q < 1) is concerned, and t cool glyph[lessmuch] t dyn remains true, and so even at ∼ 10 4 Kgas can and should physically fragment efficiently. Therefore, this produces little or no dynamical effect in the simulations. The Z = Z glyph[circledot] run noticeable suppresses the hot gas at T ∼ 10 6 Kand densities n ∼ 10 -4 -10 -2 cm -3 - much of this is gas in the extended CGM and nearby IGM (around R vir) which has been shocked (by SNe & O-star/AGB winds but also the halo virial shock), and is highly sub-solar in metallicity, so the toy model over-estimates the cooling rates by a factor ∼ 10 -30. However, this does not contribute much to the galaxy accretion rate. The 'primordial' ( Z = 0) model, on the other hand, is closer to correct in the hot halo at ∼ R vir, but significantly under-estimates the cooling rate for SNe bubbles in the gas disk and inner halo (while the shape of the warm/hot gas distribution is similar, there is much less gas in the warm phase compared to the Default cooling run).  \n<!-- image -->  \nysics (much more extreme than typical model differences in cooling) have relatively modest impacts on galaxy formation.  \nFor dwarf galaxies, with either Z = 0 or Z = Z glyph[circledot] (Fig. 28), we see there is no large effect on the galaxy formation history. There is some enhanced early cooling and star formation at Z = Z glyph[circledot] as expected (this is much higher than the mean metallicity in m10q ). For MW-mass galaxies, we again find only modest changes to the dynamics in both full cosmological simulations (Fig. 28) or restarts near z = 0 with fixed initial conditions (Fig. 27), if we adopt approximately the 'correct' metallicity in and around the galaxy, Z = Z glyph[circledot] . However for MW-mass systems, the primordial coolingonly ( Z = 0) run does differ more dramatically. We have also examined the visual galaxy morphology, and wind outflow rates and covering factors, in each run and find similar results.  \nFig. 29 shows the z = 0 temperature-density phase diagram of the m12i runs: clearly, details of the phase structure differ, as expected. With the toy model, there is no dense, cold gas ( T glyph[lessmuch] 10 4 K), by construction. For the toy model with Z = Z glyph[circledot] , the hot, denser gas has also cooled more efficiently - much of this was gas in the halo with lower metallicities which now radiates more rapidly than it should. But this does not appear to have a strong dynamical effect (there is some increase in the late-time SFR, stellar mass, and circular velocity evident in Fig. 28 in this run, owing to additional cooling, but the effect is weak).  \nThis weak dependence of galaxy properties on the more subtle details of cooling physics at T glyph[lessorsimilar] 10 4 K (for approximately the correct metallicity) is consistent with previous studies. For example, both Hopkins et al. (2011) and Hu et al. (2016) added/removed photo-electric heating, and Hopkins et al. (2012c) added/removed molecular cooling at low temperatures, and all found this had essentially no effect on SFRs, wind outflow rates, and other galaxy properties. A more detailed study on sub-galactic scales by Glover & Clark (2012) considering each microphysical cooling mechanism in turn showed that as long as some cooling channel exists, the details of that cooling are largely irrelevant in the dense ISM.  \nThe reason for this is simple: for almost all gas (by mass) in galaxies, the cooling time is much shorter than the dynamical time  \n( t cool glyph[lessmuch] t dyn), and the gas is super-sonically turbulent. Under these conditions, the details of the cooling physics have little effect on dynamics - cooling is not a 'rate-limiting' step. We find our predicted SFRs are similar so long as the cooling model captures two key effects: (1) t cool glyph[lessorsimilar] t dyn for gas with T glyph[lessmuch] 10 7 K ( n / 0 . 01cm -3 ) (although for SNe blastwaves, the more relevant comparison is to the remnant expansion time, giving T glyph[lessorsimilar] 10 6 K), and (2) gas in this rapidly-cooling regime has equilibrium temperatures glyph[lessorsimilar] 10 4 K or colder (the colder gas has no significant thermal pressure support, so assigning it T ∼ 10 4 does not appreciably change its large-scale dynamics).  \nBasically, in gas with temperatures glyph[lessorsimilar] 10 6 -10 7 K, details of the cooling physics are not important. Above this threshold, t cool glyph[greaterorsimilar] t dyn and cooling can be important; this can matter for (1) super-bubbles and hot galactic winds, and (2) 'hot halos' of virial-shocked gas around massive galaxies. For winds (1), thermal pressure-driven winds generally cool adiabatically as they expand (Chevalier & Gardner 1974), so the exact radiative cooling is not important so long as the transition from t cool glyph[lessmuch] t dyn to t cool glyph[greatermuch] t dyn happens at more or less the correct temperature (see Paper II). In hot halos (2), the cooling rate in the halo center determines accretion rates onto galaxies. This is irrelevant in dwarfs ( T vir glyph[lessmuch] 10 7 K), but important in MW-mass halos, hence we see significant effects in m12i if we adopt the Z = 0 toy model. We have re-run this turning on/off each piece of cooling physics in turn, and verified that almost the entire change in Fig. 28 owes to the Z = 0 toy model having no high-temperature ( T > 10 4 K) metal-line cooling. This suppresses the SFR around z ∼ 1 -2, as the hot halo forms. In the 'default cooling' run, accretion onto the galaxy at this time is dominated by recycled wind material from earlier episodes, mixed with other less metal-rich CGM material (leading to additional cooling in that gas) - the additional cooling from metal-rich wind material is absent in the toy model. Similar conclusions have also been reached in previous studies (e.g. Choi & Nagamine 2009; Piontek &Steinmetz 2009; Schaye et al. 2010).  \nFigure 30. Comparison of predicted galaxy formation histories as Fig. 4, with and without an explicit numerical 'metal diffusion' term. Our 'core physics only' runs using the MFM hydrodynamic solver follow fixed-mass elements (so there is no advection between elements and metals are strictly 'locked' to the element they were injected into via stellar mass-loss). The 'metal diffusion' runs add an explicit passive-scalar diffusion term to account for sub-resolution diffusion and turbulent mixing of metals between neighboring elements (see Appendix F3). This is flux-limited so that the diffusivity is never larger than would be obtained from simply numerical diffusivity in our MFV (moving-mesh or AMR-like) hydrodynamic solver. There is no systematic effect at any mass scale. In addition to the runs shown we have also compared m09 and m12b,c,f,q and find the same. Su et al. (2017) examine non-cosmological disks and our m10q and m12i runs and show the same, also for the gas and stellar morphology, ISM density distribution in different phases, galactic wind mass-loadings and velocity distributions, and ISM turbulent velocities.  \n<!-- image -->", '7.2 Effects of Numerical Metal Mixing': "In § 2.6 we note that, in our 'core physics only' simulations, using our default MFM hydrodynamic solver, passive scalars such as metals are 'locked' to the gas resolution elements into which they are initially deposited (e.g. the nearest few elements to each SN). This occurs because the MFM method is a finite-mass method, i.e. advective fluxes vanish identically . But since we have finite resolution, this fails to capture e.g. small-scale turbulent eddies and microphysical diffusion of metals between the boundaries of two neighboring resolution elements. In finite-volume methods such as grid-based Eulerian codes, or moving-mesh codes, or the meshless finite-volume 'MFV' method in GIZMO (which is similar to MFM, differing only in how the face is assumed to deform when solving the Riemann problem between cells) there is an explicit advective (mass) flux, hence there is always some un-avoidable mixing at the boundary (in fact, it is well known that such methods will tend to over-mix, relative to a converged solution, at finite resolution).  \nUnder some circumstances, it may be desirable to explicitly model this microphysical metal transport between resolution elements. In Appendix F3 we describe how this is implemented as an 'additional physics' option in FIRE-2. Essentially this follows the Smagorinsky (1963) approximation, where we explicitly solve a diffusion equation between cells with the diffusivity κ turb set to the resolution-scale 'eddy diffusivity' ( ∼ glyph[lscript] | ∆ v turb ( glyph[lscript] ) | ; i.e. assuming the mixing time scales with the eddy turnover time, based on the local spatial resolution and velocity field). This is tested directly in  \nGIZMO with our same MFM solver in idealized, converged 'turbulent box' simulations in Colbrook et al. (2017) and Rennehan et al. (in prep.). An upper limit to the mixing between cells is enforced, equal to what would be obtained by the GIZMO MFV solver for the same cells and same timestep - in other words, the maximum possible diffusivity of the added term in MFM is simply the un-avoidable diffusivity inherent to the MFV method (itself very similar to that in moving-mesh codes). This allows us to study how this particular numerical diffusion term, which is always present in hydrodynamic methods with explicit mass fluxes, alters our predictions.  \nFig. 30 shows the effects of adding this term, in several of our simulations at both dwarf and MW mass scales, on the galaxy formation history and z = 0 mass profiles. Fig. 31 shows the effects on the z = 0 ISM and CGM gas phase distribution and visual morphologies in a MW-mass galaxy (there is no detectable difference in dwarfs, consistent with other properties we have surveyed). We see essentially no systematic effect in any gross galaxy property (star formation history, stellar mass, mass profile, rotation curves, visual morphology of gas or stars), including the mean metallicity itself (both in gas and stars). This is true at all mass scales we have explored and at all resolution levels. There is slightly more cool gas at z = 0 in the disk of the m12f simulation with metal diffusion, which could owe to metal mixing promoting additional cooling from the hot halo, but we caution that this effect is similar in magnitude to stochastic variations between runs. It is worth noting that the cold gas within the disk ( T glyph[lessorsimilar] 10 4 K, nH glyph[greaterorsimilar] 0 . 1 cm -3 ) forms a slightly  \nFig. 32 further demonstrates this by considering variations to the standard metal-mixing term, at both dwarf and MW mass  \n<!-- image -->  \nFigure 31. Top: Visual morphology (as Fig. 9), at z = 0, for the m12f ( mi , 1000 = 56) simulations run with our 'core physics only' ( left ) or core physics with the additional metal diffusion term as Fig. 30 ( right ). Bottom: Phase diagram of the ISM & CGM (as Fig. 26) for the same. The diagrams are qualitatively very similar but there is slightly more cool gas in the ISM at z = 0 in the runs with metal diffusion, owing to its allowing marginally enhanced cooling from the hot halo that forms at late times in massive galaxies. There is no detectable difference at dwarf masses.  \n<!-- image -->  \ntighter (more visually obvious) cooling sequence; this is because the metal diffusion smooths local particle-to-particle variations in gas-phase metallicity in e.g. the same GMCs, and the metals are the dominant coolant in this regime.  \nSu et al. (2017) explore in greater detail the effects of this metal-mixing term in the ISM of both isolated (non-cosmological) disks and cosmological simulations with the FIRE-2 physics. They specifically compare the bivariate distribution of gas temperatures and densities in the ISM (at several redshifts); galactic wind outflow rates, phases, density distributions, and velocity distributions; turbulent velocity distribution functions in the ISM; and star formation histories. They conclude in all cases that the metal-mixing terms have weak or negligible effects. In Escala et al. (2018) we present the actual abundance ratio distribution functions within the galaxies here; there, it is clear that the mixing terms do have important effects (as expected). However because the mixing terms conserve (by construction) the mean metallicity and total metal mass, and we have already shown that the metallicity enters only relatively weakly into dynamical effects via its effect on the cooling rates (see § 7), the dynamical effect of these mixing terms on other galaxy properties is weak.  \nFigure 32. Comparison of the effects of explicit numerical metal diffusion, as Fig. 30, at dwarf ( m10q ) and MW ( m12i ) mass scales, as a function of the diffusivity. We compare: (1) Core Physics Only : metals strictly locked to original injection element. (2) Metal Diffusion : the default numerical implementation in Appendix F3. (3) Metal Diffusion x10 : as 'Metal Diffusion' but with the numerical diffusion coefficient increased by a factor 10. (4) Metal Diffusion x30 : as 'Metal Diffusion' but with the numerical diffusion coefficient increased by a factor 30. (5) Fixed Metal Diffusivity : replacing the adaptive numerical diffusion prescription with a pure isotropic diffusion equation of fixed diffusivity κ turb = 10 24 cm 2 s -1 , comparable to the largescale eddy diffusivity for ISM turbulence with σ ∼ 10kms -1 and driving scale ∼ 100pc. (6) MFV Metal Fluxes : explicitly adding the cell-to-cell metal advection fluxes that would appear if we solved the hydrodynamic equations with GIZMO's 'meshless finite volume' (MFV) method, instead of our default MFM method (which has vanishing advective fluxes). This is approximately the 'numerical metal diffusion' that would un-avoidably appear in a moving-mesh type code. In all cases, the effects are weak.  \n<!-- image -->  \nscales. Specifically we consider (in addition to the 'core physics only' and 'metal diffusion' runs) cases where we add the metal diffusion term but arbitrarily multiply the diffusivity κ turb by an additional factor of 10 or 30. We also consider a case where we replace the standard adaptive eddy diffusivity (defined in Appendix F3) with a constant κ turb = 10 24 cm 2 s -1 (i.e. simply treat the metals as if they obey a strict, constant-diffusivity diffusion equation, and we also remove the MFV-based limiter above). This value of κ turb corresponds crudely to typical values of the turbulent diffusivity estimated by Eq. F5 in Appendix F3 on the maximal turbulent scales (for glyph[lscript] ∼ 100pc and | ∆ v turb ( glyph[lscript] ) | ∼ 10kms -1 ). Finally, we consider a case where we include no diffusion equation, but instead we solve the metal flux exactly as we would if we used the MFV hydrodynamic method (with explicit mass fluxes), and include the resulting metal fluxes (even though all other fluxes use the MFM solution, identical to our 'standard' runs). In other words we always transfer  \nthe metals that would be 'numerically diffused' in a finite-volume code. In each of these cases, we see no effect on the properties in Fig. 32.", '7.3 Yields and Explicit Abundance-Ratio Tracking': 'As noted above, we track 11 species on-the-fly with yield tables for SNe (Ia & II) and stellar mass-loss rates given in Appendix A. But nucleosynthetic yields (even IMF-averaged) have very large uncertainties, especially at progenitor metallicities far from solar. For example, the Woosley & Weaver (1995) yields differ from those in Nomoto et al. (2006) even at Z = Z glyph[circledot] by ∼ 0 . 4dex for Mg and Ne. However, given the weak dependence of our results on the details of the cooling curve shape (especially at low temperatures; § 7.1) and detailed metal mixing (§ 7.2), we do not expect this to have large dynamical effects on other galaxy properties. Moreover, even high-temperature metal-line cooling (which is important to include) is dominated by O, which constitutes most of the metal mass and therefore is better constrained (differing by < 10% between the Woosley & Weaver 1995 and Nomoto et al. 2006 models, for example). We have verified this directly by re-running our m10q ( mi , 1000 = 0 . 25) and m12i ( mi , 1000 = 56) simulations, replacing the default Nomoto et al. (2006) yields for core-collapse SNe with those from Woosley & Weaver (1995). In all properties surveyed here, we see no detectable difference.  \nGiven this and our comparisons in § 7.1, it appears that as long as the total metallicity is approximately correct, differences in the detailed abundance ratios introduce no major dynamical effects. We have therefore re-run the same simulations, ignoring all detailed abundance information and simply tracking the total metal abundance, then assuming solar abundance ratios in the cooling and other relevant routines. We find this produces nearly identical results to our default simulations (in which the cooling is explicitly solved species-by-species).  \nTo leading order, then, detailed abundance patterns are important as tracers for predictions of various observables, but do not introduce important dynamical effects.', '8 EFFECTS OF THE STAR FORMATION ALGORITHM': "It has been extensively demonstrated in the literature that, provided the Toomre scale is resolved and stellar feedback is treated explicitly (as in our simulations), the exact resolution-scale SF model has essentially no effect on predicted galaxy-scale SFRs. For the sake of completeness, we demonstrate this here, but refer readers to Hopkins et al. (2011, 2013d,e) and Kim et al. (2013) for more extensive discussions (see also § 4.1.3).  \nFigs. 33-34 show the effects of varying the resolution-scale star formation model in the simulations, in both full cosmological runs at different mass scales and resolution, and re-starts of our MW-mass simulation at late times (guaranteeing an identical initial condition for the comparison). We compare our default model from § 2.4 (requiring gas be self-gravitating, self-shielding/molecular, Jeans unstable, and have density n > n crit with n crit = 1000cm -3 ), to variations with n crit = 5 -1000cm -3 ; turning on/off the self-gravity (virial), molecular, and Jeans criteria; and arbitrarily multiplying the SFR per free-fall time ˙ ρ ∗ in the gas which meets these criteria by a constant factor glyph[epsilon1] sf = 0 . 01 -100.  \nIn every case, the predicted SFR and all other galaxy properties we examine here are essentially identical. This has now been shown in many different contexts, with simulations using different detailed implementations of star formation and stellar feedback; mass resolution ranging from sub-solar to ∼ 10 6 M glyph[circledot] ; isolated (non-cosmological), galaxy-merger, and fully-cosmological  \nsimulations; circum-nuclear simulations of star formation around AGN; and simulations using very different feedback mechanisms (removing SNe, stellar mass-loss, or HII photo-heating); and in different codes (see e.g. Saitoh et al. 2008; Shetty & Ostriker 2008; Hopkins et al. 2012b, 2013c, 2011, 2013a,d,e, 2016; Federrath & Klessen 2012; Agertz et al. 2013; Kim et al. 2011, 2013; Kim & Ostriker 2015). We specifically showed the same was true in our FIRE-1 simulations in Hopkins et al. (2014); Orr et al. (2018). This is also consistent with simulations of individual GMCs (which track star formation on much smaller scales), and analytic 'multifree-fall' models of star formation, in which independent clumps collapse on their own local free-fall times (see e.g. Federrath & Klessen 2012, 2013; Faucher-Giguère et al. 2013; Guszejnov & Hopkins 2016; Grudi'c et al. 2018).  \nThis is because the rate-limiting step in star formation does not occur at the resolution limit of the simulations; rather it is the formation and efficient destruction of the largest self-gravitating objects (large GMC complexes, which have dynamical times ∼ 100Myr). Because these are resolved , and our force softening is fully-adaptive (and we do not force some artificial lower-limit to the softening), a self-gravitating cloud or sub-clump will continue to fragment, as it should, on its local free-fall time, until it either forms stars or is destroyed by feedback (see below). This means that absent star formation or sufficient feedback, the densities in such a clump will (correctly) become arbitrarily high and the internal dynamical times will become arbitrarily short, within a finite physical time (of order the original parent cloud free-fall time). So whatever value of n crit we set will eventually be exceeded, and even if glyph[epsilon1] sf glyph[lessmuch] 1 (so the cloud collapses faster than it forms stars), the actual physical SF timescale will become arbitrarily short within the clump as it collapses. Provided star formation can occur, then the simulations above have shown that feedback self-regulates the level of star formation. Stars form until sufficient numbers of young stars are present that their feedback destroys the parent clouds, which may occur after just a small fraction of the cloud turns into stars (e.g. Grudi'c et al. 2018), giving a low time-averaged efficiency per freefall time. This self-regulating behavior depends on the strength of feedback - changing the strength of feedback immediately changes the equilibrium SFR in a galactic disk - but is independent of how individual stars form in dense gas.  \nOf course, the ratio of SFR to dense gas at the resolution limit depends, by construction, on the SF model; comparison to dense gas tracers within galaxies (e.g. HCN as opposed to CO) in previous work appears to favor the 'default' normalization and SF criteria we adopt here (see the comparison in Hopkins et al. 2013e). But this is only relevant within the densest resolved gas clumps.", '9 STELLAR FEEDBACK': 'Many previous studies have argued that feedback is the most important determinant of star formation in galaxies (see references in § 1). In this section we therefore explore how basic feedback physics changes our predictions; in Paper II-Paper III, we examine more subtle physical and numerical aspects of the implementation of said feedback.', "9.1 Feedback 'Strength'": "In Fig. 35, we simply re-start our m12i simulation near z = 0 as Fig. 10, but arbitrarily multiply all feedback rates (e.g. SNe rates, wind kinetic luminosities and mass loss rates, stellar luminosities) by a factor of ∼ 3 ('strong') or ∼ 1 / 3 ('weak'), relative to the  \nFigure 33. Comparison of galaxy properties in cosmological FIRE-2 simulations (as Fig. 4) as a function of the resolution-scale assumptions for individual star particle formation. We compare m09 , m10v , m10q , and m12i (Table 1). For each, we compare a run with a minimum density for star formation of n crit = 10 or n crit = 1000cm -3 ; recall, in all cases the SF is still only allowed in self-gravitating, self-shielding gas. We have repeated this experiment for m10q and m12i at lower resolution, and for m11q , m11v , m10y with n crit = 100 and = 1000, and reach the same conclusion in every case. For m10v and m10q , we also compare a simulation in which we arbitrarily multiply the SFR in the 'eligible' gas (sufficiently dense, self-shielding/molecular, locally self-gravitating) by a factor glyph[epsilon1] SF = 0 . 01 (i.e. artificially 'slow down' the SF in the collapsing gas by a factor of ∼ 100, relative to its dynamical time). Consistent with extensive studies in the FIRE-1 and previous simulations, these choices have no effect on our predictions.  \n<!-- image -->  \n'default' values (which, recall, are taken directly from stellar evolution models without adjustment). We stress that there is no physical motivation for this - this is a larger shift than most uncertainties in stellar evolution - but we consider it purely for illustrative purposes. Exactly as expected in 'feedback-regulated' scenarios, this produces a corresponding direct shift in the SFR - the SFR is selfregulating at the level where feedback offsets gravitational collapse, so if feedback is 3 × stronger 'per star,' then 3 × fewer stars form when the equilibrium is realized. This is almost exactly what we see (the shift in SFR is very slightly sub-linear, owing to non-linear effects). This is consistent with many previous simulation and analytic studies (see Shetty & Ostriker 2008; Kim et al. 2011; Hopkins et al. 2011, 2012c; Faucher-Giguère et al. 2013; Orr et al. 2018, and references therein).", '9.2 The Role of Different Feedback Mechanisms: Galactic Winds & Cosmological Timescales': "In Fig. 36 we consider how turning off each feedback mechanism in turn alters galaxy evolution in fully cosmological simulations. In both dwarf and MW-mass galaxies, we clearly see SNe have the most dramatic effect. This is because what ultimately regulates galactic star formation efficiencies over a Hubble time is the  \ncompetition between inflow and galactic outflows, and the highspeed outflows that are not simply recycled quickly are predominantly SNe-driven (e.g. Davé et al. 2011a, and references therein). Without SNe, galaxies form far too many stars (factor > 100 in our m10q dwarf; factor ∼ 2 in the MW-mass system, going from ∼ 30% to ∼ 60% baryon-to-star 'conversion' efficiency); the stars form too early (everything collapses rapidly into early dwarf halos, and has mostly turned into stars by redshift z ∼ 3 -4); the resulting galaxy is too dense (rotation curves peak strongly, even in small dwarfs, and are never 'gently rising'); and the metallicities are essentially constant across stellar mass because galaxies have similar formation histories and baryon conversion efficiencies (latetime metallicities are lower for massive galaxies because there is no SNe enrichment, only RSG/AGB enrichment). 30 Likewise properties of the ISM (gas phases) and CGM (covering factors of outflows and metals, especially) are grossly discrepant with observations. We also note that because all of our galaxies are star-forming, the  \nFigure 34. Effects of the resolution-scale assumptions for individual star formation, in re-starts of a MW-mass galaxy ( m12i ) at late times, as Fig. 10, keeping fixed mass resolution mi , 1000 = 56 but changing the SF algorithm. In our 'default' model gas which is self-gravitating (sub-virial), selfshielding/molecular, dense ( n > n crit = 1000cm -3 ), and Jeans unstable forms stars at a rate ˙ ρ ∗ = ρ mol / t ff. We compare: (1) removing the selfgravity & molecular restrictions (all dense gas can form stars), (2-3) lowering/raising n crit = 5 -1000cm -3 , and (4-5) multiplying ˙ ρ ∗ by an arbitrary factor glyph[epsilon1] sf = 0 . 01 -100. As in Fig. 33 and our previous studies, this has no effect on the galaxy-scale SFR, because it is feedback-regulated.  \n<!-- image -->  \nSNe Type II rate dominates over the Ia rate at all times. All of these results are consistent with many previous studies (see references in § 1). We discuss the effects of SNe in much greater detail in Paper II - our purpose here is only to confirm they are indeed critical.  \nTurning off continuous stellar mass-loss (OB/AGB winds), the most obvious effect is that the late-time SFRs of both dwarfs and MW-mass galaxies are suppressed , relative to their 'default' values. Despite the fact that OB winds carry (roughly) comparable energy and momentum flux to SNe, and AGB winds can do the same if the relative star-gas velocity is sufficiently large (see Appendix D), their net effect seems to be increasing the supply of gas mass within the galaxy which can eventually cool and form stars at late times. This is also consistent with previous studies, many of which have noted that OB winds tend to not be sufficiently well-confined to build up comparable momentum (after coupling) to SNe (Jungwiert et al. 2001; Ciotti & Ostriker 2001; Ciotti et al. 2009; Schaye et al. 2010; Leitner & Kravtsov 2011; Novak et al. 2011; Gan et al. 2014; Choi et al. 2015).  \nRemoving all radiative feedback - including the effects of the meta-galactic UV background - has a dramatic effect on dwarf stellar masses, even at halo masses glyph[greaterorsimilar] 10 10 M glyph[circledot] with V max glyph[greaterorsimilar] 40kms -1 that should be well above the threshold for UVB 'quenching' star formation. This will be explored in more detail in future work (Wheeler et al., in prep.). In short, the stellar masses increase by a factor ∼ 10 (consistent with our FIRE-1 results; Hopkins et al. 2014). Unlike removing SNe, the qualitative shape of the SFH is similar to the 'full feedback' case - the stars do not form uniformly early, and the SFH is still 'bursty,' but the bursts reach much larger amplitude, owing to the much larger gas supply. Removing only local radiative feedback (keeping the UVB, but removing radiative feedback from explicit star particles) has more subtle effects: the SF becomes much more 'bursty,' the stellar mass changes but by a smaller factor, and the metallicity is suppressed; we discuss these in detail in Paper III. In MW-mass halos, the UVB has weak effects,  \nFigure 35. Effects of arbitrarily changing the strength of stellar feedback on a re-start of an identical m12i simulation (with mi , 1000 = 56) at low redshift (i.e. guaranteeing identical initial galaxy properties), as Fig. 10. Here, we compare our 'default' model (in which all feedback physics, rates and energetics are taken from stellar evolution models without re-adjustment), and compare it to two models where we arbitrarily multiply all feedback rates (e.g. SNe rates, wind kinetic luminosities and mass loss rates, stellar luminosities at all frequencies) by a factor of ∼ 3 ('strong') or ∼ 1 / 3 ('weak'). These are much larger than actual physical uncertainties in these quantities, we simply show it for illustrative purposes. Clearly SF is instantaneously feedback-regulated: the steady-state SFR is inversely proportional to the feedback strength, as expected.  \n<!-- image -->  \nbut removing local radiative feedback produces a clear, significant increase in the central Vc (perhaps expected, as the galaxy center is exactly the region where we expect massive, dense GMCs which must be pre-processed by radiative feedback in order to reach low densities where subsequent SNe can efficiently expel the material, see Thompson et al. 2005; Murray et al. 2005; Hopkins et al. 2011). The dependence is similar in FIRE-1, although there we saw a stronger dependence of stellar mass on radiative feedback; the difference owes to the improvements in the treatment of SNe which allow, in massive GMCs, for the first SNe to play a similar role to radiative feedback and 'pre-process' the cloud for subsequent SNe.  \nAs Fig. 36 clearly demonstrates, our conclusions about the systematic effects of different feedback mechanisms are independent of the hydrodynamic method (MFM vs. SPH), even where the methods themselves produce significant systematic differences. Likewise, we have checked the conclusions above are independent of resolution, star formation prescription, and details of the cooling functions.", '9.3 The Role of Different Feedback Mechanisms: Self-Regulation Within Galaxies and the Kennicutt-Schmidt Relation': "In Fig. 37 we repeat the experiment from Fig. 36 turning on and off different feedback mechanisms, but in restarts of the z ∼ 0 MWmass galaxy as Fig. 10. This allows us to separate non-linear, longtimescale cosmological effects on the SFR (e.g. galactic wind generation and recycling), from the instantaneous self-regulation of star formation within a galaxy.  \nAs expected, removing SNe produces a systematically higher SFR. Removing OB/AGB winds leads, initially, to higher SFRs, as the additional feedback (gas heating via shocked winds) is no longer present (but the gas mass of the disk is still essentially fixed at its initial value); however the SFR then declines as the miss-  \n1  \nFigure 36. Effects of different stellar feedback physics on cosmological galaxy formation histories, as Fig. 4. We compare a dwarf ( m10q ) and MW-mass galaxy ( m12i ), both run with our default, higher-order accurate hydrodynamic solver (MFM; left ); but we also show a set of simulations of m12i with the same physics variations using the smoothed-particle hydrodynamics (SPH; right ) method, to demonstrate that even though the hydrodynamic solvers produce differences in massive galaxies, the qualitative effects of different feedback mechanisms are identical (independent of the hydrodynamic method). (1) Default: Our default models include all standard stellar evolution processes. (2) No OB/AGB mass-loss: Removing continuous stellar mass loss (both OB and AGB-star winds) produces slightly lower metallicities (owing to the lack of recycling) and significantly lower late-time SFRs - it appears the primary role of stellar mass-loss is to provide an additional source of gas fueling late-time SF in both dwarfs and MW-mass systems. (3) No Radiative feedback: This removes all radiative feedback (radiation pressure as well as photo-ionization and photo-electric heating by local particles and the meta-galactic background). In dwarfs (even with V max ∼ 40kms -1 , shown here), removing the photo-ionization heating (dominated by the UVB) produces ∼ 10 × larger SFRs and stellar masses (producing large bursts that make a core and lower Vc ). In massive galaxies, the effects are weaker but removing radiation pressure produces significantly higher central densities (more strongly-peaked rotation curves in the central ∼ 5kpc). (4) No Supernovae: SNe clearly dominate on cosmological scales, as removing them produces orders-of-magnitude higher SFRs at early times, giving rise to runaway collapse to extremely high densities until the gas is depleted.  \n<!-- image -->  \ning mass-loss is unable to 're-supply' gas lost to star formation. Most interesting, without radiative feedback, the SFR systematically rises with time, comparable to or even exceeding the SFR of the no-SNe test! This is consistent with a number of studies of isolated (non-cosmological) galaxy simulations (see Hopkins et al. 2011, 2012c; Kim et al. 2013; Agertz et al. 2013; Benincasa et al. 2016).  \nWhy does the SFR appear more sensitive to radiative feedback in Fig. 37 and these isolated galaxy simulations, as opposed to Fig. 36? On cosmological spatial and temporal scales, a sim-  \nple steady-state argument (e.g., Bouché et al. 2010; Davé et al. 2012; Feldmann 2013; Lilly et al. 2013; Feldmann 2015; Mitra et al. 2017) implies the sum of the time-average galactic SFR 〈 ˙ M ∗ 〉 and superwind mass-loss rates 〈 ˙ M wind 〉 ≡ η 〈 ˙ M ∗ 〉 (where η is the wind mass-loading, which can depend on arbitrary galaxy properties) should equal the inflow rate ˙ M in (we assume gas mass in the disk is steady-state; see Davé et al. 2011a). Thus 〈 ˙ M ∗ 〉 ∼ ( 1 + η ) -1 ˙ M in. However, on scales within the galactic disk (and timescales of order the galaxy dynamical time, much shorter than the Hubble time), the SFR is set by self-regulation via feedback,  \n2  \nFigure 37. Effects of different stellar feedback physics on a re-start of an identical MW-mass simulation ( m12i ; as Fig. 35); we disable different mechanisms in turn as in Fig. 36. (1) Default: All physics included. (2) No OB/AGB mass-loss: Since the simulations start from identical initial conditions (identical gas supply), removing stellar mass-loss initially leads to a higher SFR (the OB winds cannot act as a source of feedback) but by ∼ 500Myr later, the net effect of removing the winds is to slightly decrease the SFR, as the winds cannot act continuously over this time to supply new gas mass to the disk. (3) No (Local) Radiative feedback: Here we still include a UV background, but remove only local (star particle) radiative feedback. We see a steady rise in SFR until this is forming more stars than the no-SNe case! Two effects occur: (a) GMCs are no longer efficiently disrupted by radiative feedback before SNe explode, therefore the SFR needed to self-regulate within a disk of a given surface density rises (see e.g. Hopkins et al. 2011), and (b) the outer atomic disk, which was stabilized by photo-heating from young stars in the galaxy, can now cool efficiently and form stars. (4) No Supernovae: The SFR is systematically larger, as expected, but by a relatively modest factor ∼ 2. This is because other feedback mechanisms instantaneously can regulate the SFR of gas in the galactic disk. However, removing SNe eliminates almost all galactic winds - on cosmological timescales, this means the disk will increase in gas mass and the SFR will run away, as we saw in Fig. 36. Critically, this demonstrates that how feedback regulates SF instantaneously within the disk (e.g. the Kennicutt-Schmidt relation) is not the same as how feedback regulates SF integrated over cosmological growth timescales (e.g. the stellar mass-halo mass relation).  \n<!-- image -->  \nat the point where the momentum flux per unit area from feedback, ∼ ( p ∗ / m ∗ ) ˙ Σ ∗ (where p ∗ / m ∗ ∼ 3000kms -1 is the time and IMF-averaged flux of momentum into the dense, star forming gas per unit stellar mass formed) offsets the gravitational force per unit area ∼ G Σ disk Σ gas (an identical scaling is obtained assuming feedback must destroy star-forming clouds, and/or offset turbulent dissipation, in a Q ∼ 1 disk; see Shetty & Ostriker 2008; Kim et al. 2011; Faucher-Giguère et al. 2013; Hayward & Hopkins 2017). Thus 〈 ˙ M ∗ 〉 ∼ ( p ∗ / m ∗ ) -1 G Σ disk 〈 M gas 〉 . Of course, these two scalings must be consistent: if the 'instantaneous' ˙ M ∗ is 'too low' compared to the cosmological scaling above, then ˙ M in is not offset by star formation+outflow and gas 'piles up,' raising M gas until the SFR'catches up.' If the radiative feedback contributes significantly to p ∗ / m ∗ (either directly or via pre-processing GMCs before SNe), then removing it lowers p ∗ / m ∗ and in turn increases the SFR given fixed Σ disk and M gas, as in Fig. 37. However, if η is less sensitive to radiative feedback, the winds from higher ˙ M ∗ would be stronger, and so the galaxy will self-regulate at the same 〈 ˙ M ∗ 〉 . Thus, the instantaneous self-regulation (the Kennicutt-Schmidt relation) is not  \nthe same as self-regulation of cosmological growth (e.g. the stellar mass-halo mass relation).", '10.1 Overview': "We present the FIRE-2 simulations (Table 1), a suite of cosmological simulations of galaxy formation using the FIRE physics modules in the GIZMO code. The FIRE-2 suite represents an update over FIRE-1, primarily in the use of a newer, more accurate hydrodynamic method, together with other numerical improvements to the physics algorithms and resolution - but no significant change in the actual physics simulated (§ 2). This includes high-resolution simulations run to z = 0 with full baryonic physics of the cold, multi-phase ISM, at both the dwarf and MW mass scales (reaching sub-pc spatial and ∼ 100 M glyph[circledot] mass scales; see Table 3 for details).  \nIn all properties investigated here, our primary conclusions from FIRE-1 appear qualitatively robust to these improvements in numerical accuracy (see § 3, Figs. 5-7), with typically less than factor of a few changes in global galaxy properties such as stellar mass, metallicity, and star formation rates. With explicit treatment of stellar feedback from SNe Types Ia & II, OB & AGB stellar winds, radiation pressure (UV, optical, and IR), photo-ionization and photo-electric heating, and explicit resolution of the multiphase ISM, galactic winds emerge naturally, producing galaxies with morphologies (Figs. 1 & 3), internal ISM structure (Fig. 2), flat or rising star formation histories and flat rotation curves (Figs. 4-5), metallicities (Fig. 4), and stellar masses (Fig. 7) apparently consistent with observations. Of course, it is impossible to consider an exhaustive list of galaxy properties here, and many more remain to be investigated in future work. Still, this is remarkable, considering that there is no fine-tuning or direct calibration of any parameters in the simulations to match these observations. Each of these properties has been investigated in more detail in previous work, described in § 1.  \nIt is worth noting that, as previous papers have pointed out (Hopkins et al. 2014; Ma et al. 2017b; Wetzel et al. 2016), there is no tension between even very thin stellar disks at z = 0 (e.g. Fig. 1), and strong stellar feedback which drives bursty star formation and strong galactic outflows at high redshift, strongly suppressing the stellar masses at z = 0 relative to no-feedback models. This owes to a combination of resolved venting of hot winds through a multi-phase ISM, which remove mass without disturbing galaxy morphologies, and a rapid transition from the bursty starforming mode and more 'gentle' steady-state thin-disk mode as gas fractions and specific star formation rates decline in massive, low-redshift halos (see Muratov et al. 2015; Hayward & Hopkins 2017; Faucher-Giguère 2018; Ma et al. 2017b).  \nWe have provided a complete algorithmic description of all aspects of the FIRE simulations (see § 2 & Appendices A-F). We also make our hydrodynamics+gravity code GIZMO, and all initial conditions for the runs here, public.  \nWe have considered an extensive study varying numerical and physical aspects of our simulations, to identify the most important 'ingredients' and physical effects in the simulations. Table 2 summarizes our conclusions; below we discuss the most important.", '10.2 Resolution': "(i) Mass resolution is most important for the physics that can be captured in Lagrangian methods; in § 4.1 we present a number of new resolution criteria relevant for simulations of resolved  \ngalaxy formation, phase structure, and galactic winds. For example, for dwarf galaxies, stellar masses are robust to resolution once the galaxies contain just a few star particles (Fig. 8). Metallicity and the shape of the SFH converge more slowly, requiring ∼ 100 star particles to sample the self-enrichment history and a baryonic particle mass glyph[lessorsimilar] 10 -6 M halo (Eq. 7) to avoid numerically enhanced burstiness (from single star particles representing 'too many' SNe at once). Accurately capturing phase structure in the ISM and gravitational fragmentation requires resolving the Toomre mass in the ISM, i.e. the largest GMCs, which dominate the galactic SFR (Eq. 5; Fig. 13). Morphology is very robust to resolution for dwarfs; for massive galaxies, it appears to depend on galaxy formation history - some massive systems are disky (with similar sizes and rotation curves) at all resolution levels, others (where the disk forms somewhat later) are more sensitive (Fig. 9). Details of galactic winds, in particular hot gas 'venting' and recycling in the CGM, are most sensitive to resolution of the properties we consider, with true 'convergence' likely requiring at least the ability to resolve the cooling radii of individual SNe (Eq. 6).  \n(ii) Force softenings for collisionless (DM+stellar) particles have little effect on our predictions (including galaxy baryonic properties, DM halo mass profiles, sub-structure mass and velocity distribution functions), provided they are chosen within a broad range (§ 4.2, Figs. 15-19). For DM this is between ∼ ( 0 . 005 -0 . 5 ) r 0 . 06, where r 0 . 06 is approximately the radius containing ∼ 200 DM particles at the halo center at z = 0 (§ 4.2.2-4.2.5). Our 'default' DM and stellar softening choices maximize the integration accuracy and convergence of our runs; larger softenings can artificially suppress central densities of DM halos.  \n(iii) We show that the radius r 0 . 06 is the radius to which DM profiles are converged to ∼ 10% or better - this is smaller than the commonly quoted 'convergence radius' of Power et al. (2003), owing to more accurate integration and timestepping, and smoother kernels (Figs. 11-12). At these spatial and mass resolution scales, N -body heating is negligible (Eq. 10-11). Owing to this insensitivity, it makes little difference whether we adopt constant softenings or fully-adaptive softenings (set to a fixed multiple of the local inter-particle separation) for collisionless particles.  \n(iv) Force softenings for gas should be set adaptively, so that the same mass distribution is being treated in the hydrodynamic equations and gravity equations (§ 4.2). With such a choice, turbulent fragmentation can be resolved down to small spatial scales, and no additional resolution criteria need be applied; of course, the meaning of 'spatial resolution' is adaptive, so Table 3 summarizes the effective spatial resolution of our simulations in different regimes. In dense, star-forming gas, the effective resolution reaches ∼ 1 -10 pc in MW-mass systems and ∼ 0 . 1 -1 pc in dwarfs. Because the medium is super-sonically turbulent, accurately capturing turbulent fragmentation and galactic star formation does not require that the thermal Jeans mass/length be resolved.  \n(v) If, for some reason, constant gravitational force softening for gas is desired in simulations with a resolved multiphase ISM, we show (Figs. 15-16) that it should be chosen sufficiently small that (1) the Toomre length scale (molecular gas disk scale height) is resolved ( glyph[epsilon1] glyph[lessmuch] 100pc), and (2) the minimum 'threshold' density for star formation is resolved gravitationally , so that actual self-gravitating regions are what form stars ( glyph[epsilon1] glyph[lessorsimilar] 7pc m 1 / 3 i , 1000 ( n SF , min / 100cm -3 ) -1 / 3 ). If the latter criterion is violated then star formation is driven (incorrectly) entirely by global (galaxy-wide) contraction or shock-compression of gas, and has nothing to do with local fragmentation and collapse of GMCs.  \n(vi) In addition to standard timestep criteria, guaranteeing  \nproper time-resolution of feedback requires a limiter (Eq. 12) such that star particles cannot 'skip' stages of stellar evolution or experience glyph[greatermuch] 1 SNe per particle per timestep (§ 4.3).  \n(vii) Moreover, if one adopts adaptive gravitational softening for collisionless particles, we show that standard timesteps based on the acceleration are numerically unstable. The reason is that the accelerations and potential now depend explicitly on the local particle configuration which changes as particles move through one another, not just on the total potential gradient; this means a Courantlike condition must be used to ensure stability (Fig. 20; Eq. 13).", '10.3 Hydrodynamic Methods': "(i) Comparing our new mesh-free finite-volume Godunov hydrodynamic methods ('MFM') to SPH (§ 5), we show that at all resolution levels, the properties of dwarf and high-redshift galaxies are insensitive to the details of the hydrodynamic solver (Fig. 21). However, for MW-mass and larger galaxies at late times, when a 'hot halo' develops, the details of the hydrodynamics become important for the late-time steady-state SFR of the galaxies, and their gas disk sizes (and subsequent circular velocity profiles). Similar conclusions have been reached in previous comparisons of movingmesh and fixed-grid codes to SPH, but with simpler treatments of ISM cooling/phases and feedback (Torrey et al. 2012; Sijacki et al. 2012; Kereš et al. 2012; Vogelsberger et al. 2012). We show that the gross morphology and instantaneous SFR given some ISM gas supply are not particularly sensitive (Figs. 23-24). Instead, the difference owes to the mixing and recycling of galactic winds in the CGM and initial accretion and shocking of hot gas, and its subsequent non-linear effect on cooling and re-accretion (Figs. 25-26). SPH - even in our state-of-the-art implementation - can produce some spurious non-mixing 'blobs' in the CGM, smear out accretion shocks, and suppress sub-sonic turbulence, leading to more cold and less hot, dense gas in the CGM. We explicitly show that by adding stronger 'artificial conductivity' to our SPH runs, we can essentially reproduce our MFM results (Fig. 22), indicating that details of numerical diffusion are important in the CGM of 'hot halos.'  \n(ii) We note that it is not appropriate to use 'artificial pressure floors' when a sink-particle model (like ours) for star formation (which identifies self-gravitating regions to turn into stars) is already present (§ 6). However, for modest values often used in the literature, we show such a floor has weak effects - but making the pressure floors too large produces clear numerical artifacts and prevents proper physical convergence in treating turbulent fragmentation (Fig. 27).  \n(iii) In finite-mass numerical methods without explicit advective fluxes (MFM), passive scalars (e.g. metals) remain locked (by default) to gas elements where they are first injected; in methods with such fluxes (e.g. AMR) there is a significant numerical mixing between elements. At finite resolution, these tend to under and overestimate (respectively) true physical mixing by un-resolved turbulent eddies and microphysical diffusion between neighbor cells. For AMR codes reducing the error requires higher resolution; in MFM one can add an explicit metal diffusivity. However, while this can be important for predictions of detailed abundance pattern distributions within galaxies, we show that this has essentially no effect on any properties studied here (§ 7.2; Fig 30-31), regardless of the numerical implementation (Fig. 32).  \n(iv) Given the sensitivity of the CGM in hot halos to fluid mixing details, it is likely that additional physics - magnetic fields, anisotropic Spitzer-Braginskii conduction and viscosity, cosmic ray  \ntransport - may produce larger changes in the late-time cooling rates from hot halos than the difference between hydro solvers.", '10.4 Cooling, Star Formation, & Feedback': "- (i) Consistent with previous work (Piontek & Steinmetz 2009; Hopkins et al. 2011, 2012c; Glover & Clark 2012; Hu et al. 2016), we find that most details of radiative cooling and gas chemistry generally have little effect on galaxy dynamics, star formation, and galactic winds (§ 7; Fig. 27-28). This is because almost all gas within the ISM has cooling times much shorter than dynamical times, so the exact cooling time and/or temperature is dynamically irrelevant. This is especially true in the 'cold' and 'cool' phases of the ISM/CGM/IGM ( T glyph[lessorsimilar] 10 5 K). Detailed variations in yields, which species are tracked, and numerical metal-mixing produce correspondingly weak effects.\n- (ii) However, in hot gas ( T glyph[greatermuch] 10 6 K) in the CGM or SNeheated bubbles, the cooling time can be longer than the dynamical time, and so high-temperature metal-line cooling has a significant effect on the phase structure of the CGM (Fig. 29) and cooling onto (hence SF in) the galaxy from the halo (Fig. 28; see also Choi & Nagamine 2009; Schaye et al. 2010). To leading order the total metallicity is the important quantity here, while detailed abundance-ratio variations have a much smaller effect.\n- (iii) As we have shown in previous work, provided the largest (Toomre) scales of fragmentation are resolved ( gravitationally and hydrodynamically ), and star formation occurs in that fragmenting gas, the resolution-scale star formation prescription has essentially no detectable effects on our predictions for galaxy-integrated SFRs, stellar masses, sizes, positions on the Schmidt-Kennicutt law (both galaxy-wide and spatially-resolved), galactic winds, and more (see § 8, Figs. 33-34, and Saitoh et al. 2008; Shetty & Ostriker 2008; Hopkins et al. 2011, 2013d; Kim et al. 2011; Agertz et al. 2013). This includes orders-of-magnitude variations in the density threshold and resolution-scale 'rate per free-fall time.' The key criteria are that these densities can (of course) be resolved, and that the SF occur only in self-gravitating gas above the mean galaxy density, so that the SF is automatically clustered (since the natural clustering of star formation non-linearly influences the likelihood of e.g. SNe bubbles overlapping and driving superwinds). Of course, predictions within the dense gas (corresponding to the densities where our sub-grid SF models take over) will be sensitive to the exact choices made; our default model is chosen based on previous work arguing it best reproduced observations of dense gas tracers (HCN and CO(3-2)) in idealized (non-cosmological) simulations (Hopkins et al. 2013e).\n- (iv) Provided the above criteria are met, SF is instead feedbackregulated . Uniformly increasing/decreasing the strength of feedback directly shifts both the integrated stellar mass, and position of galaxies on the KS law (Fig. 35).\n- (v) Stellar mass-loss (OB/AGB winds) is primarily important as a late-time fuel source for star formation (§ 9; Figs. 36-37); its net effect is to increase , not decrease, late-time SFRs, especially in massive galaxies.\n- (vi) Radiative feedback has a strong effect regulating the instantaneous SFRs of galaxies (i.e. the Kennicutt-Schmidt law), though less so on cosmologically-averaged SFRs (which are also regulated by the availability of fresh gas from IGM accretion).\n- (vii) The most important form of radiative feedback for dwarfs is external radiative feedback, i.e. the UV background. In Paper III, we show that removing the UVB has a much larger effect on the stellar masses of dwarfs than removing internal radiative feedback,  \neven at mass scales as large as M ∗ glyph[greaterorsimilar] 10 7 M glyph[circledot] . We will explore these effects further in future work.  \n(viii) On cosmological scales, SNe (specifically core-collapse and prompt Ia) are the most important feedback mechanism regulating galaxy properties. Without SNe (even given other feedback mechanisms), galaxies (especially dwarfs) form too many stars, and form these stars far too early. In Paper II, we therefore investigate numerical SNe treatments in more detail.", '10.5 Summary of Ingredients': 'To summarize, we find that the following criteria are essential for physically realistic high-resolution (multi-phase ISM) zoom-in simulations of galaxy formation:  \n- (i) The Toomre mass/length is mass and force-resolved; disk scale heights and at least the largest scales of fragmentation are resolved, with gravitational softening able to follow the mass.\n- (ii) In addition to standard (optically thin, primordial) cooling, some accounting for self-shielding and low-temperature ( T glyph[lessmuch] 10 4 K) cooling is included (to enable fragmentation), and hightemperature metal-line cooling is included to account for faster cooling in metal-enriched halo gas.\n- (iii) Star formation is restricted to gas at (force-resolved) densities significantly larger than the mean (ideally, to self-gravitating gas with a sink-particle type approach), so that it will occur in said fragments and therefore naturally be clustered, as observed.\n- (iv) Stellar mass loss is included, to provide a continuous additional fuel supply to the galaxy.\n- (v) Radiative feedback is included, particularly heating from the meta-galactic UV background, as well as photo-ionization heating and single-scattering photon momentum (radiation pressure) from stars in the simulation.\n- (vi) SNe (Ia & II) are included, with a coupling algorithm that carefully ensures manifest conservation of energy, mass, and momentum. More importantly, a careful accounting of both the energy and momentum budget of the coupled terms that properly treats whether the coupled radii are inside or outside the cooling radius, is necessary to obtain converged solutions (see Paper II).\n- (vii) All feedback quantities follow standard stellar evolution models for a standard IMF.', '10.6 Future Work': "In companion papers, we will study in more detail how mechanical feedback (SNe & stellar mass-loss; Paper II) and radiative feedback (photo-heating & photon momentum; Paper III) influence our predictions - both the details of their physics and numerical implementations. Because we argue that feedback is more important than many numerical details, it is extremely important to treat it as accurately as possible.  \nAs discussed above, more work is clearly warranted to investigate how additional fluid microphysics (e.g. magnetic fields) alters fluid mixing and subsequent cooling from the CGM. A major motivation of our switch to the new hydrodynamics solver in FIRE-2 is that it allows us to incorporate such physics in future work. A preliminary investigation of some of these physics is presented in Su et al. (2017). However this was primarily focused on the ISM inside galaxies; detailed, higher-resolution CGM studies are clearly necessary.  \n- In this paper, we study only systems with halo masses glyph[lessorsimilar] 2 × 10 12 M glyph[circledot] . This is because it is widely believed that feedback from supermassive black holes (not included in our simulations  \nhere) is critical to explain the properties (especially the quenching of star formation, and further quiescence) of more massive 'red and dead' systems (see Croton et al. 2006; Hopkins et al. 2006). Our preliminary studies from FIRE-1 support the idea that stellar feedback alone cannot explain all the observed properties of the most massive galaxies (Hopkins et al. 2014; Feldmann et al. 2017). In future work, we study the effects of black hole feedback on galaxy properties (for preliminary results see Anglés-Alcázar et al. 2017b).  \nFinally, we focus here almost exclusively on numerical studies. A series of papers will present the scientific predictions of these simulations for current and future observations.", 'ACKNOWLEDGMENTS': "We thank our editor and referee, Matthieu Schaller, for a number of helpful suggestions. Support for PFH and co-authors was provided by an Alfred P. Sloan Research Fellowship, NASA ATP Grant NNX14AH35G, and NSF Collaborative Research Grant #1715847 and CAREER grant #1455342. ARW was supported by a Caltech-Carnegie Fellowship, in part through the Moore Center for Theoretical Cosmology and Physics at Caltech. CAFG was supported by NSF through grants AST-1412836 and AST-1517491, and by NASA through grant NNX15AB22G. The Flatiron Institute is supported by the Simons Foundation. MBK and AF were partially supported by the NSF through grant AST-1517226. MBK also acknowledges support from NASA through HST theory grants (programs AR-12836, AR-13888, AR-13896, and AR-14282) awarded by STScI. DK was supported by NSF Grant AST1412153 and a Cottrell Scholar Award from the Research Corporation for Science Advancement. RF was supported by the Swiss National Science Foundation (grant No. 157591). Numerical calculations were run on the Caltech compute cluster 'Wheeler,' allocations TG-AST130039 & TG-AST150080 granted by the Extreme Science and Engineering Discovery Environment (XSEDE) and PRAC NSF.1713353 supported by the NSF, and the NASA HEC Program through the NAS Division at Ames Research Center and the NCCS at Goddard Space Flight Center.", 'REFERENCES': '- -. 2011a, MNRAS, 416, 1354', 'APPENDIX A: APPROXIMATE STELLAR EVOLUTION TABULATIONS': "Here we present simple fits to the stellar evolution models and yields used in the FIRE simulations. We note that in several cases the simulations in this paper utilize a more detailed look-up table; however for all practical purposes the fits here are sufficiently accurate that the differences are negligible. Stellar evolution results are obtained from STARBURST99 (Leitherer et al. 1999) assuming a Kroupa (2001) IMF. SNe Ia rates follow Mannucci et al. (2006) including both prompt and delayed populations. Yields for corecollapse SNe are IMF-averaged for the same IMF, from the tables (including hypernovae) in Nomoto et al. (2006). Yields for SNe Ia follow Iwamoto et al. (1999). Yields for OB/AGB winds are taken from the synthesis of the models from van den Hoek & Groenewegen (1997), Marigo (2001), and Izzard et al. (2004) as synthesized in Wiersma et al. (2009b), appropriately re-averaged over the IMF.  \n- (i) SNe Ia: These occur with rate-per-unit-stellar mass R Ia = dN Ia / dt = 0 for t Myr < 37 . 53 (where t Myr is the age of the star particle in Myr), then R Ia / ( SNeMyr -1 M -1 glyph[circledot] ) = 5 . 3 × 10 -8 + 1 . 6 × 10 -5 exp {-[ ( t Myr -50 ) / 10] 2 / 2 } for t Myr ≥ 37 . 53. Since the rate is per-stellar-mass, the expectation value of the number of SNe for a star particle of mass mb and timestep ∆ tb is = R Ia ( t b Myr ) mb ∆ tb , and the trial for 'success' (an explosion) is drawn from a binomial distribution. Each SN Ia has ejecta mass M ej = 1 . 4 M glyph[circledot] and energy E ej =( 1 / 2 ) M ej v 2 ej = 10 51 erg. The ejecta yield mass for the tracked species is: (Z, He, C, N, O, Ne, Mg, Si, S, Ca, Fe) = ( 1 . 4 , 0 , 0 . 049 , 1 . 2 × 10 -6 , 0 . 143 , 0 . 0045 , 0 . 0086 , 0 . 156 , 0 . 087 , 0 . 012 , 0 . 743 ) M glyph[circledot] .  \n(ii) SNe II: The core-collapse SNe rate can be surprisingly wellfit by a simple piecewise-constant function, R II = 0 for t Myr < 3 . 401 or t Myr > 37 . 53; R II / ( SNeMyr -1 M -1 glyph[circledot] ) = 5 . 408 × 10 -4 for 3 . 4 < t Myr < 10 . 37; and R II / ( SNeMyr -1 M -1 glyph[circledot] ) = 2 . 516 × 10 -4 for 10 . 37 < t Myr < 37 . 53. The IMF-averaged ejecta mass per explosion is M ej = 10 . 5 M glyph[circledot] , with ejecta energy E ej = ( 1 / 2 ) M ej v 2 ej = 10 51 erg. These are normalized so that the ejecta mass and energy exactly match the integrated totals from STARBURST99; the timeaveraged energy injection rate ( 〈 R II E ej 〉 ) is within ∼ 10% of the tabulated STARBURST99 rate at all times. The IMF-averaged yields  \nare (He, C, N, O, Ne, Mg, Si, S, Ca, Fe) = ( 3 . 87 , 0 . 133 , 0 . 0479 ˜ N , 1 . 17 , 0 . 30 , 0 . 0987 , 0 . 0933 , 0 . 0397 , 0 . 00458 , 0 . 0741 ) M glyph[circledot] , where ˜ N = MAX ( Z / Z glyph[circledot] , 1 . 65 ) accounts for the strongly progenitormetallicity dependent N yields. The total metal yield Z is given by summing the explicitly-followed species, with an additional ∼ 2% added to account for un-tracked species ( Z = 1 . 02 ∑ Zi , followed); the remaining ejecta is H. 31  \n(iii) OB/AGB Mass-Loss: We include all non-SNe mass-loss channels here, but this is dominated by OB/AGB-star winds. The IMF-integrated mass-loss rate for a stellar population/particle of mass M ∗ is ˙ M w = M ∗ f w Gyr -1 with f w = 4 . 763 ( 0 . 01 + Z / Z glyph[circledot] ) for t Myr < 1; f w = 4 . 763 ( 0 . 01 + Z / Z glyph[circledot] ) t 1 . 45 + 0 . 8 ln ( Z / Z glyph[circledot] ) Myr for 1 < t Myr < 3 . 5; f w = 29 . 4 ( t Myr / 3 . 5 ) -3 . 25 + 0 . 0042 for 3 . 5 < t Myr < 100; and f w = 0 . 42 ( t Myr / 1000 ) -1 . 1 / ( 19 . 81 -ln ( t Myr )) for t Myr > 100. The total (IMF-averaged) kinetic luminosity of the massloss is given by L kinetic = ( 1 / 2 ) ˙ M w 〈 v 2 w 〉 = ψ × 10 12 ergg -1 ˙ M w, with ψ =( 5 . 94 × 10 4 ) / ( 1 +( t Myr / 2 . 5 ) 1 . 4 +( t Myr / 10 ) 5 . 0 )+ 4 . 83 for t Myr < 100 and ψ = 4 . 83 for t Myr > 100. The yields are given by the maximum of either the progenitor stellar surface abundances or, for the light species (He, C, N, O), mass fractions = ( 0 . 36 , 0 . 016 , 0 . 0041 ˜ N , 0 . 0118 ) .  \n(iv) Radiation: Wedefine the light-to-mass ratio in a given band Ψ band, with units L glyph[circledot] / M glyph[circledot] . Then the bolometric Ψ bol = 1136 . 59 for t Myr < 3 . 5, and Ψ bol = 1500 exp[ -4 . 145 x + 0 . 691 x 2 -0 . 0576 x 3 ] with x ≡ log 10 ( t Myr / 3 . 5 ) for t Myr > 3 . 5. For the bands used in our radiation hydrodynamics, we have the following intrinsic (before attenuation) bolometric corrections. In the mid/far IR, Ψ IR = 0. In optical/NIR, Ψ opt = f opt Ψ bol with f opt = 0 . 09 for t Myr < 2 . 5; f Opt = 0 . 09 ( 1 + [ ( t Myr -2 . 5 ) / 4] 2 ) for 2 . 5 < t Myr < 6; f Opt = 1 -0 . 841 / ( 1 + [ ( t Myr -6 ) / 300] ) for t Myr > 6. For the photo-electric FUV band Ψ FUV = 271[1 +( t Myr / 3 . 4 ) 2 ] for t Myr < 3 . 4; Ψ FUV = 572 ( t Myr / 3 . 4 ) -1 . 5 for t Myr > 3 . 4. For the ionizing band Ψ ion = 500 for t Myr < 3 . 5; Ψ ion = 60 ( t Myr / 3 . 5 ) -3 . 6 + 470 ( t Myr / 3 . 5 ) 0 . 045 -1 . 82 ln t Myr for 3 . 5 < t Myr < 25; Ψ ion = 0 for t Myr > 25. The remaining UV luminosity, Ψ bol -(Ψ IR +Ψ opt + Ψ FUV +Ψ ion ) is assigned to the NUV band Ψ NUV. The fluxmean dust opacities adopted are ( κ FUV , κ NUV , κ opt , κ IR ) = ( 2000 , 1800 , 180 , 10 )( Z / Z glyph[circledot] ) cm 2 g -1 . The photo-ionization rate (and corresponding κ ion) is calculated from the neutral hydrogen density as described in Appendix E and Appendix B below.", 'APPENDIX B: APPROXIMATE COOLING FUNCTIONS': "Here we provide simple fitting-function approximations to the complete set of cooling functions used in our FIRE simulations. Note that for several of these, we use somewhat more accurate lookup tables in the simulations (as a function of temperature, density, and metallicity), but we provide functions accurate to ∼ 10% over the relevant dynamic range in the simulations ( ∼ 10 -10 9 K), so  \n31 For reference, the yields from Woosley & Weaver (1995), used in FIRE1, are (He, C, N, O, Ne, Mg, Si, S, Ca, Fe) = ( 4 . 03 , 0 . 117 , 0 . 0399 ˜ N , 1 . 06 , 0 . 169 , 0 . 0596 , 0 . 0924 , 0 . 0408 , 0 . 00492 , 0 . 0842 ) M glyph[circledot] . In both FIRE1 and FIRE-2, we choose to omit the progenitor stellar metallicity dependence of the predicted yields (using the yields for solar-metallicity progenitors instead) if the standard deviation of the metallicity dependence | dMZ ( species ) / dZ progenitor | between the models of Chieffi & Limongi (2004), Woosley & Weaver (1995), and Nomoto et al. (2006) is larger than the magnitude of the actual predicted trend | dMZ ( species ) / dZ progenitor | from Nomoto et al. 2006. The only tracked species which passes this test is N. We do set the progenitor metallicity as a 'floor' to the yields, following exactly the algorithm in Wiersma et al. (2009b).  \nthat interested readers can reproduce our full cooling physics treatment.  \nThe instantaneous cooling+heating rate per unit volume is given by the sum over all processes,  \nde a thermal dt = -n 2 H , a Λ a net = -n 2 H , a ∑ i Λ i (B1)  \nwhere e a thermal is the thermal energy density of a gas element a . Here we will use Λ to denote both cooling and heating rates, but with opposite signs (a positive sign here denotes cooling). Also, define ˜ nx ≡ nx / n H as the number of species x per hydrogen nucleus (e.g. ˜ ne , ˜ n H0, and ˜ n HeI denote the free electron, neutral hydrogen, and neutral helium numbers). Below, all units are cgs ([ T ] = K, [ n H] = cm -3 , [ Λ ] = ergs -1 cm 3 ). The processes we track include:  \n- (i) Collisional Excitation : from Katz et al. (1996) (incorporating earlier fits from Cen 1992):  \nΛ CE = ( β H0 ˜ n H0 + β HeI ˜ n HeI) ˜ ne (B2)  \nβ H0 = 7 . 50 × 10 -19 τ 5 exp ( -118348 T ) (B3)  \nβ HeI = 5 . 54 × 10 -17 τ 5 T -0 . 397 exp ( -473638 T ) (B4)  \nτ 5 ≡ [ 1 + ( T 10 5 ) 1 / 2 ] -1 (B5)  \n(ii) Collisional Ionization : also from Katz et al. (1996):  \nΛ CI = 10 -11 ˜ ne × (B6)  \n(2  \n.  \n18  \nγ  \nH0 ˜ n  \nH0  \n+  \n3  \n.  \n94  \nγ  \nHe0 ˜ n  \nHe0  \n+  \n8  \n.  \n72  \nγ  \nHeI ˜  \nn  \nHeI)  \nγ H0 = 5 . 85 × 10 -11 T 1 / 2 τ 5 exp ( -157809 . 1 T ) (B7)  \nγ He0 = 2 . 38 × 10 -11 T 1 / 2 τ 5 exp ( -285335 . 4 T ) (B8)  \nγ HeI = 5 . 68 × 10 -12 T 1 / 2 τ 5 exp ( -631515 . 0 T ) (B9)  \n(iii) Recombination : from Verner & Ferland (1996):  \nΛ  \nRec = 1 . 036 × 10 -16 T ˜ ne × (B10) ( α HI ˜ n HI + [ α HeI + 629922 . 78 T α di ] ˜ n HeI + α HeII ˜ n HeII ) α HI = 7 . 982 × 10 -11 ( 1 . 774 T 0 . 5 ) × (B11) ( 1 + T 0 . 5 1 . 774 ) -0 . 252 ( 1 + T 0 . 5 838 . 81 ) -1 . 748 α HeI = 9 . 356 × 10 -10 ( 0 . 2065 T 0 . 5 ) × (B12) ( 1 + T 0 . 5 0 . 2065 ) -0 . 2108 ( 1 + T 0 . 5 6063 . 0 ) -1 . 7892 α HeII = 1 . 5964 × 10 -10 ( 2 . 5092 T 0 . 5 ) × (B13) ( 1 + T 0 . 5 2 . 5092 ) -0 . 252 ( 1 + T 0 . 5 1677 . 6 ) -1 . 748 α di = 1 . 9 × 10 -3 T 1 . 5 exp ( -470000 T ) × (B14) ( 1 + 0 . 3 exp [ -94000 T ])  \nNote the α di term here comes from dielectric recombination.", '(iv) Free-free emission : from Rybicki & Lightman (1986):': 'Λ FF = β ff ( T ) ( ˜ n HI + ˜ n HeI + 4 ˜ n HeII) ˜ ne (B15)  \nβ ff ( T ) = 1 . 43 × 10 -27 T 1 / 2 × (B16)  \n[ 1 . 1 + 0 . 34 exp { -( 5 . 5 -log 10 [ T ] ) 2 / 3 }]  \n(v) High-Temperature Metal-Line Cooling : this refers to metal-line cooling processes in mostly ionized gas, with temperatures glyph[greaterorsimilar] 10 4 K. We use the public look-up tables from Wiersma et al. (2009a), for which:  \nΛ Metal = ∑ species i Λ i Metal = ∑ i ˜ ne ξ i ( n H , T , z ) Z i Z i glyph[circledot] (B17)  \nwhere this refers to the sum over all tracked metal species i (here C, N, O, Ne, Mg, Si, S, Ca, Fe), and Z i / Z i glyph[circledot] is the abundance of species i relative to solar. We adopt solar abundances (Z, He, C, N, O, Ne, Mg, Si, S, Ca, Fe) = ( 0 . 02 , 0 . 28, 3 . 26 × 10 -3 , 1 . 32 × 10 -3 , 8 . 65 × 10 -3 , 2 . 22 × 10 -3 , 9 . 31 × 10 -4 , 1 . 08 × 10 -3 , 6 . 44 × 10 -4 , 1 . 01 × 10 -4 , 1 . 73 × 10 -3 ), which give abundance ratios matching Asplund et al. (2009) scaled to the total metallicity = 0 . 02 (used because it matches the assumed values in the cooling computations and stellar evolution models). The functions ξ i depend on density, temperature, and redshift z (because they assumes photo-ionization by a redshift-dependent UV background); they are taken from the look-up tables provided by Wiersma et al. (2009a), at abundances Z = Z glyph[circledot] (defined above). 32  \n(vi) Low-Temperature Metal Line, Fine-Structure, & Molecular Cooling : this combines the gas-phase low-temperature cooling (including molecular and atomic processes) in mostly neutral gas below glyph[lessorsimilar] 10 4 K. From our compilation of CLOUDY runs (Ferland et al. 1998), fitting the resulting look-up tables, we obtain approximately:  \n-4 . 9202  \nΛ Cold = 2 . 896 × 10 -26 { ( T 125 . 215 ) + (B18) ( T 1349 . 86 ) -1 . 7288 + ( T 6450 . 06 ) -0 . 3075 } -1 × ( 1 +( Z / Z glyph[circledot] ) 1 + 0 . 00143 n H ) (1 -f selfshield) × ( 0 . 001 + 0 . 10 n H 1 + n H + 0 . 09 n H 1 + 0 . 1 n H + ( Z / Z glyph[circledot] ) 2 1 + n H ) × exp ( -[ T 158000 ] 2 )  \nwhere f selfshield accounts for the local radiation environment by applying a simple (fitted) local shielding correction for UV/ionizing photons, f selfshield ≡ exp ( -˜ τ ion a ) with ˜ τ ion a ≡ σ H ν 0 n H , a glyph[lscript] fit a where σ H ν 0 ≡ 6 × 10 -18 cm -2 and glyph[lscript] fit a ≡ 4 . 4pc ( T / 10 4 K ) -0 . 173 Γ -2 / 3 -12 ( Γ -12 is the ionization rate in units of 10 -12 s -1 , including both the UV background and local sources assuming they have the same spectral shape, as defined below for photo-ionization heating).  \n(vii) Dust Collisional Heating/Cooling : from Meijerink & Spaans (2005):  \nΛ Dust = 1 . 12 × 10 -32 ( T -T dust) T 1 / 2 × (B19) ( 1 -0 . 8 exp [ -75 T ]) ( Z Z glyph[circledot] )  \nwhere we take T dust = 30K here, and the Z / Z glyph[circledot] term comes from assuming a constant dust-to-metals ratio.  \n(viii) Compton Heating/Cooling : from the CMB, gives (Rybicki & Lightman 1986):  \nΛ Compton = 5 . 65 × 10 -36 ˜ ne ( T -T CMB[ z ] )( 1 + z ) 4 n -1 H (B20)  \n(ix) Photo-Ionization Heating : from the UVB and local (insimulation stellar sources) gives a heating rate, hence negative Λ , of  \nΛ ion = -˜ f ( glyph[epsilon1] H0 ˜ n H0 + glyph[epsilon1] He0 ˜ n He0 + glyph[epsilon1] HeI ˜ n HeI) n -1 H (B21)  \n˜ f ≡ ( 1 + e local ν , ion e UVB ν , ion ) f selfshield (B22)  \nlog 10 ( glyph[epsilon1] H0 ) ≈-24 . 6 + 1 . 62 x + 14 . 9 x 2 -45 . 5 x 3 (B23)  \n+ 46 . 2 x 4 -16 . 7 x 5 -exp[50 ( x -1 . 05 ) ]  \nlog 10 ( glyph[epsilon1] He0 ) ≈ log 10 ( glyph[epsilon1] H0 ) -0 . 0366 + 0 . 376 x  \nlog 10 ( glyph[epsilon1] HeI ) ≈-26 . 3 + 0 . 816 x + 78 . 2 x 2 -837 x 3 + 4770 x 4  \n-15600 x 5 + 29600 x 6 -32400 x 7 + 18900 x 8 -4550 x 9 x ≡ log 10 ( 1 + z ) (B24)  \nwhere glyph[epsilon1] H0, glyph[epsilon1] He0, and glyph[epsilon1] HeI are pre-tabulated for the assumed UV background magnitude and shape, as a function of redshift in Faucher-Giguère et al. (2009) 33 - values above are simple polynomial fits good to ∼ 10% up to z ∼ 10. The factor ˜ f accounts for both self-shielding (reducing the effective incident flux by f selfshield) and the contribution from local sources, where e local ν , ion is the ionizing band radiation energy density calculated explicitly from the radiation-hydrodynamic treatment in the code (Appendix E), and e UVB ν , ion is the meta-galactic UV background (UVB) energy density integrated in the same band for the same wavelength range (Hionizing frequencies). Note this means we assume the spectral slope of escaping, ionizing radiation from resolved stars in the simulation is the same as the UVB.', '(x) Cosmic Ray Heating : from Guo & Oh (2008):': "Λ CR = -1 . 0 × 10 -16 ( 0 . 98 + 1 . 65 ˜ ne XH ) e CR n -1 H (B25)  \nwhere we assume an approximately uniform MW-like cosmic ray background, e CR ≈ 9 . 0 × 10 -12 f CR. Here f CR = n H / ( 0 . 01 + n H ) when n H exceeds 1000 × the mean baryonic density of the Universe, and f CR = 0 otherwise, to avoid an artificially high CR heating rate in extremely low-density regions (e.g. outside galaxies) or at very high redshifts (before star formation).  \n(xi) Photo-Electric Heating : from Wolfire et al. (2003):  \nΛ PE = -1 . 3 × 10 -24 ˜ e pe ν n -1 H ( Z Z glyph[circledot] ) × (B26) ( 0 . 049 1 +( x pe / 1925 ) 0 . 73 + 0 . 037 ( T / 10 4 ) 0 . 7 1 +( x pe / 5000 ) ) x pe ≡ ˜ e pe ν T 0 . 5 0 . 5 ˜ ne n H (B27)  \nwhere ˜ e pe ν is the photon energy density in the photo-electric band, normalized to the Habing (1968) MW units, ˜ e pe ν ≡ e pe ν / ( 3 . 9 × 10 -14 ergcm -3 ) . The Z / Z glyph[circledot] term comes from assuming a constant dust-to-metals ratio. Here the field e pe ν is the FUV band radiation energy density calculated explicitly from the radiation-hydrodynamic treatment in the code, described in Appendix E.  \n(xii) Magneto-Hydrodynamic Work & Shocks : from the MHDequations, we obtain some fluid-dynamic change to the temperature (owing to compression/expansion, shocks, etc). We include this self-consistently in the fully-implicit temperature update:  \nΛ MHD = -µ ∂ u thermal ∂ t ∣ ∣ ∣ MHD n -1 H (B28)  \nwhere u thermal is the specific internal energy (internal energy per unit mass).  \n(xiii) Optically-Thick Cooling : lacking a full radiative transfer solution for cooling radiation, we approximate the effects of optically-thick cooling using the method from Rafikov (2007), which captures the most important effects by approximating each element as a 'slab' with column density estimated via the Sobolev approximation and integrating a vertical atmosphere through to its photosphere to determine the net photon escape. This amounts to first summing the contributions above to determine the net heating/cooling rate Λ Net, and then restricting this to the cooling rate of said slab:  \n| Λ Net | < Λ BB (B29)  \nΛ BB ≡ 5 . 67 × 10 -5 T 4 ( µ Σ eff ) 1 1 + κ eff Σ eff n -1 H (B30)  \nwhere Σ eff = 〈 Σ a , Sobolev column 〉 φ , θ = ρ a ( ha + ρ a / |∇ ρ a | ) uses the local Sobolev approximation to estimate the column density to infinity and is defined in Appendix E, and κ eff is the effective opacity. 34  \nAs noted in the text, the actual heating/cooling step is solved fully implicitly for each gas element on its own timestep.", 'APPENDIX C: ALGORITHMIC IMPLEMENTATION OF STAR FORMATION': '(i) Self-Gravitating: First, following standard sink-particle approaches, we calculate the virial parameter α (ratio of thermal plus kinetic energy to potential energy) in a resolution element and allow only star formation in bound particles with α < 1. From Hopkins et al. (2013d),  \nα ≡ ‖∇⊗ v ‖ 2 a +( cs , a / ha ) 2 8 π G ρ a (C1)  \nwhere ‖∇⊗ v ‖ a the Frobenius norm ( ‖ A ‖ 2 ≡ ∑ αβγ ... A 2 αβγ ... ) of the velocity gradient tensor ( ⊗ is the outer product), cs , a is the sound speed, ρ a the density, and ha the usual resolution scale (interparticle spacing). This has the advantage that it converges to an explicitly resolution-independent expression in the super-sonic turbulence limit (see Hopkins et al. 2013d). Note that in GIZMO, we always use the higher-order accurate, matrix-based gradient estimators described in Hopkins (2015), which remain second-order accurate, consistent, and robust despite arbitrary particle configurations within the stencil (see also Maron & Howes 2003; Luo et al. 2008; Lanson & Vila 2008; Mocz et al. 2014; Pakmor et al. 2016).  \n34 We take κ eff for dust ( T < 1500) from the detailed tables in Semenov et al. (2003), assuming the dust, gas, and radiative temperatures are in equilibrium (true in the optically thick limit at these temperatures), which is approximately well-fit by κ eff = 5 for 150 ≤ T ≤ 1500 and κ eff = 0 . 0027 T 3 / 2 for T < 150. At higher temperatures the system is rarely optically thick, but for completeness we compute κ eff from the gas-phase using standard approximations for stellar atmospheres: κ -1 eff = κ -1 rad + κ -1 cond , with κ cond = 2 . 6 × 10 -7 ˜ ne T 2 ρ -2 , κ rad = κ mol + 1 / ( κ -1 H -+ [ κ e + κ Kr] -1 ) , κ mol = 0 . 1 Z , κ e = 0 . 2 ( 1 + XH ) , κ H -= 1 . 1 × 10 -25 ( Z ρ ) 1 / 2 T 7 . 7 , κ Kr = 4 . 0 × 10 25 ( 1 + XH ) Z ρ 3 . 5 .  \nUsing the zeroth-order inaccurate SPH gradient estimator, in contrast, gives similar results statistically, but makes identification of individual physically collapsing clouds much more noisy.  \n(ii) Self-Shielding: If α < 1, we next calculate the shielded/molecular fraction f ( sf ) shielded , which is the fraction of the mass that should be able to self-shield and so cool efficiently (hence fragment to stellar mass scales). The expression for the shielded fraction from Krumholz & Gnedin (2011) is:  \nf ( sf ) shielded ≡ 1 -3 1 + 4 ˜ ψ a (C2)  \n˜ ψ a ≡ 0 . 6 ˜ τ a ( 0 . 01 + Za / Z glyph[circledot] ) ln ( 1 + 0 . 6 ˜ φ a + 0 . 01 ˜ φ 2 a ) (C3)  \n˜ φ a ≡ 0 . 756 ( 1 + 3 . 1 Za / Z glyph[circledot] ) 0 . 365 (C4)  \n˜ τ a ≡ 434 . 8cm 2 g -1 ρ a ( ha + ρ a |∇ ρ | a ) (C5)  \nWe require f ( sf ) shielded > 0 for star formation.  \n(iii) Jeans-Unstable: If α < 1 and f ( sf ) shielded > 0, we calculate the Jeans mass mJ , and only allow star formation in Jeans-unstable particles, specifically those where mJ < mJ , crit ≡ MAX ( ma , 10 3 M glyph[circledot] ) , where ma is the particle mass. We calculate the Jeans mass as  \nmJ = 2 M glyph[circledot] ( cs , a 0 . 2kms -1 ) 3 ( na 10 3 cm -3 ) -1 / 2 (C6)  \nwhere na ≡ ρ a /µ a is the gas number density.  \n(iv) Sufficiently-Dense: If α < 1 , f ( sf ) shielded > 0, and mJ < mJ , crit, we check if na > n crit, where n crit = 1000cm -3 is a minimum density (and na ≡ na , H = XH ρ a / mp ), to prevent spurious triggering of the above criteria in low-density gas.  \nWe then assign the gas particle a volume-integrated SFR:  \n˙ m a ∗ =Θ( α a , f shielded , a , na , mJ , a ) f ( sf ) shielded m gas , a t freefall , a (C7)  \nt freefall , a ≡ √ 3 π 32 G ρ a (C8) Θ= { 1 ( α < 1 , f ( sf ) shielded > 0 , mJ < mJ , crit , na > n crit ) 0 otherwise (C9)  \nBecause we wish to maintain equal stellar and gas element masses, at each timestep ∆ ta weassign the gas particle a probability pa = 1 -exp ( -˙ m a ∗ ∆ ta / m a gas ) of turning into a star particle that timestep; we draw a uniform random variable 0 < x < 1 and if x < pa , we convert the gas particle to a star particle. It inherits all relevant properties of its parent particle.', 'APPENDIX D: ALGORITHMIC IMPLEMENTATION OF MECHANICAL FEEDBACK': "Here we describe our implementation of mechanical feedback, used for SNe (Types Ia & II) and stellar mass-loss. This algorithm was first developed and presented in a series of papers, beginning with Hopkins et al. (2012c), and the version used in FIRE-1 (which contains most of the important features here up to some specific numerical improvements for FIRE-2) was presented in detail in Hopkins et al. (2014). Similar aspects of that algorithm - in particular the treatment of SNe momentum accounting for the terminal momentum - have been recently developed for other codes by Kimm & Cen (2014), Martizzi et al. (2015), and Rosdahl et al. (2017).  \nIn Paper II, we discuss each piece of this algorithm in detail,  \nand consider a large suite of idealized test problems and cosmological simulations, to test and validate each and show how it influences our predictions. However for the sake of completeness, we include the full algorithm here.  \n(i) Every timestep ∆ ta , for each star particle a (at position x a ), we first determine whether an 'event' occurs: a SN Ia, SN II, and/or non-zero stellar mass-loss. This follows the rates and algorithms in Appendix A. If an event occurs, it has some initial ejecta (or wind) mass m ej, metal mass mZ , ej (defined for each species we track), momentum p ej = m ej v ej, and energy E ej. These also are given in Appendix A (for winds, m ej =∆ t ˙ Mw from the star, for SNe it is the ejecta mass).  \n(ii) Identify gas elements surrounding the star particle: in a grid code this is straightforward, but in our mesh-free method, we define an effective neighbor number N ngb =( 4 π/ 3 ) H 3 a ¯ na ( Ha ) in the same manner as for hydrodynamics, where W is the kernel function, ¯ na = ∑ W ( x ba ≡ x b -x a , Ha ) , and Ha is the search radius. The equation for N ngb ( Ha ) is non-linear so is solved iteratively in the neighbor search; see Springel (2005). Thus we obtain all gas elements b within a radius | x ba | < Ha (where ' a sees b '); we also identify all neighbors with | x ba | < Hb (i.e. ' b sees a '). We show in Paper II that this is important to ensure the 'effective faces' close and the resulting distribution of ejecta is isotropic, in regions with highly disordered gas element positions.  \n(iii) Assign 'vector weights' to each neighbor, by first boosting to the rest-frame of the star ( x a = 0 , v a ≡ d x a / dt = 0 ; in which the ejecta should be isotropic), then calculating the 'effective face' that would be seen by the star particle (using the same definitions of inter-cell faces used in the hydrodynamics) and integrating the ejecta over solid angle through to each face. This amounts to defining the vector weight function ¯ w ba  \n¯ w ba ≡ w ba ∑ c | w ca | (D1)  \nw ba ≡ ω ba ∑ + , -∑ α (ˆ x ± ba ) α ( f α ± ) a (D2)  \n( f α ± ) a ≡ { 1 2 [ 1 + (∑ c ω ca | ˆ x ∓ ca | α ∑ c ω ca | ˆ x ± ca | α ) 2 ]} 1 / 2 (D3)  \nω ba = ∆Ω ba 4 π ≡ 1 2 ( 1 -1 √ 1 +( A ba · ˆ x ba ) / ( π | x ba | 2 ) ) (D4)  \nwhere A ba is the effective vector face between elements b and a used in the finite-volume hydrodynamic calculations, 35 and the ˆ x ± ca  \n35 For our MFM hydrodynamic method, the face A ba is defined as (see Hopkins 2015):  \nA ba ≡ ¯ n -1 a ¯ q b ( x a ) + ¯ n -1 b ¯ q a ( x b ) (D5)  \n¯ q b ( x a ) ≡ E -1 a · x ba W ( x ba , Ha ) (D6)  \nE a ≡ ∑ c ( x ca ⊗ x ca ) W ( x ca , Ha ) (D7)  \nFor SPH, the face is defined by the simpler relation A ba = [ ¯ n -2 a ∂ W ( | x | ba , Ha ) /∂ | x | ba +¯ n -2 b ∂ W ( | x | ba , Hb ) /∂ | x | ba ] ˆ x ba . In movingmesh or fixed-grid finite-volume codes, the face A ba is the explicit geometric mesh face between cells.  \nare the positive or negative (singly-signed) projection vectors:  \nˆ x ba ≡ x ba | x ba | = ∑ + , -ˆ x ± ba (D8)  \n(ˆ x + ba ) α ≡| x ba | -1 MAX ( x α ba , 0 ) ∣ ∣ ∣ α = x , y , z (D9)  \n(ˆ x -ba ) α ≡| x ba | -1 MIN ( x α ba , 0 ) ∣ ∣ ∣ α = x , y , z (D10)  \nThese expressions are complicated but are derived in detail in Paper II. Their important properties are (1) they maintain manifest conservation of mass, momentum, and energy (see below). (2) They give fluxes which are statistically isotropic in the rest frame of the star, i.e. the ejecta are not systematically biased in one direction or another, even if there is a global density gradient such that there are, on average, more gas elements in one direction. We demonstrate this in numerical tests explicitly in Paper II, and show that many simpler prescriptions lead to systematic, unphysical biases in the ejecta deposition, e.g. if there is a thin, dense disk such that more gas neighbors are 'in the disk,' simple weighting ¯ w ba proportional to, say, the SPH kernel, leads to almost all the ejecta being coupled in the disk, driving an expanding ring, with almost no ejecta going into the vertical direction - when in fact the converged solution to this problem is exactly the opposite (hot gas 'vents' in the vertical direction). (3) They approximate, as closely as possible without an expensive numerical quadrature, the exact integral of the ejecta through and into the 'domains' of each gas neighbor determined by the hydrodynamic volume partition.  \n(iv) Assign initial fluxes in the rest-frame of the star:  \n∆ mb = | ¯ w b | m ej (D11)  \n∆ mZ , b = | ¯ w b | mZ , ej (D12)  \n∆ Eb = | ¯ w b | E ej (D13)  \n∆ p b = ¯ w b p ej (D14)  \nIt is easy to see that our definition of ¯ w ba guarantees exact conservation of mass, energy, and linear momentum, and that the correct total radial (outward) momentum is assigned, e.g.:  \n∑ ∆ mb = m ej (D15)  \n∑ ∆ mZ , b = mZ , ej (D16)  \n∑ ∆ Eb = E ej (D17)  \n∑ | ∆ p b | = p ej (D18)  \n∑ ∆ p b = 0 (D19)  \n(v) Boost back to the simulation ('lab') frame: if the star is moving with velocity v a , then this boost transforms the momentum and energy fluxes:  \n∆ p ' b ≡ ∆ p b +∆ mb v a (D20)  \n∆ E ' b ≡ ∆ Eb + 1 2 ∆ mb ( | ∆ p ' b | 2 -| ∆ p b | 2 ) (D21)  \n(the mass fluxes are unchanged, ∆ m ' b =∆ mb , ∆ m ' Z , b =∆ mZ , b ). Of course this maintains manifest conservation: the total momentum added to the neighbors via the ∆ mb v a term exactly cancels that lost by the star, since its mass decreases by ∑ ∆ mb = m ej.  \n(vi) Account for PdV (mechanical) work: consider that we have a particle b representing a volume domain with mass mb around our source, with some mean distance in the volume element | x ba | (which we call the 'coupling radius'). The ejecta, in order to reach that point, must sweep up the mass mb (in e.g. a shock or shell) -  \nit cannot simply 'spread uniformly' throughout the volume. This means some PdV work must have been done, converting thermal energy into kinetic energy. Thus the correct momentum ( ∆ p '' b ) to couple into the domain b is not the initial ejecta momentum ∆ p ' b . Rather, if the shock is energy-conserving (neglecting second-order terms in the ratio of particle velocity to ejecta velocity, discussed in Paper II), it is trivial to show that the correct momentum is ∆ p ' b ( 1 + mb / ∆ mb ) 1 / 2 . In the early stages of SNe expansion, the shocks are indeed energy-conserving to high accuracy. Of course, at sufficiently long times (or equivalently large radii and/or large entrained masses), the shock becomes radiative, the residual thermal energy is lost, and the shock asymptotically reaches a final 'terminal momentum' p t (and becomes momentum, rather than energy conserving). Therefore we must impose an upper limit ∆ p ' b p t / p ej. We therefore have:  \n∆ p '' b ≡ ∆ p ' b MIN [√ 1 + mb ∆ mb , p t p ej ] (D22)  \np t M glyph[circledot] kms -1 ≈ 4 . 8 × 10 5 ( E tot , ej 10 51 erg ) 13 14 ( nb cm -3 ) -1 7 f ( Zb ) 3 2 (D23)  \nf ( Z ) = { 2 ( Z / Z glyph[circledot] < 0 . 01 ) ( Z / Z glyph[circledot] ) -0 . 14 ( 0 . 01 ≤ Z / Z glyph[circledot] ) (D24)  \nThe expression for p t comes from high-resolution simulations of individual SNe explosions (see e.g. Cioffi et al. 1988; Draine & Woods 1991; Slavin & Cox 1992; Thornton et al. 1998; Martizzi et al. 2015; Walch & Naab 2015; Kim & Ostriker 2015; Haid et al. 2016; Iffrig & Hennebelle 2015; Hu et al. 2016; Li et al. 2015), in media with different densities and metallicities. We discuss this at length in Paper II, and show that (1) it is the correct expression for a single SN explosion in a homogeneous background, given the same cooling functions and all other physics implemented in FIRE, (2) it appears to be remarkably robust, across many numerical studies, and (3) our conclusions are robust to variations in the exact value of p t much larger than its actual physical uncertainty. It is easy to verify, given the form of Eq. D22, that at sufficiently high resolution ( mb glyph[lessmuch] 1000 M glyph[circledot] ), the p t term simply never enters our equations - in other words, the SNe cooling radii are always resolved. This motivates our SNe explosion resolution criteria in the text. However the design of the expressions here is such that our coupling scheme automatically correctly treats each of e.g. the ejecta free-streaming, Sedov-Taylor, and snowplow phases. Note that since ∆ E represents the total energy, this is not directly modified by changing ∆ p (the correct thermal-kinetic breakdown is automatic).  \nAlso note that Eq. D22 is an approximation if the gas surrounding the star particle is moving at a non-uniform velocity (with nonnegligible velocities relative to the ejecta); in this limit the exact expression is given in Paper II (Appendix E).  \n(vii) Add final fluxes to the neighboring gas elements, in a fullyconservative manner:  \nm new b = mb +∆ m ' b (D25)  \n( Zmb ) new = Z new m new b =( Zmb ) +∆ m ' Z , b (D26)  \np new b = m new b v new b = p b +∆ p '' b (D27)  \nE new b = E new kinetic + U new internal = Eb +∆ E ' b (D28)  \nSo (like our hydrodynamic update), we add conserved quantities ( m , p , E ) and from those update primitive quantities ( Z , v , internal energy U , etc.). We check that any residual momentum or mass (from e.g. round-off error) is re-assigned to the star so conservation is always machine-accurate.  \nTo be fully consistent with the radiative losses described above (when the cooling radius is un-resolved), we must also modify U new internal . Following Thornton et al. (1998), the thermal post-shock energy outside R cool decays rapidly as ∝ ( r / R cool ) -6 . 5 ; so we estimate the effective R cool via the requirement that, at the end of the energy-conserving phase, ( 1 / 2 )( m ej + m swept[ R cool] ) v 2 f = ( 1 / 2 ) m ej v 2 ej and p t = m swept[ R cool] vf (where m swept is the enclose mass 'swept up' by the shell), giving R cool ≈ 28 . 4pc ( nb / cm -3 ) -3 / 7 ( E tot , ej / 10 51 erg ) 2 / 7 f ( Zb ) for p t in Eq. D23. If | x ba | < R cool, we leave U new internal un-modified. If | x ba | > R cool, we calculate the increase in internal energy from shockheating, ignoring cooling: ∆ Ub ≡ U new internal -U old internal , and then modify it to determine the correct internal energy: U internal = U old internal + ∆ Ub ( | x ba | / R cool ) -6 . 5 . We show in Paper II that this extra step has a negligible effect, since (by definition), when we couple the ejecta to a size/mass scale larger than R cool, it will radiate its energy rapidly, so in practice we find that if we simply leave ∆ Ub un-modified, the residual energy is (correctly) radiated away in the next timestep. But for the sake of physical consistency and accuracy, we adopt the full expression here.", 'APPENDIX E: ALGORITHMIC IMPLEMENTATION OF RADIATIVE FEEDBACK': "Now we describe the implementation of radiative feedback, used for radiation pressure (in all wavebands, UV-through-IR), photoionization, and photo-electric heating. The algorithm here was first developed in Hopkins et al. (2012c), and the version used in FIRE1, which is almost exactly identical to that here, 36 was presented in Hopkins et al. (2014). As noted in the text, we for convenience denote the radiative transport algorithm as the 'LEBRON' (Locally Extincted Background Radiation in Optically-thin Networks) approximation.  \nWe emphasize that this is not the same as the algorithm used in some earlier work (e.g. Hopkins et al. 2011). That algorithm was developed for very specific simulations which followed only infrared multiple-scattering radiation pressure (ignoring single-scattering radiation pressure, photo-heating, SNe, and OB/AGB winds), with much higher resolution than the FIRE simulations here (following star formation down to protostellar cores with densities ∼ 10 6 cm -3 and size scales < 0 . 1pc).  \nIn Paper III, we discuss each piece of the FIRE radiative feedback algorithm in detail, and consider a suite of both idealized test problems and cosmological simulations, to test and validate each and show how they influence our predictions. We also compare to a set of radiation-hydrodynamics simulations using alternative approximations to the radiation-hydrodynamics equations, specifically the flux-limited diffusion (FLD), optically-thin variable Eddington tensor (OTVET), first-moment (M1), and full Monte Carlo methods. We show that these give similar conclusions, provided care is taken with the alternative methods to ensure the radiation pressure terms are not artificially suppressed.  \nThe complete algorithm is:  \n36 The only difference between the radiation algorithms in FIRE-1 and FIRE-2, as noted in the main text, is that in FIRE-2 we allow ionizing photons to propagate outside of the numerical domain boundaries. In FIRE1, for numerical convenience, their propagation was 'truncated' at these (large-scale) boundaries. The effects of this are negligible in our simulations here, because almost all the 'work' done by ionizing photons is on nearby gas.  \n(i) Determine Background Radiation (Source Luminosities): every timestep ∆ ta for each star particle a , we take the luminosity L a ν = Ψ a ν m ∗ , a as a function of the star particle age, metallicity, and mass, directly from the stellar evolution models. This follows the tabulation given in Appendix A, for each of five broad bands we follow: ionizing ( L ion, λ < 912Å), far-UV ( L FUV, 912Å <λ< 1550Å), near-UV ( L UV, 1550 <λ< 3600Å), optical/near-IR ( L Opt, 3600Å < λ < 3 µ ), and mid/far-IR ( L IR, λ > 3 µ ). The FUV band is used for photo-electric heating, while NUV/optical bands dominate the single-scattering radiation pressure (bolometric luminosity), and the IR band is reserved for light re-radiated by dust.  \n(ii) Locally Extinct: We now process the absorption/extinction in the vicinity of each source. Along a sightline, the optical depth seen by the source is τ a ν = κ ν Σ a column , where κ ν is the flux-weighted opacity in each band (given in Appendix A; these are calculated for the mean un-obscured spectrum in the stellar populations models, integrated over each band). Since we are interested in local extinction, we approximate Σ a column using the Sobolev approximation for the isotropic (angle-averaged) column density integrated outward from the source:  \nτ a ν = 〈 κ 〉 ν ×〈 Σ a , Sobolev column 〉 φ , θ (E1)  \n〈 κ 〉 ν ≡ ∫ band κ ν 〈 L ν 〉 unabsorbed d ν ∫ band 〈 L ν 〉 unabsorbed d ν (E2)  \n〈 Σ a , Sobolev column 〉 φ , θ ≡ ρ a [ ha + ρ a |∇ ρ a | ] (E3)  \nHere ρ a , ∇ ρ a , and ha are the gas density, density gradient and interelement spacing, evaluated at location x a (with the same algorithm as our usual hydrodynamics). The ρ a / |∇ ρ a | term accounts for the gas column integrated to infinity - it is exact for e.g. a density distribution which declines exponentially with distance from the source - while the ρ a ha term is just the column through the local cell. Since this is isotropic, the absorbed luminosity (in a narrow band) is just L a abs , ν = ( 1 -exp ( -τ a ν )) L a ν and the surviving un-absorbed luminosity is:  \nL a emergent , ν = exp ( -τ a ν ) L a ν (E4)  \nRecall our spectral templates include negligible 'initial' luminosity in the mid/far IR band: but we assume the luminosity absorbed by dust - that from the FUV, NUV, and optical/near-IR bands - is immediately re-radiated in the mid/far-IR band, giving an emergent IR luminosity:  \nL a IR = ∑ ν = FUV , UV , Opt L a abs , ν = ∑ ν = FUV , UV , Opt ( 1 -exp ( -τ a ν )) L a ν (E5)  \nFor the ionizing band, the opacity comes from neutral hydrogen, and we must jointly solve for the ionizing state and photon absorption; we therefore treat this separately using a simple Stromgren approximation. From Step (i) above (the stellar evolution models) we have ˙ N ion a , the rate of production of ionizing photons ( ∝ L ion). Now, we take all gas elements b in the vicinity of a and sort them by increasing distance | x ba | . Beginning with the closest, we test whether it is already ionized (either because Tb > 10 4 K, or because it is already tagged as a member of another HII region), and if so we move on to the next-closest particle. If it is not ionized, we calculate the ionization rate needed to fully ionize it as: ∆ ˙ Nb = N ( H ) b β ne , b (where N ( H ) b = XH mb /µ mp is the number of H atoms in b , β ≈ 3 × 10 -13 cm 3 s -1 is the recombination coefficient, and ne , b is the electron density assuming full ionization). If ∆ ˙ Nb ≤ ˙ N ion a , then particle b is tagged as being within an HII region,  \nand the photons are 'consumed,' so ˙ N ion a → ˙ N ion a -∆ ˙ Nb . We then proceed to the next particle and repeat. If we reach a particle which is not ionized but for which ∆ ˙ Nb > ˙ N ion a , we determine whether or not to ionize it randomly, with probability p = ˙ N ion a / ∆ ˙ Nb , and consume the remaining photons (guaranteeing the correct total mass is ionized, on average). Any particle tagged as 'within an HII region' is fully ionized and not allowed to cool to temperatures lower than < 10 4 K within that same timestep. If we reach the end of the local computational domain, or distance from the source where the optically-thin flux density falls below the meta-galactic ionizing background, we stop the iteration and the remaining photons are 'emergent' as in Eq. E4. Tests of this algorithm (both static but also dynamic tests of D-Type ionization front expansion) are also shown in Hu et al. (2017).  \n(iii) Account for Momentum of Locally-Absorbed Photons: Over a timestep ∆ t , the absorbed photons in these bands impart their single-scattering momentum to the surrounding gas, radially directed away from the star particle, with total momentum  \n∆ p = L a abs c ∆ t = ∆ t c ∑ ν L a abs , ν (E6)  \nThis momentum flux is distributed among the neighbors, directed radially away from the star particle, as described in Appendix D.  \n(iv) Transport the Locally-Extincted Radiation via an Opticallythin Network: We now have an 'emergent' spectrum after local attenuation around each star, L emergent , ν . Since we assume the absorption is dominated by the gas/dust local to the star, and the emission (from the star) is isotropic, the incident flux F b ν and photon energy density e b ν at a distant gas element b are just  \nF b ν = ∑ a F a ν , emergent = ∑ a L a emergent , ν 4 π | x b -x a | 2 x b -x a | x b -x a | (E7)  \ne b ν = ∑ a L a emergent , ν 4 π c | x b -x a | 2 (E8)  \nThis is identical in form to the equation for gravity, so is computed in the same pass in the gravity tree (we 'soften' the sources identical to how we soften gravity, in fact, to prevent a 1 / r 2 divergence and reflect the physical fact that each star particle really represents many stars distributed within the softening length). 37  \n(v) Calculate incident radiative acceleration from long-range fluxes: For a gas element b with effective face area Ab and mass mb (hence surface density Σ b ≡ mb / Ab ) seeing an incident flux F b ν ,  \n37 The softening kernel, following Hopkins (2015), is given by replacing Eqs. E7-E8 with F ν , emergent = ∑ a ( 1 / 4 π ) L a emergent , ν ( x b -x a ) H -3 a F s ( uba ) and glyph[epsilon1] b ν = ∑ a ( 1 / 4 π c ) L a emergent , ν | x b -x a | H -3 a F s ( uba ) where uba ≡| x b -x a | / Ha , Ha =( 24 /π ) 1 / 3 ha =( 24 /π ) 1 / 3 glyph[epsilon1] a is the maximum kernel search radius, and  \nF s ( u ) ≡                        32 15 [ 5 + 3 u 2 ( 5 u -6 ) ] ( u ≤ 1 2 ) 32 15 [ 10 -45 2 u + 18 u 2 -5 u 3 -1 32 u 3 ] ( 1 2 < u < 1 ) 1 u 3 ( u ≥ 1) (E9)  \nThis becomes exactly inverse-square at r > H , and but prevents a 1 / r 2 divergence as r → 0 (with flux ∼ 1 / h 2 a instead of 1 / r 2 ).  \nthe exact radiative acceleration is given by  \nd v b dt ∣ ∣ ∣ ν = 1 mb d p b dt ∣ ∣ ∣ ν = F b ν c Ab mb [ 1 -exp ( -κ ν mb Ab )] (E10)  \nIn the optically thin limit ( κ ν Σ b glyph[lessmuch] 1), this reduces to the common expression a = κ ν F ν / c , but in the optically thick limit ( κ ν Σ b glyph[greatermuch] 1), the force saturates at mb a =( F ν Ab ) / c , i.e. the element absorbs all the flux across its subtended area (but no more). For simplicity here we take the effective area to be that of a sphere with the same volume as the element ( = mb /ρ b =( 4 π/ 3 ) h 3 b ), i.e. π h 2 b ; using the more complicated hydrodynamic face areas introduces negligible ( glyph[lessmuch] 10%) differences in the accelerations here. We adopt κ ν = 〈 κ 〉 ν for each band.  \n(vi) Self-shield and pass incident fluxes to cooling routines: Because we have accounted for shielding around the emitter, but not the absorber, we include an additional shielding pass at absorption for the photo-heating terms: at a gas element b , we take the photon energy density 〈 e b ν 〉 = e b ν exp ( -τ b ν ) , where e b ν is the photon energy density given by Eq. E8, and τ b ν is the optical depth estimated using the Sobolev approximation in Eq. E1, but now at the location of the absorbing gas element (instead of around the emitting star). The resulting, shielded e b ν are passed to the cooling/heating routines, to compute photo-ionization and photo-electric heating rates as in Appendix B.", 'APPENDIX F: ADDITIONAL FLUID PHYSICS: MAGNETIC FIELDS, CONDUCTION, VISCOSITY, TURBULENT DIFFUSION': "Here, we describe the numerical implementations of additional physics not included in the 'core physics only' FIRE simulations, but studied either here or in companion papers (e.g. Su et al. 2017) which take standard FIRE-2 simulations and add, e.g. magnetic fields. We emphasize again that these physics are not used in the 'default' or 'core physics' runs in this paper. However, because we wish to present a complete, thorough, and fully-consistent set of numerical methods, we summarize them here, referring to the appropriate methods papers for more details.", 'F1 Magnetic Fields': 'In simulations with magnetic fields, we solve the equations of ideal magnetohydrodynamics (MHD) as implemented in GIZMO in Hopkins & Raives (2016). The exact numerical formulation of the equations is presented there along with an extensive series of several dozen test problems, as well as tests of full galaxy simulations using our FIRE physics. The tests demonstrate that the implementation in our MFM solver is accurate and converges at second order. In particular, Hopkins & Raives (2016), Hopkins (2016), and Hopkins & Lee (2016) show that our implementation correctly captures traditionally difficult phenomena such as the magneto-rotational instability (MRI), magnetic jet launching in disks, magnetic fluidmixing instabilities, and sub-sonic and super-sonic MHD turbulent dynamos. The accuracy and convergence order appears comparable to state-of-the-art grid codes (e.g. ATHENA) on the problems of interest and greatly superior to the P-SPH implementation in GIZMO, especially in problems where angular momentum, super-sonic advection, strong shocks, and fluid mixing instabilities appear (typical of cosmological simulations). Non-ideal MHD effects (Ohmic resistivity, ambipolar diffusion, and the Hall effect) are also implemented and well-tested in GIZMO (Hopkins 2017), but these are not expected to be important on galactic scales.  \nReaders interested in more details of our numerical implementation of MHD should consult Hopkins & Raives (2016) and the public GIZMO source code.', 'F2 Anisotropic Spitzer-Braginskii Viscosity and Conduction': "The implementation of anisotropic Spitzer-Braginskii viscosity and conduction in GIZMO is described and tested in Hopkins (2017). In addition to the usual MHD fluxes, this adds an anisotropic viscous stress-energy tensor Π to the momentum flux ( Fp = Π ) and energy flux ( F e = Π · v ), and a conductive energy flux F e = K · ∇ T . The appropriate anisotropic tensor expressions for MHD are given by (Spitzer & Härm 1953; Braginskii 1965):  \nK ≡ κ cond ˆ B ⊗ ˆ B (F1)  \nκ cond = 0 . 96 fi ( kBT ) 5 / 2 kB m 1 / 2 e e 4 ln Λ ( 1 + 4 . 2 glyph[lscript] e /glyph[lscript] T ) -1 (F2)  \nΠ ≡ 3 ν visc ( ˆ B ⊗ ˆ B -1 3 I )[( ˆ B ⊗ ˆ B -1 3 I ) : ( ∇⊗ v ) ] (F3)  \nν visc = 0 . 406 fi m 1 / 2 i ( kBT ) 5 / 2 ( Zi e ) 4 ln Λ ( 1 + 4 . 2 glyph[lscript] i /glyph[lscript] v ) -1 (F4)  \nwhere ⊗ denotes the outer product; ˆ B is the direction of the magnetic field vector; I is the identity matrix; v the velocity; ':' denotes the double-dot-product ( A : B ≡ Trace ( A · B ) ); ln Λ ≈ 37 . 8 is the Coulomb logarithm (Sarazin 1988); me , e , mi , Zi e are the electron mass and charge and ion mass and charge; fi the ionized fraction (calculated self-consistently in our cooling routines); kB the Boltzmann constant; glyph[lscript] e ( glyph[lscript] i ) is the electron (ion) mean-free path, and glyph[lscript] T = T / |∇ T | ( glyph[lscript] v = | v | / ||∇⊗ v || ) is the temperature (velocity) gradient scale length (this term correctly accounts for saturation of κ cond or ν visc when electrons/ions have long mean free paths, by not allowing the gradient scale length to be shorter than glyph[lscript] e , i ). In these equations, κ cond is the conductivity, and ν visc the viscosity. Additional details of the coefficients, and a study of their effects, are in Su et al. (2017). In Hopkins (2017), we show that the numerical implementation of these fluxes is accurate, able to handle arbitrarily large anisotropies, converges comparably to higher-order fixedgrid codes, and is able to correctly capture complicated non-linear instabilities sourced by anisotropic diffusion such as the magnetothermal and heat-flux bouyancy instabilities.  \nReaders interested in more details of our numerical implementation of anisotropic diffusion should consult Hopkins (2017) and the public GIZMO source code.", 'F3 Sub-Grid Turbulent Eddy Diffusivity': "In some models for turbulence (e.g. mixing-length theory), the effects of unresolved (small-scale) eddies and microphysical processes transporting passive scalars (such as metals) are treated as diffusive processes. The implementation and tests of generic diffusion operators in GIZMO are presented in Hopkins (2017); the solver is an explicit finite-volume scheme which converges at second-order accuracy (comparable to higher-order grid codes) and manifestly conserves metal mass. For a passive scalar, the transport equation is: ∂ ( ρ Z ) /∂ t = ∇· ( κ turb ρ ∇ Z ) , where Z is the abundance per unit mass of the scalar (i.e. the metallicity) and the 'eddy diffusivity' κ turb ∼ λ eddy v eddy is the product of the scale length and rms velocity of the largest un-resolved eddies (those at the resolution scale), which dominate the transport on unresolved scales (larger eddies are, of course, resolved). In other words, one simply assumes that the diffusion or mixing time at scale λ eddy scales with the eddy turnover time. Following the common Smagorinsky  \n(1963) approximation, we can approximate the 'eddy diffusivity' as  \nκ a turb ≡ √ 2 C 2 ‖ S a ‖ h 2 a (F5)  \nwhere C ∼ 0 . 05 -0 . 15 is a constant calibrated to numerical simulations, motivated by a Kolmogorov cascade in Smagorinsky (1963), ha is the grid scale (for our MFM method, this is equal to the rms inter-element separation), and S a ≡ [ ( ∇⊗ v ) a +( ∇⊗ v ) T a ] -Trace ( ∇⊗ v ) a / 3 is the symmetric shear tensor (and ‖ S ‖ denotes the Frobenius norm). Note that we use our higher-order matrixbased gradient formalism from Hopkins (2015) to calculate S a ; this is much more accurate and less noisy compared to common SPH or pure 'face-based' mesh gradient estimators, which is important to reduce artificial numerical diffusivity (see e.g. Maron & Howes 2003; Luo et al. 2008; Lanson & Vila 2008; Mocz et al. 2014; Pakmor et al. 2016).  \nThe key assumption here - namely, the assumption that the diffusion timescale scales with eddy turnover time - has been verified in many experiments on ISM scales (see Pan & Scannapieco 2010; Petit et al. 2015, and references therein), and the scaling in Eq. F5 has been used in many applications in galaxy simulations (e.g. Wadsley et al. 2008). In Colbrook et al. (2017), we have performed our own study of the turbulent mixing, using 3D, highresolution supersonic turbulent box simulations (with and without magnetic fields and/or shear), and verify that this prescription, with C ≈ 0 . 05, is reasonable specifically in our identical MFM code with the definitions of h and S here (although such simple prescriptions do fail to capture some potentially important non-Gaussian features which emerge from real, resolved turbulent mixing). An independent, more extensive study (including a range of more complex problem setups) will be presented in Rennehan et al. (in prep.), but also finds C ≈ 0 . 03 -0 . 05. We therefore adopt C = 0 . 05. However, in the main text (§ 7.2) and in Escala et al. (2018), we show that order-of-magnitude variation in C produces no significant effects on our predictions.  \nWe stress that a term like Eq. F5 is 'built into' many numerical hydrodynamic methods. Specifically, it is well known that in finite-volume methods with advective mass fluxes (e.g. traditional grid-based methods or moving-mesh codes), an intrinsic numerical diffusivity in advection with magnitude ∼ ha ∆ v ( ha ) appears; this automatically produces scalar/metal diffusion via 'numerical mixing.' It is straightforward to show that if the necessary assumptions of the Smagorinsky model (Eq. F5) are true, then the artificial numerical mixing in these methods is always larger than the 'true' effective turbulent diffusivity. If we used such a method, it would therefore not be necessary to explicitly solve Eq. F5. However, in our default MFM hydrodynamic method, we follow fixed-mass elements (i.e. there are no advective mass fluxes, hence no artificial 'numerical diffusivity' of passive scalars). While the methods will converge to an identical solution at sufficiently high resolution (Hopkins 2015), the concern is that at fixed resolution, MFM will under-estimate the metal-mixing owing to un-resolved smallscale eddies that should mix between the boundaries of neighboring resolution elements. If we resolved individual stars, the stars would draw mass from many resolution elements (each with their own abundances) and this would still not be a problem, but at our resolution single star particles inherit the abundances of their (single) parent gas particle, so this effect can artificially introduce 'shot noise' in the abundances of stars forming from neighboring gas elements if we do not include an explicit numerical mixing term. One therefore can view Eq. F5 as a purely numerical term which 're-  \nes' the desirable aspect of the numerical diffusivity present in certain numerical methods.  \nWe caution, however, that simple diffusion prescriptions such as Eq. F5, naively applied, can substantially over-estimate the diffusivity. The critical assumption is that the resolution-scale motion ∆ v ∼‖ S ‖ h ∼ v eddy is entirely due to turbulence; if there is any bulk motion included in S , this will over-estimate κ turb. This can be particularly problematic if e.g. differential rotation in a disk or shear in CGM outflows is poorly-resolved, in which case the naively inferred κ turb can over-estimate by an order of magnitude the true turbulent motion. For example, if the disk scale height H disk ≈ σ turb / Ω is unresolved, ha glyph[greatermuch] H disk, then Eq. F5 will return κ turb ∼ ha ( ha Ω) instead of the correct maximum diffusivity for disk-scale eddies, κ turb ∼ H disk σ turb ∼ H 2 disk Ω . There is no obvious universal 'switch' to cure these pathologies; however we can limit the magnitude of the errors. In Hopkins (2015) we develop two mesh-free finiteelement hydrodynamic methods, our default MFM method here, and a 'meshless finite-volume' (MFV) method, the latter of which includes advective mass fluxes (more similar to a moving-mesh code), but otherwise is identical to MFM. The MFV method therefore includes the inherent numerical diffusivity described above. In the diffusion step, therefore, we can first calculate the absolute value of the metal flux that would have been calculated in MFV (owing simply to the advection term; see Hopkins 2015 for the exact values of these terms), and then impose this as an upper limit to the diffusive flux. Since MFV is a second-order, quasi-Lagrangian method, this eliminates the most egregious errors in Eq. F5. We find that in idealized test problems, this correction is negligible, but in realistic cosmological simulations it prevents the most severe pathological situations. Essentially, then, our implementation of 'unresolved turbulent diffusion' (Eq. F5) is guaranteed to at most - produce the same metal mixing we would have obtained had we simply run our simulations using a finite-volume (MFV or moving-mesh) hydrodynamic method. A more detailed study of various (more sophisticated) turbulent and numerical mixing models will be presented in Rennehan et al. (in prep.); preliminary results indicate that the alternative methods give similar results in galaxy-scale simulations.  \nReaders interested in more details of our numerical implementation of eddy diffusivity should consult Hopkins (2017) and the public GIZMO source code.", 'APPENDIX G: COMPUTATIONAL SCALING & RUNTIME REQUIREMENTS': "The work here was made possible by extensive optimizations and improvements to the code scaling for 'zoom-in' simulations. The challenge in these high-dynamic range problems is that small, dense regions (e.g. dense GMCs or star clusters), which occupy an extremely small fraction of the total mass and volume of the simulation (and therefore cannot be 'broken up' over too many processors) require extremely small timesteps. But, especially with strong feedback present, the rest of the simulation (which is free to take much larger timesteps) cannot advance until they 'catch up.' We have addressed this with several improvements.  \n(i) We have made optimizations to the feedback sub-routines (for example, long-range radiation forces) which require significant neighbor communication but can be efficiently included in other operations such as tree construction, reducing their portion of the runtime from ≈ 40% in previously-published work to < 10% in the current code. More generally, the entire GIZMO code has been lineby-line optimized, manually unrolling or (where possible) vector-  \n]  \nFigure G1. Code scalings of GIZMO in full production-quality FIRE-2 simulations, at our production resolution, with the full, identical physics of gravity, hydrodynamics, cooling, star formation, and feedback to our published simulations. Left: Strong scaling for a zoom-in of a MW-mass ( m12i ) or dwarf galaxy ( m10q ) halo, each using 1 . 5 × 10 8 particles, run to 25% of the age of the Universe, using the optimal MPI+OpenMP hybrid configuration at each core number. Our optimizations allow us to maintain near-ideal strong scaling to ∼ 14 , 000 cores per billion particles (2048 for the problem shown). Center: Weak scaling, for a full cosmological box, populated with the same high-resolution particles, run for a short fraction of the age of the Universe ( z ∼ 10). Here, we keep resolution fixed at baryonic particle mass 7000 M glyph[circledot] , identical to our high-resolution MW simulation at left , but we increase the cosmological volume from 2 to 10 4 comoving Mpc 3 . The weak scaling of GIZMO's gravity+MHD algorithm is near-ideal (actually slightly better at intermediate volume, owing to fixed overheads and increasing statistical homogeneity of the box at larger sizes), to greater than a million threads (here 2 20 threads, 2 18 cores, 2 17 MPI tasks, 2 14 nodes). Right: 'Weak scaling' test, increasing the resolution instead of the problem size (specifically, increasing the particle number for the same MW-mass galaxy). Because the resolution increases (hence timestep decreases) with particle number here, the ideal weak scaling for a converged solution is wall-clock time ∝ t Hubble / ∆ t ∝ (∆ x / cs ) -1 ∝ N 1 / 3 , shown. Our achieved scaling is only slightly worse than this, because new, dense structures such as star clusters appear at higher resolution.  \n<!-- image -->  \nizing certain expensive operations, replacing expensive functions in neighbor loops with look-up tables, eliminating redundant operations, and pre-computing additional quantities outside of loops; this has produced an additional factor ∼ 2 speed improvement.  \n(ii) We have optimized the structure of the domain decomposition to make it more spatially flexible and separately parallelize each level of the timestep hierarchy (increasing the memory imbalances by factors of ∼ 2 -3, but extending the strong scaling to ∼ 2 times as many cores at fixed resolution). We have also more aggressively implemented a problem-specific particle weighting scheme, where e.g. dense, star-forming gas, and young stars (as opposed to old stars) are given larger weights in the domain decomposition so that their future cost (via gravitational collapse and/or feedback) is more accurately predicted. This allows for a further factor of ∼ 2 -3 reduction in load imbalances at the smallest timesteps.  \n(iii) We employ a hybrid tree-particle mesh gravity solver, following GADGET-3, to efficiently reduce the cost of the gravity solution for the low-resolution regions outside of the zoom-in area.  \n(iv) We use the adaptive individual-timestep integration scheme from GADGET-3 with a hierarchical power-of-two subdivision, updated such that in each timestep we calculate pairwise updates of all fluxes of conserved quantities at interfaces (maintaining exact conservation, and eliminating all redundant pair-wise interations; see Springel 2010).  \n(v) We adopt adaptive gravitational softenings as described above. This imposes essentially no cost for gas (since hi must be computed already for hydrodynamics), but allows us to take much larger timesteps for low-density particles, and (more important for this problem) avoid over-softening particles in dense regions (where a too-large softening radius might encompass thousands of neighbor particles, which imposes a substantial cost in the tree-gravity calculation).  \n(vi) We have developed and use a hybrid OpenMP-MPI parallelization of the code, which allows us to extend the weak scaling of the code considerably further as we go to large processor number (where we previously found communication costs between neighbors, which are alleviated by the shared-memory structure of OpenMP, were beginning to dominate the run time).  \n(vii) Our MFM hydrodynamic solver is actually slightly faster than the P-SPH algorithm used for FIRE-1. Hopkins (2015) compare run-times on three problems: a 3D Kelvin-Helmholtz instability, an isolated (non-cosmological) galaxy with star formation and feedback, and the cosmological but strictly ideal-gas (no star formation or cooling) Santa Barbara cluster test. They show that compared to the P-SPH formulation of SPH from Hopkins (2013a) (which incorporated improvements in artificial diffusion terms as well as a larger neighbor number needed in SPH to capture certain instabilities), the speedup on these tests ranged from a factor ∼ 1 . 3 -2 . 5, mostly owing to the larger neighbor number needed in P-SPH to achieve comparable accuracy.  \n(viii) We also stress that the inclusion of realistic feedback itself greatly speeds up the code (for galaxy formation simulations) - perhaps more than any purely numerical optimization. Dense regions which would otherwise slow down the computation (as described above) tend to be quickly destroyed by stellar feedback. Without feedback, it would be impossible to run the simulations here with the same resolution and cooling physics below redshift z ∼ 2, because the extremely dense relic star clusters would require constant, extremely short timesteps.  \nReaders interested in more details should consult the public source code.  \nFigure G1 demonstrates the scaling of the code GIZMO on a production quality set of FIRE-2 simulations, including all the physics of our production runs (full cosmological integration with self-gravity, baryonic physics including cooling, star formation, and stellar feedback, etc.). This is a 'real world' test, as opposed to the scaling on idealized test problems (which can, of course, be much better). All runs were run with an otherwise identical version of the code; at each CPU number the optimal OpenMP-MPI configuration was used. Strong and weak-resolution scaling tests were run on the XSEDE Stampede machine, weak-problem size runs on the DOE ALCF Mira machine.  \nOur optimizations allow us to extend good strong scaling, at our modest 'typical FIRE-1' resolution, to ∼ 1024 cores. Even more strikingly, the optimizations we have made allow us to maintain good weak scaling up to ∼ 4096 cores for a simulation with 3 × 10 8 particles - our 'Latte' resolution - and ∼ 16 , 384 cores  \nFigure H1. Believeland.  \n<!-- image -->  \nfor a zoom-in simulation with > 10 9 particles (and as many as ∼ 10 6 cores for large-volume simulations with ∼ 10 10 particles, which naturally exhibit superior weak scaling). This is especially non-trivial for these sorts of problems, since the gravitational softening is adaptive, so higher particle number implies smaller force softening and hence smaller timesteps in the dense regions.  \nFor comparison, if we compare similar runs to the strongscaling tests in Fig. G1 with the public GADGET-2 code (using the simpler Springel et al. (2005) sub-grid model for ISM physics and feedback, with wind mass-loading set by-hand to produce a similar mass as our GIZMO runs), we find the scaling saturates at ∼ 64 -128 cores.", 'APPENDIX H: CONGRATULATIONS!': 'You made it to the end of this paper. Fig. H1 is your reward. Go Cleveland!'}