FIRE-2 Simulations: Physics versus Numerics in Galaxy Formation

TL;DR

FIRE-2 simulations incorporate advanced numerics and physics to study galaxy formation, showing that many galaxy properties are robust to numerical methods, with some central mass concentrations sensitive to these details.

Contribution

This paper introduces FIRE-2, an improved simulation framework with updated numerics and physics, and systematically compares it to FIRE-1 to assess robustness of galaxy formation results.

Findings

Most galaxy properties are robust to numerics if certain criteria are met.

Central mass concentrations in massive galaxies are sensitive to numerical details.

Multiple feedback mechanisms are crucial for realistic galaxy evolution.

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.

Tables

Table 1. Table 1: Initial suite of FIRE-2 simulations run to redshift z = 0 𝑧 0 z=0

	Parameters 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.
(2) $M_{\mathrm{halo}}^{\mathrm{vir}}$: Virial mass (following Bryan & Norman, 1998) of the “target” halo at $z=0$ in simulation with baryons.
(3) $R_{\mathrm{vir}}$: Virial radius (spherical) of the main halo at $z=0$.
(4) $M_{\ast }$: Stellar mass (within $\le 3R_{1/2}$) of the central galaxy in the target halo at $z=0$.
(5) $R_{1/2}$: Half-mass radius of stars in the central galaxy at $z=0$ (see § 3).
(6) $m_{i,\mathrm{\hspace{0.17em}1000}}$: the mass of baryonic (gas or star) particles, in units of $1000M_{☉}$. Dark matter (DM) particle masses are $\approx 5\times$ larger, according to the universal baryon fraction.
(7) ${ϵ}_{\mathrm{gas}}^{\mathrm{MIN}}$: Minimum gravitational force softening reached by gas in the simulation. For gas, gravitational force resolution is identical to hydrodynamic spatial resolution (the same spatial gas distribution appears in gravity and hydrodynamic equations at all times). The gravitational force softening, ${ϵ}_i$, therefore corresponds to the gas inter-particle separation, $h_i$: ${ϵ}_i=h_i$. See § 4 for other definitions of “spatial resolution”. Forces become exactly Keplerian (point-mass like) at $>1.95{ϵ}_{\mathrm{gas}}$ from a particle; the “Plummer-equivalent” softening is $\approx 0.7{ϵ}_{\mathrm{gas}}$.
(8) $r_{\mathrm{DM}}^{\mathrm{conv}}$: 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 ${ϵ}_{\mathrm{DM}}$ is much less important, as long as it is $\lesssim r_{\mathrm{DM}}^{\mathrm{conv}}$. In our default runs, ${ϵ}_{\mathrm{DM}}$ is fixed with values $=40\mathrm{pc}$ for our dwarfs and $=30\mathrm{pc}$ in our MW-mass (m12) systems.
(9) Notes: Additional information on each simulation.

Table 2. Table 2: Physics & Numerics Explored in This Paper (and Papers II &  III )

	A 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.
(1) Physics/Numerics: what we consider.
(2) §: Section where we pursue a detailed study of the effects of each numerics/physics on galaxy formation.
(3) Effects in FIRE-2 Simulations: Overall summary of the effects of variation in the relevant physics or numerics, insofar as it is relevant (or not) for the predictions of our simulations. This applies only for quantities discussed in this paper, that is, global galaxy properties such as SFRs, stellar masses, sizes, and morphologies. For example, although we show that arbitrarily removing molecular chemistry from our cooling networks has no effect on predicted galaxy properties or star formation (because other cooling channels are available and molecular gas is primarily a tracer, not a causal driver of star formation), molecular chemistry is obviously fundamentally important if one wishes to predict molecular lines. Furthermore, we do not examine detailed properties of the CGM or IGM, where different physics may dominate.
(4) Guidelines: Approximate “rules of thumb” for the relevant physics or numerics in the context of our “zoom-in” galaxy simulations. In the text, we provide more detailed guidelines. For example, for numerical resolution and other numerical parameters, we provide equations that approximately determine whether or not key physics should be resolved.

Table 3. Table 3: Illustrative Examples of “Resolution” Scales in a MW-Mass Halo ( m12i )

Several “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.

Equations

$$m_{i,1000}$$

$$m_{\mathrm{DM}}$$

$$M_{\mathrm{halo}}^{\min }$$

$$M_{\ast }^{\min }=M_{\mathrm{cloud}}^{\min }$$

$$m_{i,1000}^{\mathrm{Toomre}}\lesssim \frac{100}{N_{10}}\frac{M_{\mathrm{Toomre}}}{10^6{M}_{\text{sun}}}\sim 100\left(\frac{f_{\mathrm{gas}}^2{M}_{\mathrm{gas}}}{10^6{M}_{\text{sun}}{N}_{10}}\right)$$

$$m_{i,1000}^{\mathrm{SNe}}\lesssim 0.9(N_{\mathrm{NGB}}/32{)}^{-1/3}{N}_{\mathrm{SNe}}$$

$$m_{i,1000}^{\mathrm{burstiness}}\lesssim 5\left(\frac{M_{\mathrm{halo}}}{10^{10}{M}_{\text{sun}}}\right)$$

$$h_i^{\mathrm{gas}}$$

$${ϵ}_i^{\mathrm{gas}}\equiv h_i^{\mathrm{gas}}$$

$$\frac{Q_{\mathrm{heat}}^{gas-gas}}{\mathrm{erg}\mathrm{c}{\mathrm{m}}^3{\mathrm{s}}^{-1}}$$

$$\frac{Q_{\mathrm{heat}}^{gas-DM}}{\mathrm{erg}\mathrm{c}{\mathrm{m}}^3{\mathrm{s}}^{-1}}$$

$$\Delta t_{\ast }<\mathrm{MAX}\left(10^4\mathrm{yr},\frac{t_{\ast }}{300}\right)$$

$$\Delta t_{\mathrm{AGS}}^a$$

$$v_{\mathrm{sig},\mathrm{AGS}}^a$$

$$⟨\tilde{\nabla }\cdot \mathbf{v}{⟩}_a$$

$$\frac{d{E}_a^{\mathrm{SPH}}}{dt}=$$

$${\tilde{W}}_{ba}\equiv$$

$${\tilde{v}}_{ab}^{\mathrm{sig}}\equiv$$

$${\alpha }_{ab}\equiv$$

$$\Theta ({\tilde{v}}_{ab}^{\mathrm{sig}})\equiv$$

$$\frac{d{e}_{\mathrm{thermal}}^a}{dt}=-n_{\mathrm{H},a}^2{\Lambda }_{\mathrm{net}}^a=-n_{\mathrm{H},a}^2\sum _i{\Lambda }_i$$

$${\Lambda }_{\mathrm{CE}}=$$

$${\beta }_{H0}$$

$${\beta }_{\mathrm{HeI}}$$

$${\tau }_5$$

$${\Lambda }_{\mathrm{CI}}=$$

$$\left(2.18{\gamma }_{H0}{\tilde{n}}_{H0}+3.94{\gamma }_{He0}{\tilde{n}}_{He0}+8.72{\gamma }_{\mathrm{HeI}}{\tilde{n}}_{\mathrm{HeI}}\right)$$

$${\gamma }_{H0}$$

$${\gamma }_{He0}$$

$${\gamma }_{\mathrm{HeI}}$$

$${\Lambda }_{\mathrm{Rec}}=$$

$$\left({\alpha }_{\mathrm{HI}}{\tilde{n}}_{\mathrm{HI}}+\left[{\alpha }_{\mathrm{HeI}}+\frac{629922.78}{T}{\alpha }_{\mathrm{di}}\right]{\tilde{n}}_{\mathrm{HeI}}+{\alpha }_{\mathrm{HeII}}{\tilde{n}}_{\mathrm{HeII}}\right)$$

$${\alpha }_{\mathrm{HI}}$$

$${\left(1+\frac{T^{0.5}}{1.774}\right)}^{-0.252}{\left(1+\frac{T^{0.5}}{838.81}\right)}^{-1.748}$$

$${\alpha }_{\mathrm{HeI}}$$

$${\left(1+\frac{T^{0.5}}{0.2065}\right)}^{-0.2108}{\left(1+\frac{T^{0.5}}{6063.0}\right)}^{-1.7892}$$

$${\alpha }_{\mathrm{HeII}}$$

$${\left(1+\frac{T^{0.5}}{2.5092}\right)}^{-0.252}{\left(1+\frac{T^{0.5}}{1677.6}\right)}^{-1.748}$$

$${\alpha }_{\mathrm{di}}$$

$$\left(1+0.3\exp \left[-\frac{94000}{T}\right]\right)$$

$${\Lambda }_{\mathrm{FF}}=$$

$${\beta }_{\mathrm{ff}}(T)$$

$$\left[1.1+0.34\exp \left\{-(5.5-{\log }_{10}[T]{)}^2/3\right\}\right]$$

$${\Lambda }_{\mathrm{Metal}}=$$

Full text

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 of 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).

However, the observed Kennicutt-Schmidt (KS) relation implies that gas consumption timescales are long ($\sim 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 incorporate 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 $\dot{M}_{\rm 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 (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).

Until 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., 2011; 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 CGM-scale properties that depend explicitly on e.g. the phase-structure of feedback-driven outflows (see Hummels et al., 2013).

Accurate treatment of star formation and galactic winds ultimately 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 that 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, b, 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.

Motivated by the success and predictive power of these simulations, in Hopkins et al. (2014) we introduced the Feedback In Realistic Environments (FIRE) project.111See the FIRE project website:

http://fire.northwestern.edu

For additional movies and images of FIRE simulations, see:

http://www.tapir.caltech.edu/~phopkins/Site/animations 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\sim 10-10^{10}\,$K, star formation restricted only to self-gravitating, self-shielding, molecular, high density ($n_{\rm H}\gtrsim 5-50\,{\rm cm^{-3}}$) gas, resolution reaching $\sim 250\,M_{\sun}$ or $\sim 0.5\,$pc, and (most importantly) explicit treatment of stellar feedback including the energy, momentum, mass, and metal fluxes from SNe Types Ia & II, stellar mass-loss (O-star and AGB), radiation pressure (UV and IR), and photo-ionization and photo-electric heating. All stellar evolution and feedback inputs are taken directly from stellar evolution models, without subsequent “parameter tuning.”

In 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.

For 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, cosmic 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).

We 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).

The 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).

Table 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 of 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.

Before 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, 2012b, 2012a, 2013a, 2013b, 2013c, 2013e, 2013d). These developments and improvements were then synthesized into the physics implemented in FIRE-1, described in Hopkins et al. (2014).

In 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.

As 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 smoothed-particle 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).

We therefore employ the meshless finite-mass (MFM) magneto-hydrodynamics solver in GIZMO.222A public version of this code is available at http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html 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 has been tested extensively (involving $\sim 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), moving-mesh codes (e.g. AREPO; Springel 2010), and “modern” SPH methods (e.g. P-SPH; Hopkins 2013a; Rosswog 2015; Hu et al. 2014).

We 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).

As 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.

For 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.

Readers 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.

By 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.333Hopkins (2015) show this is the leading-order accurate expression for the potential if we integrate Poisson’s equation using the exact differential mass distribution assumed in the hydrodynamic equations. Because the separations $h_{i}^{\rm gas}$ change, we must be careful to maintain energy and momentum conservation correctly when elements interact inside the softening; Price & Monaghan (2007) show how the relevant expressions can be rigorously derived from the softened $N$-body particle Lagrangian and the expressions used in GIZMO are presented in Hopkins (2015), Appendix H2. This means that, for gas, the gravitational force softening resolution $\epsilon_{i}^{\rm gas}$ identically follows the inter-particle/cell separation $h_{i}^{\rm gas}=\Delta x$, where $\Delta x$ is the equivalent cell-length in a fixed, Cartesian mesh (i.e. the average distance between particle centers, around element $i$, is $h_{i}^{\rm gas}$), so $\epsilon_{i}^{\rm gas}\equiv h_{i}^{\rm gas}=16\,{\rm pc}\,m_{i,\,1000}^{1/3}\,(n_{\rm H}/10\,{\rm cm^{-3}})^{-1/3}$ (where $m_{i,\,1000}$ is the mass resolution in units of $1000\,M_{\sun}$). Note this definition is independent of the exact softening kernel shape or “neighbor number.”

As 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 self-gravity 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.

For 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 $\epsilon_{i}^{\rm DM}$ 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 $\epsilon_{i}^{\ast}$ 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).

Readers 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, metal-line, 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).

With 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 radiative 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 Faucher-Giguè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 $30\,$K. Cosmic ray heating follows Guo & Oh (2008) accounting for both hadronic and Compton interactions, with a uniform cosmic ray background of $\sim 5\,{\rm eV\,cm^{-3}}$. At very high densities ($\sim 10^{10}\,{\rm 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 directly into our fully-implicit solution. A $10\,$K temperature floor is enforced, but has no detectable effect on our results.

In 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 sink-particle 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:

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), $\alpha\equiv(\delta v^{2}+c_{s}^{2})\,\delta r/G\,m_{\rm gas}(<\delta r)=[\|\nabla\otimes{\bf v}\|_{i}^{2}+(c_{s,\,i}/h_{i})^{2}]/(8\pi\,G\,\rho_{i})<1$, where $\delta v=\|\nabla\otimes{\bf v}\|_{i}\,h_{i}$ and $c_{s}$ give the kinetic and thermal energy, respectively, within the resolution scale $\delta r\rightarrow h_{i}$ around the particle ($\otimes$ 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 this 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 self-gravitating 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. 2. 2.

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 ($\rho_{\rm mol}=f_{\rm mol}\,\rho$). This is specifically a requirement that the gas be self-shielding, and therefore able to cool to low temperatures. Given the high $n_{\rm crit}$ (see below), this criterion typically has negligible impact, since the high-density gas is typically all shielded and molecular anyways. 3. 3.

Jeans Unstable: We require a thermal Jeans mass below the maximum of the particle mass or $10^{3}\,M_{\sun}$ 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_{\sun}$ 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. 4. 4.

Dense: To prevent spurious application of the criteria above, we also check that $n_{\rm H}>n_{\rm crit}=1000\,{\rm cm^{-3}}$ (much larger than the mean galaxy density $\langle n\rangle$). This restricts star formation to dense molecular clouds fragmenting out of the background disk.

If gas meets all the criteria above, we assume it turns into stars at a rate $\dot{\rho}_{\ast}=\rho_{\rm mol}/t_{\rm ff}$ where $t_{\rm 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; Guszejnov & Hopkins, 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 $\sim 1-10\%$ per free-fall time (see Hopkins et al. 2011, 2013e, 2014) even in gas with densities $\gtrsim 100\,{\rm cm^{-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.

We 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_{\ast}=t-t_{\rm form}$, metallicity (inherited from its progenitor 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).444Of 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).

Here, 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_{\sun}$ 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.

The 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.

Supernovae (Ia & II): (For details, see Appendix D.) Every timestep $\Delta 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 $\Delta t$; our mass and time resolution is such that $p\ll 1$, i.e. we explicitly treat individual SNe explosions, rather than model their collective effects indirectly. We determine probabilistically if an event occurs within $\Delta 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_{\sun}$) 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.” 2. 2.

Continuous Stellar Mass-Loss (OB/AGB-star Winds): (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 $=\Delta t\,\dot{M}_{\rm wind}$. 3. 3.

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 re-radiate 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. 4. 4.

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_{\rm abs}\,\hat{n}/c$ (where $\hat{n}$ is the direction of ray propagation) to the gas. This automatically accounts for both direct UV/optical single-scattering, 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, multiple-scattering, 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 $\sim 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).

We 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 results indicate the effects on large scales are (unsurprisingly) small compared to our IMF-averaged approach.

For 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 photo-electric 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.

The 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 power-of-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).

The timestep is set by the minimum of various criteria. All particles obey limits $\Delta t<0.2\,(h_{i}/|{\bf a}_{i}|)^{1/2}$ (see Power et al. 2003; here $h_{i}$ is the minimum of the inter-particle separation or Plummer-equivalent force softening)555For reference, with our definitions, $\Delta t<0.2\,(h_{i}/|{\bf a}_{i}|)^{1/2}$ is equivalent to setting the parameter “ErrTolIntAccuracy” $\approx 0.01$ in GADGET-2. and $\Delta t<0.25/|\nabla\cdot{\bf v}_{i}|$, where ${\bf v}_{i}$ and ${\bf 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 $\Delta a<10^{-4}\,a$ (where $a$ is the scale-factor), but this is almost never important. Gas elements must also obey the Courant (CFL) condition: $\Delta t<0.4\,h_{i}/v_{{\rm sig},\,i}^{\rm max}$, where $v_{{\rm sig},\,i}^{\rm max}$ 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).666Some physics, such as cooling, photo-ionization, and recombination, are handled in a fully-implicit numerical scheme, which (in the limit where, say, the cooling time is much shorter than the timestep) iteratively solves for the equilibrium temperature balancing all heating and cooling physics over each timestep. This means they do not impose an additional explicit timestep criterion. In the above equations, note $h_{i}$ 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 $\sim 2$ if we replaced $h_{i}$ 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 Courant-like timestep criteria given in § 4.3.3.

For star particles, we additionally impose a restriction $\Delta t<{\rm MAX}(10^{4}\,{\rm yr},\ t_{\ast}/300)$, where $t_{\ast}$ 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.

Readers 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\sim 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 $\sim 5\,R_{\rm 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 $\sim 2\,R_{\rm 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_{\rm vir}$.

Table 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 histories 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” late-time 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 $\Lambda$CDM cosmology with $h\approx 0.70$, $\Omega_{\rm M}=1-\Omega_{\Lambda}\approx 0.27$, and $\Omega_{b}\approx 0.045$ (consistent with current constraints; see Planck Collaboration et al., 2014).777For the sake of comparison with other work, some ICs are matched to simulations which adopted very slightly different cosmological parameters. These differences are at the $\sim 1\%$ level and their effects are much smaller than standard halo-to-halo variation.

We 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.

For readers interested in more details, or reproducing our results, all initial conditions used here are publicly available.888For the MUSIC files necessary to generate all ICs here, see:

http://www.tapir.caltech.edu/~phopkins/publicICs

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).

With these optimizations, each high-resolution, production-quality simulations of a MW-mass galaxy with particle masses $\sim 7000\,M_{\sun}$ (a few $\times 10^{8}$ total particles) typically requires $\sim 10^{6}\,$cpu-hrs; for our smallest dwarfs (particle masses $\sim 250\,M_{\sun}$), their much lower star formation efficiencies and baryonic densities reduce this to $\sim 10^{4}\,$cpu-hrs. All simulations here were run on the XSEDE Stampede, Comet or NASA Pleiades machines.

Details 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.

More Accurate Hydrodynamic Solver: As described in § 2.1, FIRE-2 uses the newer, more accurate mesh-free finite-volume Godunov-type MFM method to solve the hydrodynamic equations. FIRE-1 used the older “pressure-energy” SPH (“P-SPH”) 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. 2. 2.

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 resolution-independent, it actually can influence galaxy morphologies more dramatically at the highest resolutions. 3. 3.

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, photo-electric 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 specifically, 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. 4. 4.

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 FIRE-2, 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. 5. 5.

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 $\sim 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. For 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. 6. 6.

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_{\rm halo}\sim 10^{9}-10^{12}$, similar to our FIRE-1 simulations. All the simulations here have been run to redshift $z=0$.

Fig. 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 ($m_{i,\,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 $\sim 10\,$kpc from the 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.

Fig. 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.

It 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\gtrsim 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.

Fig. 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.

In 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,999We 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_{\rm vir}$ and radius $R_{\rm 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\times 10^{4}$ DM, $<100$ baryonic, or $<10$ star particles), and subhalos (any halo within $<2\,R_{\rm 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_{\rm vir}$), then taking all stars within $3\times$ 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_{\rm vir}$ gives similar results, except one case with a $z\approx 0$ merger; but for massive galaxies the $0.1\,R_{\rm vir}$ cut includes stars that are clearly part of the halo. which is (as expected) always smaller than the total mass in $R_{\rm 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_{\ast}\gtrsim 10^{9}\,M_{\sun}$ (so it is extrapolated here), but this extrapolation has been shown to reproduce well the observed local group dwarf luminosity functions to $\sim 10^{4}\,M_{\sun}$ (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_{\ast}>10^{9}\,M_{\sun}$) with local group and field dwarf constraints. Garrison-Kimmel et al. (2016) perform a similar exercise, allowing the scatter below $M_{\rm halo}<10^{11}\,M_{\sun}$ to vary; we plot their best-fit median relations for a constant scatter below this mass $\approx 1\,$dex (but note that for any scatter in the range $0.5-2\,$dex, 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.2\,$dex above $M_{\rm halo}>10^{11.5}\,M_{\sun}$ (the value favored by Behroozi et al. 2012 and Moster et al. 2013), and varies linearly with halo mass as $\sigma_{\rm dex}=0.2-0.2\,\log_{10}(M_{\rm halo}/10^{11.5}\,M_{\sun})$ at lower masses (rising to $\approx 0.5\,$dex at $M_{\rm halo}\approx 10^{10}\,M_{\sun}$).101010At ultra-faint stellar masses $\lesssim 1000\,M_{\sun}$, it is likely that details of first star (Pop III) formation and primordial molecular (metal-free) cooling, not treated explicitly in FIRE, become important. This will be the subject of future study, but we consider these stages un-resolved given our mass resolution here, and simply initialize a metallicity “floor” of $10^{-4}\,Z_{\sun}$.

We 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 $\sim 10^{4}-10^{9}\,M_{\sun}$. 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 $\sim 10^{6}-10^{11}\,M_{\sun}$ to observations at $z\gtrsim 5$. However more quantitative comparison to the SMF of massive galaxies at lower redshifts ($z\sim 0-2$) will require larger statistical volumes.

In 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 FIRE-2 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.

As 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_{\rm halo}\sim 10^{11}\,M_{\sun}$, 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 function 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 FIRE-2 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 (Faucher-Giguere et al., 2015; Faucher-Giguère et al., 2016), although this will be studied in more detail in the future.

There 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 lower-mass 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 $\sim 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 $\sim 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_{\ast}-M_{\rm halo}$ relation at $M_{\rm halo}\leq 10^{11}\,M_{\sun}$. 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 $\sim 2-5$ systematic uncertainty in the calibration of observed galaxy metallicities. For MW-mass galaxies, the stellar masses are slightly higher in FIRE-2, 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.

Throughout 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 $m_{i}$ in the high-resolution Lagrangian region.111111The 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.”

To 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 $\sim 10^{6}$). Averaged over our runs at $z=0$, $\sim 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 $m_{i}$ in this paper.

Dark matter particles are always more massive by the universal baryon fraction, $m_{\rm dm}\approx 5\,m_{i}$; 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.

We define, for convenience, the baryonic particle mass in units of $1000\,M_{\sun}$, and note that this specifies both baryonic and dark matter particle masses:

[TABLE]

4.1.1 Requirements for “Resolved” Self-Gravitating Structures

The 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 $\sim 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:

[TABLE]

where $m_{i,\,1000}\equiv m_{i}/1000\,M_{\sun}$ is the baryonic particle mass, to which we reference all quantities, and $N_{10}=N_{\rm desired}/10$ reflects the “desired” number of particles ($N_{\rm desired}$).

We 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.

Of 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_{\ast}\sim 10^{4}\,M_{\sun}$ typically live in halos of mass $M_{\rm halo}\sim 10^{9}\,M_{\sun}$; using the $M_{\ast}-M_{\rm halo}$ relation at low masses from Moster et al. (2013), we estimate the minimum halo mass with a “believable” baryonic relic is $\sim 10^{9}\,M_{\sun}$ at $m_{i,\,1000}\sim 1$. This is actually quite well-resolved in dark matter, with $\gtrsim 10^{5}$ particles. This also means that $\gtrsim 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 feedback had a real effect (it had to prevent $>10^{5}$ particles from cooling and forming stars!).

Moreover, 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_{\rm halo}\gtrsim 10^{9}\,M_{\sun}$, independent of the number of stars (given a realistic $M_{\ast}-M_{\rm halo}$ relation).

We 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 $\sim 200$ elements. This is much more demanding than “existence” of halos or gaseous structures (which requires just a few elements), as expected.

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 $\sim 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 $\sim 2$ of the exact solution. This means that a single star particle generating SNe in a halo and blowing out $\sim 100-1000\,$ surviving gas particles (total baryonic mass $\gtrsim 10^{6}\,M_{\sun}\,m_{i,\,1000}$ or, at the Universal baryon fraction, halo mass $\gtrsim 10^{7}\,M_{\sun}\,m_{i,\,1000}$) is at least hydrodynamically well-behaved, if not converged.121212For reference, m10v, the lowest-baryonic mass system in Fig. 8, contains $\sim 10^{4}$ ($10^{7}$) gas elements inside $<1$ kpc ($R_{\rm vir}$) at $z=0$, at our fiducial resolution. 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.

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, 2012b, 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_{\rm GMC}^{\rm max}$ is set by the Toomre mass, $M_{\rm Toomre}\sim f_{\rm gas}^{2}\,M_{\rm gas}$ where $f_{\rm gas}$ is the gas fraction inside the relevant effective radius of the baryonic galaxy. So, if we require $M_{\rm cloud}^{\rm min}\lesssim M_{\rm Toomre}$, we obtain the desired resolution criterion $m_{i,\,1000}\lesssim m_{i,\,1000}^{\rm Toomre}$ where

[TABLE]

This defined the target resolution of the original FIRE-1 massive-galaxy simulations ($10^{12}\,M_{\sun}$ halos). 2. 2.

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_{\rm cool}\sim 3000\,M_{\sun}$; if a number $N_{\rm SNe}$ occur before the cooling time expires, the “cooling mass” simply scales $M_{\rm cool}\propto N_{\rm SNe}$. If SN energy is injected in a kernel-weighted fashion over $N_{\rm NGB}$ elements as we do here, then almost all of the energy goes into the nearest $N_{\rm NGB}^{1/3}$ neighbors; requiring they have a total mass $<M_{\rm cool}$ in turn requires:

[TABLE]

This set the target mass resolution for the original FIRE-1 dwarf galaxies ($<10^{10}\,M_{\sun}$ halos). As we show in Paper II, unless $m_{i,\,1000}\ll m_{i,\,1000}^{\rm SNe}$, it is necessary to properly account for the conversion of energy into momentum in the unresolved Sedov-Taylor (S-T) phases. Failure to do so will significantly under-estimate the effects of feedback.131313Note 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_{\rm SNe}\sim m_{0}/100\,M_{\sun}$ go off in an overlapping resolution element within a cooling time, then they can heat a mass $\sim 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. 3. 3.

Dwarf Galaxy “Burstiness”: In previous papers we have shown that for small dwarf galaxies, star formation is robustly “bursty”;141414For the sake of quantitative comparison, we define a specific measure of “burstiness” as the standard deviation in the quantity $\log{(\dot{M}_{\ast}(\Delta t_{1})/\dot{M}_{\ast}(\Delta t_{0}))}$ where $\dot{M}_{\ast}(\Delta t)$ is the SFR averaged over a time interval $\Delta t$, and we compare a short interval $\Delta t_{1}=10\,$Myr and longer interval $\Delta t_{0}=$ Gyr. 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 co-eval 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 $\sim 10$ particles (since a single lone gas particle cannot be at much higher density than its neighbors), then taking the IMF-average $N_{\rm SNe}\sim m_{i}/100\,M_{\sun}$ per particle, the entire baryonic mass in the galaxy can be heated to $T\gg 10^{6}\,$K (“super-bubble” temperatures) if that baryonic mass is below $\sim 600\,m_{i}$. So resolving “venting” of individual super-bubbles (much smaller than the entire galaxy) requires $\gtrsim 1000$ gas particles in the galaxy. Assuming further a gas fraction of $\sim 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:

[TABLE]

Equivalently, we can say that for a simulation with fixed mass resolution, once halos exceed a mass $M_{\rm halo}\gtrsim 2\times 10^{9}\,M_{\sun}\,m_{i,\,1000}$ the numerical “excess clustering/burstiness” is not important. 4. 4.

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\lesssim R_{\rm sonic}\sim 0.1\,{\rm pc}$ and mass scale $\sim 1\,M_{\sun}$). We discuss this below, but the turbulent Jeans mass of a super-sonically turbulent GMC with virial parameter $\alpha\sim 1$ is of order the GMC mass itself ($\sim 10^{7}\,M_{\sun}$), whereas the thermal Jeans mass is $\sim 0.1\,M_{\sun}$. 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 ($\sim{\rm 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; Guszejnov & Hopkins, 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 typical values for mass, spatial, and time resolution for m12i runs at different resolution.

Fig. 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.

For the dwarfs, the total stellar mass (and circular velocity profile) becomes robust to better than a factor $<2$, with just $\gtrsim 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 ($m_{i,\,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 galaxy. The metallicity and burstiness of the SFH appear to robust to a factor $\sim 2$ ($\sim 10\%$) when the number of stars reaches $>100$ ($>500$) at $m_{i,\,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 $m_{i,\,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 $\sim 5-10$ star particles per galaxy; but the metallicities of the satellites with $\lesssim 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 $\sim 10-20\%$ larger at the lowest resolutions here). In Paper II, we show the star formation histories of dwarf galaxy simulations reaching $\sim 30\,M_{\sun}$ resolution (which will be studied in more detail in Wheeler et al., in prep), and show that they agree well with the $\sim 250\,M_{\sun}$-resolution runs here.

For the massive, MW-mass system, even the lowest-resolution ($m_{i\,1000}=450$) run has $\gg 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 $\sim 10^{6}\,M_{\sun}$ – this means much lower “launch velocities” and temperatures, 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.

As 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 $m_{i,\,1000}=7$ run is $\sim 1.4$ times larger than the $m_{i,\,1000}=450$ run, which translates to a factor $\sim 1.8$ lower mass inside $<5\,$kpc; this in turn lowers the peak in the rotation curve from $\sim 370\,{\rm km\,s^{-1}}$ to $\sim 270\,{\rm km\,s^{-1}}$. Thus the galaxy rotation curve is noticeably “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 $m_{i,\,1000}=450$ there is still a prominent, extended thin disk. This may owe to the fact that m12m and m12f have a somewhat larger mass, but likely also owes to the specific fact that m12i has a series of mergers around $z\sim 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.

If the dependence on resolution in our MW-mass runs owes to a 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 $m_{i,\,1000}=56$, re-started at $z\approx 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 $\sim 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).

This 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.

To 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_{\rm relax}\sim 0.6\,t_{0}$ (where $t_{0}\equiv t_{\rm circ}(R_{200})=2\pi\,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/\sqrt{8})\,(N[<r]/\ln{N[<r]})\,(\rho_{\rm crit}/\bar{\rho}[<r])^{1/2}<0.6$, satisfied for the densities in Fig. 11 when $N\gtrsim 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 $\sim 1.3$ ($0.1\,$dex) of the converged (NFW) solution we find excellent agreement independent of resolution, down to radii containing just $N\gtrsim 200$ particles.

Fig. 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 $\lesssim 0.15\,$dex outside the radii enclosing $\sim 200$ particles. In m10v, the $z=0$ galaxy is sufficiently low-mass that baryons have a negligible 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_{\ast}$ and correspondingly lower central DM density). In all cases, our highest-resolution baryonic simulations reach an approximate “convergence radii” of $\sim 50-100$ pc in the dark matter. We of course achieve much higher effective spatial resolution in the baryons dense enough to form stars.

Note 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_{\rm relax}\gtrsim{\rm couple}\times t_{\rm 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 $\sim 2$ particles.

Fig. 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\sim 0$, but study the mass function of dense, cold gas structures within the ISM.151515For 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 $\leq 8000\,$K and densities $n>5\,{\rm cm^{-3}}$, within $0.1$ virial radii of the center of the main galaxy, at $\sim 30$ uniformly-spaced snapshots in time between $z=0.06-0$, and plot the time-averaged mass function of gas structures, $M\,dN/d\log_{10}{M}$ for each simulation. 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 resolution-dependence of the simulations. In each case, the shape of the mass function is, as expected, a power-law with slope $dN/dM\propto 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_{\rm 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 $\geq 3$ particles. We also compare the internal properties of the clouds, specifically the linewidth-size relation (the one-dimensional 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_{\rm cloud}$ is equivalent to their $\sigma_{r}$). The generic power-law scaling here is also predicted in turbulent fragmentation models (references above) and similar to observations.

Remarkably, 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 $\sim 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 ($m_{i\,1000}=470$), where the most massive clouds appear slightly more massive: this is because the clouds with $\sim 10^{6}\,M_{\sun}$ (not quite the peak of the MF, but just below it) contribute a significant total mass but cannot be resolved (this is $\sim 2$ particles), so that mass is “shifted” into more massive structures numerically. But even in this marginal case, the gas-mass weighted mean $M_{\rm GMC}$ is only over-estimated by a factor $\sim 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_{\sun}$, similar to what is observed in the Milky Way (Williams & McKee, 1997).

Note 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.

This further demonstrates, as we argued above, that resolving the thermal Jeans mass in cold gas is not important to cloud-or-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 ($\ll 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 $\sim 10^{3}-10^{5}$), yet they have no appreciable systematic effect on any results we measure.

Finally, 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 $\sim 10-30$ DM particles, the predictions are independent of resolution, consistent with other studies.161616Compare, for example, Fig. 14 here to Fig. 9 in Springel et al. (2008) (our runs with $m_{i,\,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 high-resolution 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). 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.

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 $h_{i}$ – it can, in principle, become arbitrarily small. In our Lagrangian method, the mean inter-particle spacing is directly related to the gas density $n_{\rm H}$, via

[TABLE]

In Hopkins (2015); Hopkins & Raives (2016) we show that the “effective” resolution for sound waves is identical to second-order 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 higher-order reconstruction, sound waves can be accurately re-constructed to scales about $\sim 1/2$ the inter-particle separation (Hopkins 2015 Fig. 2, Hopkins & Raives 2016 Fig. 1). Contact discontinuities can be captured and conserved across $\sim 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 $\sim 3$ elements, close to ATHENA (Hopkins 2015 Figs. 3, 10-15, 29-30). Intrinsically multi-dimensional 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 $\sim 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\gtrsim 1-2\,h_{i}$. In multi-dimensional super-sonic turbulence, numerical dissipation and noise truncates the inertial range at wavelengths approximately $\sim 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 $\epsilon$. 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,

[TABLE]

For clarity, we use the same definition of force softening and spatial resolution, corresponding to the inter-particle separation $h_{i}=\Delta x$ (defined such that $\rho_{i}=m_{i}/h_{i}^{3}$). For our default kernel, this means that the commonly quoted “Plummer equivalent” softening is $\epsilon_{\rm plummer}\approx(2/3)\,\epsilon=(2/3)\,h_{i}$. For the reasons in § 2.2 and § 4.2.2 below, we adopt constant force-softening parameters for DM ($\epsilon_{\rm DM}$) and stars ($\epsilon_{\ast}$), 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 $\epsilon$ will exert identical gravitational forces at all radii).

We will explore variations below, but first summarize our best estimate of “optimal” parameters here:

• •

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).171717Strictly speaking, it is only possible to solve the Poisson equation for a mass distribution accurate to the same order as the reconstruction order of the hydrodynamic method. But up to truncation errors, this is exact in matching the gas distribution. (2) Many authors have shown that this provides an optimal numerical 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 ambiguity about co-moving or physical softenings in cosmological integrations.181818Adaptive 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.

If constant gravitational softenings for gas must be used, they must be chosen sufficiently small to resolve the vertical scale-heights of the cold gas disk, and Toomre lengths of the most massive Toomre-mass objects discussed in our mass-resolution criteria ($\epsilon_{\rm gas}\ll 100\,$pc).

• •

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 ($\sim 100\times$ 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 $\nabla\cdot{\bf v}_{\rm DM}^{a}<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 particle distribution and their phase-space distribution; for discussion see Abel et al. 2012).

Therefore, 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 $\epsilon_{\rm DM}$ to the inter-particle spacing $h_{i}^{\rm DM}\sim 17\,{\rm pc}\,m_{i,\,1000}^{1/3}\,(\rho_{\rm DM}/10^{9}\,M_{\sun}\,{\rm kpc^{-3}})^{-1/3}$.

For 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_{\rm relax}\approx 0.06\,t_{\rm circ}(R_{200})$) by a factor of at least a couple, i.e. $\epsilon_{\rm DM}\lesssim 0.5\,r_{0.06}$. This is approximately ensured if $\epsilon_{\rm DM}<30\,{\rm pc}\,m_{i,\,1000}^{1/2}\,(M_{\rm vir}/10^{12}\,M_{\sun})^{-0.2}$ (we estimate this by taking $M_{\rm enc}(<r_{0.06})\approx 220\,m_{i}$, and assuming NFW profiles with concentration $c\approx 10\,(M_{\rm vir}/10^{12}\,M_{\sun})^{-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 $\epsilon_{\rm DM}>0.03\,{\rm pc}\,m_{i,\,1000}\,(M_{\rm vir}/10^{12}\,M_{\sun})^{-2/3}$ (obtained by comparing the hard-scatter $\Delta v$ to $V_{c}$; § 4.2.6). Note the large dynamic range in between (e.g. for our default runs, this gives $1\,{\rm pc}\lesssim\epsilon_{\rm DM}\lesssim 75\,{\rm pc}$). We will show we can vary $\epsilon_{\rm DM}$ by multiple orders of magnitude without changing out results.191919Note that it becomes impossible to satisfy both criteria for particle masses $m_{i}$ corresponding to $\lesssim 40$ DM particles in the halo; of course this is because the halos then cannot be internally spatially resolved. To be conservative we choose softenings about $\sim 2-3$ times smaller than the upper limit above (giving $\epsilon_{\rm DM}\sim h_{i}^{\rm DM}$ 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).202020Another common criteria for DM softening is that the maximal 2-body acceleration ($a_{i}\sim G\,m_{i}/\epsilon_{\rm DM}^{2}$) not exceed the bulk acceleration at some radius ($a_{\rm bulk}\sim G\,M_{\rm enc}(<r)/r^{2}$). If we assume an NFW profile, and require $a_{i}<\langle a_{\rm bulk}\rangle$ where $\langle a_{\rm bulk}\rangle$ is the mass-averaged mean acceleration of DM within the halo, we obtain $\epsilon_{\rm DM}\gtrsim 6\,{\rm pc}\,(c/10)^{-0.45}\,m_{i,\,1000}^{1/2}\,(M_{\rm vir}/10^{12}\,M_{\sun})^{-1/2}$ (we approximate the exact scaling with the halo concentration $c$ by a power-law, good to $\sim 10\%$ over the entire range of interest). Our fixed-$\epsilon_{\rm 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_{\rm vir}$, using order-of-magnitude smaller $\epsilon_{\rm DM}$. This is because this acceleration criterion is not meaningful when $\epsilon_{\rm DM}$ is below the inter-particle separation, because only very rarely will particles approach within separations $\ll\epsilon_{\rm DM}$, and when they do, the net velocity imparted by the encounter will be given by the scaling in § 4.2.6: $(\Delta v/\sigma_{v,\,{\rm DM}})\sim 10^{-4}\,m_{i,1000}\,(20\,{\rm pc}/\epsilon_{\rm DM})\,(\sigma_{v,\,{\rm DM}}/100\,{\rm km\,s^{-1}})^{-2}\ll 1$. Equivalently, the timescale required for $N$-body heating at $\sim R_{\rm vir}$ to perturb a particle orbit by $\sim 10\%$ is much greater than $t_{\rm Hubble}$ for all $\epsilon_{\rm DM}\gtrsim 1\,{\rm pc}$.

• •

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 $\epsilon_{\ast}$ 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” (self-gravitating, 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 $\epsilon_{\ast}$ fixed in time but variable particle-to-particle, to equal the softening of the gas from which the particle formed; while this does not change our results, it leads to wildly variable $\epsilon_{\ast}$ 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 $\epsilon_{\ast}$ 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 $\epsilon_{\ast}$ similar to the gas softening at the mean density of star formation, $\epsilon_{\ast}\sim 3.4\,{\rm pc}\,m_{i,\,1000}^{1/3}\,(\langle n\rangle_{\rm SF}/10^{3}\,{\rm 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 “effective 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 $h_{i}$) at:

Maximum Density/Minimum Softening: This is simply the minimum $\epsilon_{i}^{\rm gas}\equiv h_{i}^{\rm gas}$ 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_{\rm max}\sim 5\times 10^{6}\,{\rm cm^{-3}}$, which are reached by the simulation only briefly in the galactic nucleus in an intense, high-redshift starburst. 2. 2.

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 $h_{i}$. We show the mean density, $\langle n_{\rm SF}\rangle$, weighted by total star formation, integrated over the entire galaxy history to $z=0$, and corresponding $h_{i}^{\langle SF\rangle}$. We emphasize that because star formation is only allowed in self-gravitating gas which is also self-shielding and molecular, the mean $\langle n_{\rm SF}\rangle$ is always significantly larger than the minimum density at which star formation is allowed ($n_{\rm SF,\,min}$), even when the latter is set to $n_{\rm SF,\,min}=1000\,{\rm cm^{-3}}$. 3. 3.

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, 2013c, 2013b; Hopkins & Christiansen, 2013; Guszejnov & Hopkins, 2015; Guszejnov & Hopkins, 2016). To get some insight, however, for the warm ($\gtrsim 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 $\sim 10$ elements at these temperatures, and the corresponding resolved Jeans length (with $\sim 3\,h_{i}$) and gas density.212121Define the Jeans length $\lambda^{J}$ and mass $m_{J}=(4\pi/3)\,\rho\,(\lambda^{J}/2)^{3}$. Following § 4.1.1, a “resolved” structure has some number $N\sim 10\,N_{10}$ elements, so mass $\sim N\,m_{i}$. By definition of the inter-particle spacing, we also have $\rho=m_{i}/h_{i}^{3}$. This gives $\lambda_{J}\approx 3\,h_{i}\,N_{10}^{1/3}$, always – in other words, our adaptive softening ensures the Jeans length $\lambda_{J}$ is always spatially-resolved provided the Jeans mass $m_{J}$ is mass-resolved. Now use this and the definition of Jeans length, $\lambda^{J}=c_{s}/\sqrt{G\,\rho}$, to solve for the Jeans radius $\lambda^{J}/2\approx 9\,G\,m_{i}\,N_{10}/c_{s}^{2}$; for $T=10^{4}\,K$ (and assuming fully-ionized gas), this becomes $\approx 0.3\,{\rm pc}\,m_{i,\,1000}\,N_{10}$. Taking a minimum resolved size $N_{10}=1$, this gives the minimum resolved $\lambda_{J}$; the corresponding $h_{i}$ gives, in turn, the corresponding maximum density $n_{\rm max}^{J}$. 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 $\sim 100-200\,{\rm pc}$ and densities $\sim 10\,{\rm cm^{-3}}$), which form out of the warm gas. Note that if the gas is hotter, it becomes easier to resolve Jeans-scale structures. 4. 4.

Turbulent Jeans Mass (Cold Gas): Once the gas cools below temperatures $\sim 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 $\lambda_{J}=c_{s}/\sqrt{G\,\rho}$ by $\lambda_{\rm turb}\approx\langle v_{\rm turb}^{2}(\lambda_{\rm turb})\rangle^{1/2}/\sqrt{G\,\rho}$, where $v_{\rm turb}(\lambda_{\rm turb})$ is the rms turbulent velocity measured within the same region. A rigorous 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_{\rm turb}\propto\lambda^{1/2}$), and that the “parent” clouds have virial parameter of unity and (according to Larson’s Laws) a universal surface density $\Sigma_{\rm cloud}=\Sigma_{300}\,300\,M_{\sun}\,{\rm pc^{-2}}$, and take a minimum resolved-object mass $\sim N\,m_{i}$, then following Hopkins (2013b) gives a minimum-resolved length $\lambda_{\rm turb}\approx 4.0\,{\rm pc}\,(m_{i,\,1000}\,N_{10}/\Sigma_{300})^{1/2}$ and maximum-resolved density $n_{\rm max}^{\rm turb}\approx 1.2\times 10^{4}\,\Sigma_{300}\,(m_{i,\,1000}\,N_{10}/\Sigma_{300})^{-1/2}$. 5. 5.

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_{i}^{\rm DM,\,min}$, as well as rms value within twice the effective radius of the baryonic galaxy at $z=0$, $h_{i}^{\rm DM,\,core}$. 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 $\sim 200$ DM particles. We therefore quote this as well, for the DM.

Qualitatively, 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 ($\sim 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).222222For adaptive DM softening, we determine $\epsilon_{\rm DM}$ using the same methods as for gas, but only counting other DM as “neighbors.” In other words, we set $\epsilon_{i}^{\rm DM}=h_{i}^{\rm DM}$, where $h_{i}^{\rm DM}$ is determined from the dark matter particle neighbor distribution in the identical manner to $h_{i}^{\rm gas}$, ensuring $h_{i}^{\rm DM}$ is the mean (kernel-averaged) inter-particle spacing of dark matter particles within the kernel. This gives $\epsilon_{i}^{\rm DM}\equiv h_{i}^{\rm DM}=17\,{\rm pc}\ m_{i,\,1000}^{1/3}\,({\rho_{\rm DM}}/{10^{9}\,M_{\sun}\,{\rm kpc^{-3}}})^{-1/3}$. For adaptive softening for stars, we follow the same exercise but include all baryonic (gas+star) particles as “neighbors.”

To 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 $10\,$Myr 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\times$ larger softening $(\epsilon_{\rm DM},\,\epsilon_{\ast})=(200,\,20)\,{\rm pc}$; the “burstiness” is still obviously present but is quantitatively suppressed, and the stellar masses of both increase by a factor $\sim 1.4$. It makes sense that such large softenings suppress 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\%$).

Fig. 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\sim 200$ particles, these choices (with fixed $\epsilon_{\rm DM}$ varied from $\sim 10-1000$ pc) have almost no effect on the mass profile. Of course, if $\epsilon_{\rm DM}$ is too large, it will eventually suppress any smaller-scale structure (e.g. artificially flattening the DM cusp) – this occurs when $\epsilon_{\rm 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 $\epsilon_{\rm DM}<c^{-2}\,[\ln{(1+c)}-c/(1+c)]\,R_{200}\lesssim 0.01\,R_{200}$; this is easily satisfied by any of our simulations that satisfy $\epsilon_{\rm DM}\lesssim 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.

Fig. 17 does raise one important caveat: if $\epsilon_{\rm DM}\ll h_{i}$ (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_{\rm DM}=2.8\times 10^{6}\,M_{\sun}$ and $\epsilon_{\rm DM}=10\,$pc; 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 $\epsilon_{\rm DM}\gtrsim 0.02\,{\rm pc}\,m_{i,\,1000}$, easily satisfied in all our production simulations.

Fig. 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 $\epsilon_{\rm DM}\gtrsim r_{0.06}$. Although the effect is small, this is expected: if the profiles are converged down to $\sim r_{0.06}$ with appropriate softening, then over-softening by setting $\epsilon_{\rm DM}\gtrsim r_{0.06}$ will necessarily smear the profile out at these radii.

Fig. 19 shows the subhalo mass function, $V_{\rm 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 good 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.

In Figs. 15-16, we also considered the effects of changing the gas force softening, replacing our default self-consistent adaptive softening ($\epsilon_{\rm gas}=h_{i}$) with a fixed, constant physical softening ($\epsilon_{\rm gas}$). In both our dwarf and MW-mass simulations, we see, reassuringly, that for sufficiently small fixed $\epsilon_{\rm gas}$, the differences are essentially negligible (entirely consistent with stochastic fluctuations).

However, recall that our star formation model is based on identifying self-gravitating gas above some density threshold. At fixed mass resolution, a fixed $\epsilon_{\rm gas}$ sets a minimum inter-particle separation, hence maximum gas density, at which self-gravity will be correctly calculated: this density is $n^{\rm max}=m_{p}^{-1}\,\,m_{i}\,(\epsilon_{\rm gas}^{\rm min})^{-3}\approx 1000\,{\rm cm^{-3}}\,m_{i,\,1000}\,(\epsilon_{\rm gas}^{\rm min}/3.5\,{\rm pc})^{-3}$. So not surprisingly, when we make $\epsilon_{\rm gas}$ large enough that $n^{\rm max}\lesssim n_{\rm crit}$ (the minimum density for star formation), we can artificially suppress the SFR. Obviously, one should not therefore, with fixed $\epsilon_{\rm gas}$, choose $n_{\rm crit}\gtrsim n^{\rm max}$. If we lower $n_{\rm crit}$ sufficiently so that $n_{\rm crit}\ll n^{\rm max}$, then we recover nearly-identical behavior to our default simulations with adaptive $\epsilon_{\rm gas}$ and larger $n_{\rm crit}$.

This is important: some simulations allow $n_{\rm crit}\ll n_{\rm max}$ (e.g. Guedes et al. 2011; Shen et al. 2014 where $\epsilon_{\rm gas}\sim 100-200\,{\rm pc}$, $n^{\rm max}\sim 0.1\,{\rm cm^{-3}}$, and $n_{\rm crit}\sim 5-100\,{\rm cm^{-3}}$) – external forces (e.g. shocks) can still produce $n>n_{\rm crit}$ (if $h_{i}$ is allowed to be smaller than $\epsilon_{\rm gas}$), but this means star formation has nothing to do with self-gravity locally. We therefore consider two similar experiments (${\bf m10q}$ with $\epsilon_{\rm gas}=20\,$pc, $n^{\rm max}=10\,{\rm cm^{-3}}$, and m12i with $\epsilon_{\rm gas}=140\,$pc, $n^{\rm max}=1\,{\rm cm^{-3}}$; both with $n_{\rm crit}=100\,{\rm cm^{-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;232323Gas which can cool efficiently to $T\sim 10\,$K and has turbulent $\delta v\sim 10\,{\rm km\,s^{-1}}$ should produce an approximately lognormal density distribution with $\gtrsim 10\%$ of the mass exceeding $\sim 100\times$ the mean density (see e.g. Federrath et al., 2008). the overdensities cannot “detach” from the turbulent flow and collapse to still higher densities (like real GMCs) if $\epsilon_{\rm gas}$ is too large (which means the properties of the cold, dense gas will be incorrect), but they can still reach $n_{\rm crit}$, at which point the galaxy-averaged SFR is self-regulated by a balance of inflow and feedback. In the m10q ($\epsilon_{\rm gas}=20\,$pc) run, the too-large $\epsilon_{\rm gas}$ still suppresses the SFR and stellar mass by a factor $\sim 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 $\epsilon_{\rm 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 the 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_{\rm encounter}$. Since the duration of the encounter is short, the impulsive velocity change is $\Delta v=|\Delta{\bf v}|\sim G\,m_{i}/(h_{i}\,v_{\rm encounter})$ (a proper integration over the exact kernel shape used in the simulations gives a maximal deflection $\approx 1.3\,G\,m_{i}/(h_{i}\,v_{\rm encounter}$), which occurs for encounters with impact parameter $\approx h_{i}$). For gas, this is $\Delta v\sim 0.01\,{\rm km\,s^{-1}}\,m_{i,\,1000}^{2/3}\,(n/{\rm cm^{-3}})^{1/3}\,(v_{\rm encounter}/10\,{\rm km\,s^{-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, $\sim 100\,{\rm km\,s^{-1}}$. So $\Delta v\sim 1.3\,G\,m_{\rm DM}/(h_{\rm DM}\,v_{\rm encounter})\sim 0.017\,{\rm km\,s^{-1}}\,m_{i,\,1000}\,(h_{\rm DM}/20\,{\rm pc})^{-1}\,(v_{\rm encounter}/100\,{\rm km\,s^{-1}})^{-1}$, only slightly larger even when we adopt the smallest DM softening seen in our simulations.242424In fact, the $\Delta 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 $|\Delta\rho|/|\rho|\sim 0.005$. Thus the actual hard-scattering amplitudes are suppressed by a factor of $\sim 100$ from their already-small values.

We can translate this to an “N-body heating rate,” assuming each encounter adds to the velocity incoherently, with particle encounter rate $\sim v_{\rm encounter}/h_{i}$. Accounting for the number of baryons per gas particle ($\sim m_{i}/m_{p}$), this gives a heating rate (rate-of-change of kinetic energy) of:

[TABLE]

These should be compared to the physical heating/cooling rates of the gas, $\Lambda\sim 10^{-23}-10^{-22}\,{\rm erg\,cm^{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.252525We 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^{\rm gas-DM}_{\rm 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^{\rm gas-gas}_{\rm heat}$ is unmeasurably small.

4.3 Time Resolution

The time resolution in the simulations is set by the timestep $\Delta 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 $\Delta t<0.2\,(h_{i}/|{\bf 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.

In either case, the minimum timesteps reached in our simulations are $\approx 100\,$yr (see Table 3). This is reached regularly by some elements from redshifts $z\sim 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

We have 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 Power et al. (2003), we find that a coefficient smaller than $\alpha=0.2$ in the timestep limiter $\Delta t<\alpha\,(h_{i}/|{\bf a}_{i}|)^{1/2}$ produces no appreciable gain in accuracy, but much larger values ($\alpha\gtrsim 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 ${\bf a}_{i}$ includes all accelerations, i.e. hydrodynamic and gravitational), a Courant-type condition is still necessary, with coefficients $\alpha>0.8$ in $\Delta t<\alpha\,h_{i}/v_{{\rm sig},\,i}^{\rm max}$ giving rise to numerical instability, while values $\alpha<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).

4.3.2 Stellar-Evolution Timestep Limits

Recall, we do not allow the timestep for star particles to exceed

[TABLE]

where $t_{\ast}$ is the age of the star, to ensure stellar evolution is time-resolved 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 $\approx 3$ (which produced no measureable effect), or most radically using $\Delta t_{\ast}<20\,{\rm Myr}\,m_{i,\,1000}^{-1}$ – for $m_{i,\,1000}\gtrsim 1$, this is equivalent to allowing as many as $\sim 10\,$SNe 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 ($\sim 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 $\sim$Myr. However, if we removed the limiter entirely ($\Delta t_{\ast}\rightarrow\infty$), 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).

However, 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 $\epsilon$ over the course of their orbits.

Fig. 20 shows the results of using adaptive softening in DM-only simulations (m10v) of varying resolution, without imposing any additional timestep limiter beyond what is described above for collisionless particles with fixed softening $\epsilon$. Clearly, the central sub-kpc regions of the DM profile exhibit a spurious density “cusp” well in excess of the converged (much higher-resolution) results from constant-$\epsilon_{\rm 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 fixed-$\epsilon$ 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.

With 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 ($|{\bf 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.

For 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:

[TABLE]

Here $H_{a}\approx 2\,\epsilon_{a}$ (with $\epsilon_{a}=h_{a}$) 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 $\epsilon$ of particle $a$. By analogy to the Courant condition, we define $v_{\rm sig,\,AGS}^{a}$ as the maximum approach/recession velocity of any neighbor within this interaction kernel – this requirement is simply that two particles cannot “cross” more than $\sim 1/2$ their relative softening lengths in a single timestep. The particle-divergence $\langle\tilde{\nabla}\cdot{\bf v}\rangle_{a}$ is defined exactly the same as the traditional SPH velocity divergence,262626With adaptive softening for DM, in 3D, we define $\epsilon_{a}=h_{a}=(3\,N_{\rm eff}/4\pi)^{1/3}\,H_{a}$ such that, in 3D, $h_{a}^{3}\,\bar{n}_{a}=1$ and $4\pi/3\,H_{a}^{3}\,\bar{n}_{a}=N_{\rm eff}$ (with $N_{\rm eff}=32$ for our standard cubic spline kernel), where $\bar{n}_{a}\equiv\sum W({\bf x}_{ba},\,H_{a})$ is a kernel-averaged particle neighbor number density. With some straightforward algebra, this gives an exact discrete equation for the Lagrangian derivative of $h_{a}$, $Dh_{a}/Dt=-(h_{a}/3\Omega_{a}\bar{n}_{a})\,\nabla\bar{n}_{a}\cdot\partial({\bf x}_{b}-{\bf x}_{a})/\partial t=-(h_{a}/3)\,\langle\tilde{\nabla}\cdot{\bf v}\rangle_{a}$, where $\Omega_{a}\equiv 1+(h_{a}/3\bar{n}_{a})\,\partial\bar{n}_{a}/\partial h_{a}$ (for derivations, see Price, 2012; Hopkins, 2013a). So if we do not want the softening $h_{a}$ to change by more than a factor $|(\Delta h_{a})/h_{a}|<\alpha$ in one timestep, we require a timestep $\Delta t<3\alpha/|\langle\tilde{\nabla}\cdot{\bf v}\rangle_{a}|$. in such a manner that this requirement prevents $\epsilon_{a}$ from changing by more than $\sim 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_{{\rm 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.

Fig. 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 $\sim 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).

5 Effects of the Hydrodynamic Method

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. Fortunately, GIZMO is an inherently multi-method code, so we can compare simulations with otherwise identical physics and numerics, replacing only the hydrodynamic solver.

Figs. 21, 22, 23, & 24 compare GIZMO simulations using both our (1) default (non-SPH) MFM method,272727We 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 problems (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. (2) various “flavors” of SPH.

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:

• •

FIRE-1 SPH: This is the SPH formulation used in FIRE-1, 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 $\sim 64$ effective neighbors instead of $\sim 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-least-squares 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 fluid-mixing and shock capturing, while reducing noise and artificial numerical diffusivity away from shocks (see references above).

Readers interested in more details of the SPH method here should consult the public version of GIZMO: this is the default SPH implementation.

• •

FIRE-1 SPH + Stronger Mixing: A common 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.

With 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:

[TABLE]

where ${\bf x}_{ba}\equiv{\bf x}_{b}-{\bf x}_{a}$ is the particle separation, $m$ mass, $\bar{\rho}$ the SPH-estimated density, $W$ the SPH smoothing kernel as a function of ${\bf x}_{ba}$ and smoothing length $h_{a}$, $c_{s}$ the sound speed, ${\bf v}$ the velocity, $P$ the pressure, and $\alpha_{a}^{V}$ is the artificial viscosity coefficient defined in Hopkins 2015, Appendix F2, Eq. F16-F17. This $\alpha^{V}$ varies between $0.05-2$ depending on the velocity divergence $\nabla\cdot{\bf v}$ and its time-derivative, reaching large values when particles approach increasingly rapidly ($\nabla\cdot{\bf v}<0$ and $d[\nabla\cdot{\bf v}]/dt<0$) and decaying rapidly otherwise.

The 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 $\approx 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” $\alpha_{ab}$ which suppresses diffusion in supersonically receding flows and sharp phase discontinuities.

• •

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 $\sim 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_{\ast}\lesssim 10^{10}\,M_{\sun}$ at redshifts $z\gtrsim 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.

This is not surprising: such galaxies are in the “cold accretion” regime where there is little or no “hot halo” of virial shock-heated gas 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).

Figs. 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 MFM and SPH at low redshift, the SFRs are initially identical, and only slowly drift apart (here by just $\sim 30\%$ over $\sim 1\,$Gyr). Moreover Fig. 24 shows that the visual galaxy morphologies in SPH and MFM are very similar.

These results are consistent with a series of studies of non-cosmological 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_{\ast}\gtrsim 10^{10}\,M_{\sun}$, corresponding to halo masses $\gtrsim 10^{11.5}\,M_{\sun}$. In the D-SPH runs (which maximize the difference), the final stellar mass is suppressed relative to the MFM run by a factor of $\sim 1.5$ (low resolution) or $\sim 3$ (intermediate resolution), giving a correspondingly smaller peak circular velocity.

Similarly, Davé et al. (2016) compare low-resolution cosmological large-volume simulations using a simplified, non-FIRE sub-grid ISM and feedback model and MFM and “P-SPH” (FIRE-1 SPH) methods in GIZMO. Over the mass range simulated ($M_{\ast}\sim 10^{9}-10^{12}\,M_{\sun}$), the stellar mass functions agree very well at $z>1$, but then begin to differ, with SPH tending towards smaller masses by $\sim 0.2\,$dex 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.

Not 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és-Alcá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).

Fig. 25 shows a visual comparison of the gas morphology and phase structure in the CGM at redshifts $z\sim 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.

Fig. 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 shock-capturing, 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 $\mathcal{M}\sim 1$, and requires $\sim 8-10$ inter-particle spacings ($\sim 30-50\,$kpc at the virial radius for the low resolution in Fig. 25) to fully capture shock jumps.

The 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 ($\sim 10^{5}\,$K) and cool ($\sim 10^{4}\,$K) CGM gas, we see clearly that in D-SPH there is a large amount of cold ($T\ll 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\gtrsim 10^{6}\,$K with $n_{H}\gtrsim 10^{-4}\,{\rm cm^{-3}}$), owing to the difficulties of shock-capturing.

Together, 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).

We stress, 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.

This 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).

This also suggests that the large central $V_{c}$ “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.

For 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 self-consistently, 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_{\rm artificial}\sim G\,(N\,h\,\rho)^{2}/\gamma$ is added to the equations of motion, where $h$ is the resolution (local gravitational softening), $\rho$ 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_{\rm artificial}\propto\rho^{4/3}$).

A common misconception is that this “prevents artificial fragmentation,” but this is not accurate (see Teyssier, 2015); in cases where $P_{\rm 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 $\sim N\,h$, 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 precisely 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).282828We 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, $\sim 0.1\,M_{\sun}$, in the cold ISM), one needed to maintain $\gtrsim 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.

Furthermore, “artificial pressure” prescriptions require careful treatment, as they: (a) can violate energy conservation (providing an infinite source of “$PdV$ work”); (b) are designed for thermally-supported 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. Specifically 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\sim 1$), and (3) fully-adaptive gravitational softenings are used for gas ($\epsilon_{\rm gas}=h_{i}$), 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.

That 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\sim 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\lesssim 30\,{\rm pc}\,m_{i,\,1000}^{1/3}\,(n/10\,{\rm cm^{-3}})^{-1/3}$; so for the resolution and typical gas densities inside the central $\sim 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 low-temperature 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\gtrsim 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 Do the Details of Cooling & Chemistry Matter?

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 $\sim 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.).

Figs. 27-28 show the results of turning on or off different portions of the cooling physics in our simulations.292929For the m12i example in Fig. 28, we show results at very low resolution, which maximizes the contrast between different cooling models. At higher resolution ($m_{i,\,1000}=56$) the results are similar but with slightly weaker differences. Specifically, we compare 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 $\Lambda(T,\,Z,\,I_{\nu},\,n)\rightarrow\Lambda(T)$ for fully ionized, non-self-shielded gas (so there is essentially no cooling to $T\ll 10^{4}\,$K), with a single metallicity $Z=Z_{\sun}$ 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 cooling physics (much more extreme than typical model differences in cooling) have relatively modest impacts on galaxy formation.

For dwarf galaxies, with either $Z=0$ or $Z=Z_{\sun}$ (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_{\sun}$ 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 re-starts near $z=0$ with fixed initial conditions (Fig. 27), if we adopt approximately the “correct” metallicity in and around the galaxy, $Z=Z_{\sun}$. However for MW-mass systems, the primordial cooling-only ($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.

Fig. 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\ll 10^{4}\,$K), by construction. For the toy model with $Z=Z_{\sun}$, 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).

This weak dependence of galaxy properties on the more subtle details of cooling physics at $T\lesssim 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.

The reason for this is simple: for almost all gas (by mass) in galaxies, the cooling time is much shorter than the dynamical time ($t_{\rm cool}\ll t_{\rm 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_{\rm cool}\lesssim t_{\rm dyn}$ for gas with $T\ll 10^{7}\,{\rm K}\,(n/0.01\,{\rm cm^{-3}})$ (although for SNe blastwaves, the more relevant comparison is to the remnant expansion time, giving $T\lesssim 10^{6}\,$K), and (2) gas in this rapidly-cooling regime has equilibrium temperatures $\lesssim 10^{4}$ K or colder (the colder gas has no significant thermal pressure support, so assigning it $T\sim 10^{4}$ does not appreciably change its large-scale dynamics).

Basically, in gas with temperatures $\lesssim 10^{6}-10^{7}$ K, details of the cooling physics are not important. Above this threshold, $t_{\rm cool}\gtrsim t_{\rm 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_{\rm cool}\ll t_{\rm dyn}$ to $t_{\rm cool}\gg t_{\rm 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_{\rm vir}\ll 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\sim 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).

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).

Under some circumstances, it may be desirable to explicitly model this microphysical metal transport between resolution elements. In Appendix F.3 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 $\kappa_{\rm turb}$ set to the resolution-scale “eddy diffusivity” ($\sim\ell\,|\Delta v_{\rm turb}(\ell)|$; 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 GIZMO 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.

Fig. 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\lesssim 10^{4}\,$K, $n_{H}\gtrsim 0.1$ cm*-3*) forms a slightly tighter (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.

Su 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.

Fig. 32 further demonstrates this by considering variations to the standard metal-mixing term, at both dwarf and MW mass scales. 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 $\kappa_{\rm 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 F.3) with a constant $\kappa_{\rm turb}=10^{24}\,{\rm 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 $\kappa_{\rm turb}$ corresponds crudely to typical values of the turbulent diffusivity estimated by Eq. 96 in Appendix F.3 on the maximal turbulent scales (for $\ell\sim 100\,$pc and $|\Delta v_{\rm turb}(\ell)|\sim 10\,{\rm km\,s^{-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 the 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_{\sun}$ by $\sim 0.4\,$dex 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 ($m_{i,\,1000}=0.25$) and m12i ($m_{i,\,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.

Given 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).

To 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, 2013e) and Kim et al. (2013) for more extensive discussions (see also § 4.1.3).

Figs. 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_{\rm crit}$ with $n_{\rm crit}=1000\,{\rm cm^{-3}}$), to variations with $n_{\rm crit}=5-1000\,{\rm cm^{-3}}$; turning on/off the self-gravity (virial), molecular, and Jeans criteria; and arbitrarily multiplying the SFR per free-fall time $\dot{\rho}_{\ast}$ in the gas which meets these criteria by a constant factor $\epsilon_{\rm sf}=0.01-100$.

In 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 $\sim 10^{6}\,M_{\sun}$; isolated (non-cosmological), galaxy-merger, and fully-cosmological simulations; 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, 2013d, 2013e, 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 “multi-free-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ć et al., 2018).

This 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 $\sim 100\,$Myr). 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_{\rm crit}$ we set will eventually be exceeded, and even if $\epsilon_{\rm sf}\ll 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ć et al., 2018), giving a low time-averaged efficiency per free-fall 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.

Of 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 $\sim 3$ (“strong”) or $\sim 1/3$ (“weak”), relative to the “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 self-regulating at the level where feedback offsets gravitational collapse, so if feedback is $3\times$ stronger “per star,” then $3\times$ 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 competition between inflow and galactic outflows, and the high-speed 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 $\sim 2$ in the MW-mass system, going from $\sim 30\%$ to $\sim 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\sim 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 (late-time metallicities are lower for massive galaxies because there is no SNe enrichment, only RSG/AGB enrichment).303030Note the surprisingly lower $V_{c}$ for the no-SNe, MFM run owes to very early gas exhaustion leading to the galaxy growing only via “dry” gas-poor mergers at late times, which puff up the central dispersion. 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 SNe 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.

Turning 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).

Removing all radiative feedback – including the effects of the meta-galactic UV background – has a dramatic effect on dwarf stellar masses, even at halo masses $\gtrsim 10^{10}\,M_{\sun}$ with $V_{\rm max}\gtrsim 40\,{\rm km\,s^{-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 $\sim 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, but removing local radiative feedback produces a clear, significant increase in the central $V_{c}$ (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.

As 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\sim 0$ MW-mass galaxy as Fig. 10. This allows us to separate non-linear, long-timescale cosmological effects on the SFR (e.g. galactic wind generation and recycling), from the instantaneous self-regulation of star formation within a galaxy.

As 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 missing 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).

Why 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 simple 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 $\langle\dot{M}_{\ast}\rangle$ and superwind mass-loss rates $\langle\dot{M}_{\rm wind}\rangle\equiv\eta\,\langle\dot{M}_{\ast}\rangle$ (where $\eta$ is the wind mass-loading, which can depend on arbitrary galaxy properties) should equal the inflow rate $\dot{M}_{\rm in}$ (we assume gas mass in the disk is steady-state; see Davé et al. 2011a). Thus $\langle\dot{M}_{\ast}\rangle\sim(1+\eta)^{-1}\,\dot{M}_{\rm 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, at the point where the momentum flux per unit area from feedback, $\sim(p_{\ast}/m_{\ast})\,\dot{\Sigma}_{\ast}$ (where $p_{\ast}/m_{\ast}\sim 3000\,{\rm km\,s^{-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 $\sim G\,\Sigma_{\rm disk}\,\Sigma_{\rm gas}$ (an identical scaling is obtained assuming feedback must destroy star-forming clouds, and/or offset turbulent dissipation, in a $Q\sim 1$ disk; see Shetty & Ostriker, 2008; Kim et al., 2011; Faucher-Giguère et al., 2013; Hayward & Hopkins, 2017). Thus $\langle\dot{M}_{\ast}\rangle\sim(p_{\ast}/m_{\ast})^{-1}\,G\,\Sigma_{\rm disk}\,\langle M_{\rm gas}\rangle$. Of course, these two scalings must be consistent: if the “instantaneous” $\dot{M}_{\ast}$ is “too low” compared to the cosmological scaling above, then $\dot{M}_{\rm in}$ is not offset by star formation+outflow and gas “piles up,” raising $M_{\rm gas}$ until the SFR “catches up.” If the radiative feedback contributes significantly to $p_{\ast}/m_{\ast}$ (either directly or via pre-processing GMCs before SNe), then removing it lowers $p_{\ast}/m_{\ast}$ and in turn increases the SFR given fixed $\Sigma_{\rm disk}$ and $M_{\rm gas}$, as in Fig. 37. However, if $\eta$ is less sensitive to radiative feedback, the winds from higher $\dot{M}_{\ast}$ would be stronger, and so the galaxy will self-regulate at the same $\langle\dot{M}_{\ast}\rangle$. Thus, the instantaneous self-regulation (the Kennicutt-Schmidt relation) is not the same as self-regulation of cosmological growth (e.g. the stellar mass-halo mass relation).

10 Summary & Conclusions

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 $\sim 100\,M_{\sun}$ mass scales; see Table 3 for details).

In 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 multi-phase 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.

It 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 star-forming 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).

We 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.

We 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

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 galaxy 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 $\sim 100$ star particles to sample the self-enrichment history and a baryonic particle mass $\lesssim 10^{-6}\,M_{\rm 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). 2. 2.

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 $\sim(0.005-0.5)\,r_{0.06}$, where $r_{0.06}$ is approximately the radius containing $\sim 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. 3. 3.

We show that the radius $r_{0.06}$ is the radius to which DM profiles are converged to $\sim 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. 4. 4.

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 $\sim 1-10$ pc in MW-mass systems and $\sim 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. 5. 5.

If, for some reason, constant gravitational force softening for gas is desired in simulations with a resolved multi-phase 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 ($\epsilon\ll 100\,$pc), and (2) the minimum “threshold” density for star formation is resolved gravitationally, so that actual self-gravitating regions are what form stars ($\epsilon\lesssim 7\,{\rm pc}\,m_{i,\,1000}^{1/3}\,(n_{\rm SF,\,min}/100\,{\rm cm^{-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. 6. 6.

In addition to standard timestep criteria, guaranteeing proper time-resolution of feedback requires a limiter (Eq. 12) such that star particles cannot “skip” stages of stellar evolution or experience $\gg 1$ SNe per particle per timestep (§ 4.3). 7. 7.

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 Courant-like condition must be used to ensure stability (Fig. 20; Eq. 13).

10.3 Hydrodynamic Methods

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 moving-mesh 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.” 2. 2.

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). 3. 3.

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 over-estimate (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). 4. 4.

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 transport – 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

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\lesssim 10^{5}\,$K). Detailed variations in yields, which species are tracked, and numerical metal-mixing produce correspondingly weak effects. 2. 2.

However, in hot gas ($T\gg 10^{6}\,$K) in the CGM or SNe-heated 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. 3. 3.

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). 4. 4.

Provided the above criteria are met, SF is instead feedback-regulated. Uniformly increasing/decreasing the strength of feedback directly shifts both the integrated stellar mass, and position of galaxies on the KS law (Fig. 35). 5. 5.

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. 6. 6.

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). 7. 7.

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, even at mass scales as large as $M_{\ast}\gtrsim 10^{7}\,M_{\sun}$. We will explore these effects further in future work. 8. 8.

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:

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. 2. 2.

In addition to standard (optically thin, primordial) cooling, some accounting for self-shielding and low-temperature ($T\ll 10^{4}\,$K) cooling is included (to enable fragmentation), and high-temperature metal-line cooling is included to account for faster cooling in metal-enriched halo gas. 3. 3.

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. 4. 4.

Stellar mass loss is included, to provide a continuous additional fuel supply to the galaxy. 5. 5.

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. 6. 6.

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). 7. 7.

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.

As 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.

In this paper, we study only systems with halo masses $\lesssim 2\times 10^{12}\,M_{\sun}$. This is because it is widely believed that feedback from supermassive black holes (not included in our simulations here) 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).

Finally, 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.

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 core-collapse 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.

SNe Ia: These occur with rate-per-unit-stellar mass $R_{\rm Ia}=dN_{\rm Ia}/dt=0$ for $t_{\rm Myr}<37.53$ (where $t_{\rm Myr}$ is the age of the star particle in Myr), then $R_{\rm Ia}/({\rm SNe\,Myr^{-1}\,M_{\sun}^{-1}})=5.3\times 10^{-8}+1.6\times 10^{-5}\,\exp{\{-[(t_{\rm Myr}-50)/10]^{2}/2\}}$ for $t_{\rm Myr}\geq 37.53$. Since the rate is per-stellar-mass, the expectation value of the number of SNe for a star particle of mass $m_{b}$ and timestep $\Delta t_{b}$ is $=R_{\rm Ia}(t_{\rm Myr}^{b})\,m_{b}\,\Delta t_{b}$, and the trial for “success” (an explosion) is drawn from a binomial distribution. Each SN Ia has ejecta mass $M_{\rm ej}=1.4\,M_{\sun}$ and energy $E_{\rm ej}=(1/2)\,M_{\rm ej}\,v_{\rm ej}^{2}=10^{51}\,{\rm 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\times 10^{-6},$ $0.143,$ $0.0045,$ $0.0086,$ $0.156,$ $0.087,$ $0.012,$ $0.743)\,M_{\sun}$. 2. 2.

SNe II: The core-collapse SNe rate can be surprisingly well-fit by a simple piecewise-constant function, $R_{\rm II}=0$ for $t_{\rm Myr}<3.401$ or $t_{\rm Myr}>37.53$; $R_{\rm II}/({\rm SNe\,Myr^{-1}\,M_{\sun}^{-1}})=5.408\times 10^{-4}$ for $3.4<t_{\rm Myr}<10.37$; and $R_{\rm II}/({\rm SNe\,Myr^{-1}\,M_{\sun}^{-1}})=2.516\times 10^{-4}$ for $10.37<t_{\rm Myr}<37.53$. The IMF-averaged ejecta mass per explosion is $M_{\rm ej}=10.5\,M_{\sun}$, with ejecta energy $E_{\rm ej}=(1/2)\,M_{\rm ej}\,v_{\rm ej}^{2}=10^{51}\,{\rm erg}$. These are normalized so that the ejecta mass and energy exactly match the integrated totals from STARBURST99; the time-averaged energy injection rate ($\langle R_{\rm II}\,E_{\rm ej}\rangle$) is within $\sim 10\%$ of the tabulated STARBURST99 rate at all times. The IMF-averaged yields are (He, C, N, O, Ne, Mg, Si, S, Ca, Fe) $=(3.87,$ $0.133,$ $0.0479\,\tilde{N},$ $1.17,$ $0.30,$ $0.0987,$ $0.0933,$ $0.0397,$ $0.00458,$ $0.0741)\,M_{\sun}$, where $\tilde{N}={\rm MAX}(Z/Z_{\sun},\,1.65)$ accounts for the strongly progenitor-metallicity dependent N yields. The total metal yield $Z$ is given by summing the explicitly-followed species, with an additional $\sim 2\%$ added to account for un-tracked species ($Z=1.02\,\sum\,Z_{i,\,{\rm followed}}$); the remaining ejecta is H.313131For reference, the yields from Woosley & Weaver (1995), used in FIRE-1, are (He, C, N, O, Ne, Mg, Si, S, Ca, Fe) $=(4.03,$ $0.117,$ $0.0399\,\tilde{N},$ $1.06,$ $0.169,$ $0.0596,$ $0.0924,$ $0.0408,$ $0.00492,$ $0.0842)\,M_{\sun}$. In both FIRE-1 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 $|dM_{Z}({\rm species})/dZ_{\rm 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 $|dM_{Z}({\rm species})/dZ_{\rm 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). 3. 3.

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_{\ast}$ is $\dot{M}_{\rm w}=M_{\ast}\,f_{\rm w}\,{\rm Gyr^{-1}}$ with $f_{\rm w}=4.763\,(0.01+Z/Z_{\sun})$ for $t_{\rm Myr}<1$; $f_{\rm w}=4.763\,(0.01+Z/Z_{\sun})\,t_{\rm Myr}^{1.45+0.8\,\ln(Z/Z_{\sun})}$ for $1<t_{\rm Myr}<3.5$; $f_{\rm w}=29.4\,(t_{\rm Myr}/3.5)^{-3.25}+0.0042$ for $3.5<t_{\rm Myr}<100$; and $f_{\rm w}=0.42\,(t_{\rm Myr}/1000)^{-1.1}/(19.81-\ln{(t_{\rm Myr})})$ for $t_{\rm Myr}>100$. The total (IMF-averaged) kinetic luminosity of the mass-loss is given by $L_{\rm kinetic}=(1/2)\,\dot{M}_{\rm w}\,\langle v_{\rm w}^{2}\rangle=\psi\times 10^{12}\,{\rm erg\,g^{-1}}\,\dot{M}_{\rm w}$, with $\psi=(5.94\times 10^{4})/(1+(t_{\rm Myr}/2.5)^{1.4}+(t_{\rm Myr}/10)^{5.0})+4.83$ for $t_{\rm Myr}<100$ and $\psi=4.83$ for $t_{\rm 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\,\tilde{N},\,0.0118})$. 4. 4.

Radiation: We define the light-to-mass ratio in a given band $\Psi_{\rm band}$, with units $L_{\sun}/M_{\sun}$. Then the bolometric $\Psi_{\rm bol}=1136.59$ for $t_{\rm Myr}<3.5$, and $\Psi_{\rm bol}=1500\,\exp{[-4.145\,x+0.691\,x^{2}-0.0576\,x^{3}]}$ with $x\equiv\log_{10}(t_{\rm Myr}/3.5)$ for $t_{\rm 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, $\Psi_{\rm IR}=0$. In optical/NIR, $\Psi_{\rm opt}=f_{\rm opt}\,\Psi_{\rm bol}$ with $f_{\rm opt}=0.09$ for $t_{\rm Myr}<2.5$; $f_{\rm Opt}=0.09\,(1+[(t_{\rm Myr}-2.5)/4]^{2})$ for $2.5<t_{\rm Myr}<6$; $f_{\rm Opt}=1-0.841/(1+[(t_{\rm Myr}-6)/300])$ for $t_{\rm Myr}>6$. For the photo-electric FUV band $\Psi_{\rm FUV}=271\,[1+(t_{\rm Myr}/3.4)^{2}]$ for $t_{\rm Myr}<3.4$; $\Psi_{\rm FUV}=572\,(t_{\rm Myr}/3.4)^{-1.5}$ for $t_{\rm Myr}>3.4$. For the ionizing band $\Psi_{\rm ion}=500$ for $t_{\rm Myr}<3.5$; $\Psi_{\rm ion}=60\,(t_{\rm Myr}/3.5)^{-3.6}+470\,(t_{\rm Myr}/3.5)^{0.045-1.82\,\ln{t_{\rm Myr}}}$ for $3.5<t_{\rm Myr}<25$; $\Psi_{\rm ion}=0$ for $t_{\rm Myr}>25$. The remaining UV luminosity, $\Psi_{\rm bol}-(\Psi_{\rm IR}+\Psi_{\rm opt}+\Psi_{\rm FUV}+\Psi_{\rm ion})$ is assigned to the NUV band $\Psi_{\rm NUV}$. The flux-mean dust opacities adopted are $(\kappa_{\rm FUV},$ $\kappa_{\rm NUV},$ $\kappa_{\rm opt},$ $\kappa_{\rm IR})$ $=(2000,\,1800,\,180,\,10)\,(Z/Z_{\sun})\,{\rm cm^{2}\,g^{-1}}$. The photo-ionization rate (and corresponding $\kappa_{\rm 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 look-up tables in the simulations (as a function of temperature, density, and metallicity), but we provide functions accurate to $\sim 10\%$ over the relevant dynamic range in the simulations ($\sim 10-10^{9}\,$K), so that interested readers can reproduce our full cooling physics treatment.

The instantaneous cooling+heating rate per unit volume is given by the sum over all processes,

[TABLE]

where $e_{\rm thermal}^{a}$ is the thermal energy density of a gas element $a$. Here we will use $\Lambda$ to denote both cooling and heating rates, but with opposite signs (a positive sign here denotes cooling). Also, define $\tilde{n}_{x}\equiv n_{x}/n_{\rm H}$ as the number of species $x$ per hydrogen nucleus (e.g. $\tilde{n}_{e}$, $\tilde{n}_{\rm H0}$, and $\tilde{n}_{\rm HeI}$ denote the free electron, neutral hydrogen, and neutral helium numbers). Below, all units are cgs ($[T]={\rm K}$, $[n_{\rm H}]={\rm cm^{-3}}$, $[\Lambda]={\rm erg\,s^{-1}\,cm^{3}}$). The processes we track include:

Collisional Excitation: from Katz et al. (1996) (incorporating earlier fits from Cen 1992):

[TABLE] 2. 2.

Collisional Ionization: also from Katz et al. (1996):

[TABLE] 3. 3.

Recombination: from Verner & Ferland (1996):

[TABLE]

Note the $\alpha_{\rm di}$ term here comes from dielectric recombination. 4. 4.

Free-free emission: from Rybicki & Lightman (1986):

[TABLE] 5. 5.

High-Temperature Metal-Line Cooling: this refers to metal-line cooling processes in mostly ionized gas, with temperatures $\gtrsim 10^{4}\,$K. We use the public look-up tables from Wiersma et al. (2009a), for which:

[TABLE]

where 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_{\sun}^{i}$ 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\times 10^{-3}$, $1.32\times 10^{-3}$, $8.65\times 10^{-3}$, $2.22\times 10^{-3}$, $9.31\times 10^{-4}$, $1.08\times 10^{-3}$, $6.44\times 10^{-4}$, $1.01\times 10^{-4}$, $1.73\times 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 $\xi_{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_{\sun}$ (defined above).323232Available at http://www.strw.leidenuniv.nl/WSS08 6. 6.

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 $\lesssim 10^{4}\,$K. From our compilation of CLOUDY runs (Ferland et al., 1998), fitting the resulting look-up tables, we obtain approximately:

[TABLE]

where $f_{\rm selfshield}$ accounts for the local radiation environment by applying a simple (fitted) local shielding correction for UV/ionizing photons, $f_{\rm selfshield}\equiv\exp{(-\tilde{\tau}_{a}^{\rm ion})}$ with $\tilde{\tau}_{a}^{\rm ion}\equiv\sigma^{H}_{\nu_{0}}\,n_{{\rm H},\,a}\,\ell_{a}^{\rm fit}$ where $\sigma^{H}_{\nu_{0}}\equiv 6\times 10^{-18}\,{\rm cm^{-2}}$ and $\ell_{a}^{\rm fit}\equiv 4.4\,{\rm pc}\,(T/10^{4}\,K)^{-0.173}\,\Gamma_{-12}^{-2/3}$ ($\Gamma_{-12}$ is the ionization rate in units of $10^{-12}\,{\rm s^{-1}}$, including both the UV background and local sources assuming they have the same spectral shape, as defined below for photo-ionization heating). 7. 7.

Dust Collisional Heating/Cooling: from Meijerink & Spaans (2005):

[TABLE]

where we take $T_{\rm dust}=30\,$K here, and the $Z/Z_{\sun}$ term comes from assuming a constant dust-to-metals ratio. 8. 8.

Compton Heating/Cooling: from the CMB, gives (Rybicki & Lightman, 1986):

[TABLE] 9. 9.

Photo-Ionization Heating: from the UVB and local (in-simulation stellar sources) gives a heating rate, hence negative $\Lambda$, of

[TABLE]

where $\epsilon_{\rm H0}$, $\epsilon_{\rm He0}$, and $\epsilon_{\rm HeI}$ are pre-tabulated for the assumed UV background magnitude and shape, as a function of redshift in Faucher-Giguère et al. (2009)333333See http://galaxies.northwestern.edu/uvb – values above are simple polynomial fits good to $\sim 10\%$ up to $z\sim 10$. The factor $\tilde{f}$ accounts for both self-shielding (reducing the effective incident flux by $f_{\rm selfshield}$) and the contribution from local sources, where $e_{\nu,\,{\rm ion}}^{\rm local}$ is the ionizing band radiation energy density calculated explicitly from the radiation-hydrodynamic treatment in the code (Appendix E), and $e_{\nu,\,{\rm ion}}^{\rm UVB}$ is the meta-galactic UV background (UVB) energy density integrated in the same band for the same wavelength range (H-ionizing 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. 10. 10.

Cosmic Ray Heating: from Guo & Oh (2008):

[TABLE]

where we assume an approximately uniform MW-like cosmic ray background, $e_{\rm CR}\approx 9.0\times 10^{-12}\,f_{\rm CR}$. Here $f_{\rm CR}=n_{\rm H}/(0.01+n_{\rm H})$ when $n_{\rm H}$ exceeds $1000\times$ the mean baryonic density of the Universe, and $f_{\rm 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). 11. 11.

Photo-Electric Heating: from Wolfire et al. (2003):

[TABLE]

where $\tilde{e}_{\nu}^{\rm pe}$ is the photon energy density in the photo-electric band, normalized to the Habing (1968) MW units, $\tilde{e}_{\nu}^{\rm pe}\equiv e_{\nu}^{\rm pe}/(3.9\times 10^{-14}\,{\rm erg\,cm^{-3}})$. The $Z/Z_{\sun}$ term comes from assuming a constant dust-to-metals ratio. Here the field $e_{\nu}^{\rm pe}$ is the FUV band radiation energy density calculated explicitly from the radiation-hydrodynamic treatment in the code, described in Appendix E. 12. 12.

Magneto-Hydrodynamic Work & Shocks: from the MHD equations, 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:

[TABLE]

where $u_{\rm thermal}$ is the specific internal energy (internal energy per unit mass). 13. 13.

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 $\Lambda_{\rm Net}$, and then restricting this to the cooling rate of said slab:

[TABLE]

where $\Sigma_{\rm eff}=\langle\Sigma_{\rm column}^{a,\,{\rm Sobolev}}\rangle_{\phi,\,\theta}=\rho_{a}\,(h_{a}+\rho_{a}/|\nabla\rho_{a}|)$ uses the local Sobolev approximation to estimate the column density to infinity and is defined in Appendix E, and $\kappa_{\rm eff}$ is the effective opacity.343434We take $\kappa_{\rm 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 $\kappa_{\rm eff}=5$ for $150\leq T\leq 1500$ and $\kappa_{\rm eff}=0.0027\,T^{3/2}$ for $T<150$. At higher temperatures the system is rarely optically thick, but for completeness we compute $\kappa_{\rm eff}$ from the gas-phase using standard approximations for stellar atmospheres: $\kappa_{\rm eff}^{-1}=\kappa_{\rm rad}^{-1}+\kappa_{\rm cond}^{-1}$, with $\kappa_{\rm cond}=2.6\times 10^{-7}\,\tilde{n}_{e}\,T^{2}\,\rho^{-2}$, $\kappa_{\rm rad}=\kappa_{\rm mol}+1/(\kappa_{H^{-}}^{-1}+[\kappa_{e}+\kappa_{\rm Kr}]^{-1})$, $\kappa_{\rm mol}=0.1\,Z$, $\kappa_{e}=0.2\,(1+X_{H})$, $\kappa_{H^{-}}=1.1\times 10^{-25}\,(Z\,\rho)^{1/2}\,T^{7.7}$, $\kappa_{\rm Kr}=4.0\times 10^{25}\,(1+X_{H})\,Z\,\rho\,^{3.5}$.

As 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

Self-Gravitating: First, following standard sink-particle approaches, we calculate the virial parameter $\alpha$ (ratio of thermal plus kinetic energy to potential energy) in a resolution element and allow only star formation in bound particles with $\alpha<1$. From Hopkins et al. (2013d),

[TABLE]

where $\|\nabla\otimes{\bf v}\|_{a}$ the Frobenius norm ($\|{\bf A}\|^{2}\equiv\sum_{\alpha\beta\gamma...}\,A^{2}_{\alpha\beta\gamma...}$) of the velocity gradient tensor ($\otimes$ is the outer product), $c_{s,\,a}$ is the sound speed, $\rho_{a}$ the density, and $h_{a}$ the usual resolution scale (inter-particle 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). Using the zeroth-order inaccurate SPH gradient estimator, in contrast, gives similar results statistically, but makes identification of individual physically collapsing clouds much more noisy. 2. 2.

Self-Shielding: If $\alpha<1$, we next calculate the shielded/molecular fraction $f_{\rm shielded}^{({\rm sf})}$, 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:

[TABLE]

We require $f_{\rm shielded}^{({\rm sf})}>0$ for star formation. 3. 3.

Jeans-Unstable: If $\alpha<1$ and $f_{\rm shielded}^{({\rm sf})}>0$, we calculate the Jeans mass $m_{J}$, and only allow star formation in Jeans-unstable particles, specifically those where $m_{J}<m_{J,\,{\rm crit}}\equiv{\rm MAX}(m_{a},\,10^{3}\,M_{\sun})$, where $m_{a}$ is the particle mass. We calculate the Jeans mass as

[TABLE]

where $n_{a}\equiv\rho_{a}/\mu_{a}$ is the gas number density. 4. 4.

Sufficiently-Dense: If $\alpha<1,\,f_{\rm shielded}^{({\rm sf})}>0$, and $m_{J}<m_{J,\,{\rm crit}}$, we check if $n_{a}>n_{\rm crit}$, where $n_{\rm crit}=1000\,{\rm cm^{-3}}$ is a minimum density (and $n_{a}\equiv n_{a,\,H}=X_{H}\,\rho_{a}/m_{p}$), to prevent spurious triggering of the above criteria in low-density gas.

We then assign the gas particle a volume-integrated SFR:

[TABLE]

Because we wish to maintain equal stellar and gas element masses, at each timestep $\Delta t_{a}$ we assign the gas particle a probability $p_{a}=1-\exp{(-\dot{m}_{\ast}^{a}\,\Delta t_{a}/m^{a}_{\rm gas})}$ of turning into a star particle that timestep; we draw a uniform random variable $0<x<1$ and if $x<p_{a}$, 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).

In Paper II, we discuss each piece of this algorithm in detail, and 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.

Every timestep $\Delta t_{a}$, for each star particle $a$ (at position ${\bf 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_{\rm ej}$, metal mass $m_{Z,\,{\rm ej}}$ (defined for each species we track), momentum $p_{\rm ej}=m_{\rm ej}\,v_{\rm ej}$, and energy $E_{\rm ej}$. These also are given in Appendix A (for winds, $m_{\rm ej}=\Delta t\,\dot{M}_{w}$ from the star, for SNe it is the ejecta mass). 2. 2.

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_{\rm ngb}=(4\pi/3)\,H_{a}^{3}\,\bar{n}_{a}(H_{a})$ in the same manner as for hydrodynamics, where $W$ is the kernel function, $\bar{n}_{a}=\sum W({\bf x}_{ba}\equiv{\bf x}_{b}-{\bf x}_{a},\,H_{a})$, and $H_{a}$ is the search radius. The equation for $N_{\rm ngb}(H_{a})$ is non-linear so is solved iteratively in the neighbor search; see Springel (2005). Thus we obtain all gas elements $b$ within a radius $|{\bf x}_{ba}|<H_{a}$ (where “$a$ sees $b$”); we also identify all neighbors with $|{\bf x}_{ba}|<H_{b}$ (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. 3. 3.

Assign “vector weights” to each neighbor, by first boosting to the rest-frame of the star (${\bf x}_{a}={\bf 0}$, ${\bf v}_{a}\equiv d{\bf x}_{a}/dt={\bf 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 $\bar{\bf w}_{ba}$

[TABLE]

where ${\bf A}_{ba}$ is the effective vector face between elements $b$ and $a$ used in the finite-volume hydrodynamic calculations,353535For our MFM hydrodynamic method, the face ${\bf A}_{ba}$ is defined as (see Hopkins 2015):

$\displaystyle{\bf A}_{ba}$ $\displaystyle\equiv\bar{n}_{a}^{-1}\,\bar{\bf q}_{b}({\bf x}_{a})+\bar{n}_{b}^{-1}\,\bar{\bf q}_{a}({\bf x}_{b})$

(58)

$\displaystyle\bar{\bf q}_{b}({\bf x}_{a})$ $\displaystyle\equiv{\bf E}_{a}^{-1}\cdot{\bf x}_{ba}\,W({\bf x}_{ba},\,H_{a})$

(59)

$\displaystyle{\bf E}_{a}$ $\displaystyle\equiv\sum_{c}\,({\bf x}_{ca}\otimes{\bf x}_{ca})\,W({\bf x}_{ca},\,H_{a})$

(60)

For SPH, the face is defined by the simpler relation ${\bf A}_{ba}=[\bar{n}_{a}^{-2}\,\partial W(|{\bf x}|_{ba},\,H_{a})/\partial|{\bf x}|_{ba}+\bar{n}_{b}^{-2}\partial W(|{\bf x}|_{ba},\,H_{b})/\partial|{\bf x}|_{ba}]\ \hat{\bf x}_{ba}$. In moving-mesh or fixed-grid finite-volume codes, the face ${\bf A}_{ba}$ is the explicit geometric mesh face between cells. and the $\hat{\bf x}_{ca}^{\pm}$ are the positive or negative (singly-signed) projection vectors:

[TABLE]

These 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 $\bar{\bf 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. 4. 4.

Assign initial fluxes in the rest-frame of the star:

[TABLE]

It is easy to see that our definition of $\bar{\bf w}_{ba}$ guarantees exact conservation of mass, energy, and linear momentum, and that the correct total radial (outward) momentum is assigned, e.g.:

[TABLE] 5. 5.

Boost back to the simulation (“lab”) frame: if the star is moving with velocity ${\bf v}_{a}$, then this boost transforms the momentum and energy fluxes:

[TABLE]

(the mass fluxes are unchanged, $\Delta m_{b}^{\prime}=\Delta m_{b}$, $\Delta m_{Z,\,b}^{\prime}=\Delta m_{Z,\,b}$). Of course this maintains manifest conservation: the total momentum added to the neighbors via the $\Delta m_{b}\,{\bf v}_{a}$ term exactly cancels that lost by the star, since its mass decreases by $\sum\Delta m_{b}=m_{\rm ej}$. 6. 6.

Account for $PdV$ (mechanical) work: consider that we have a particle $b$ representing a volume domain with mass $m_{b}$ around our source, with some mean distance in the volume element $|{\bf x}_{ba}|$ (which we call the “coupling radius”). The ejecta, in order to reach that point, must sweep up the mass $m_{b}$ (in e.g. a shock or shell) – it 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 ($\Delta{\bf p}_{b}^{\prime\prime}$) to couple into the domain $b$ is not the initial ejecta momentum $\Delta{\bf p}_{b}^{\prime}$. 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 $\Delta{\bf p}_{b}^{\prime}\,(1+m_{b}/\Delta m_{b})^{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_{\rm t}$ (and becomes momentum, rather than energy conserving). Therefore we must impose an upper limit $\Delta{\bf p}_{b}^{\prime}\,p_{\rm t}/p_{\rm ej}$. We therefore have:

[TABLE]

The expression for $p_{\rm 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_{\rm t}$ much larger than its actual physical uncertainty. It is easy to verify, given the form of Eq. 75, that at sufficiently high resolution ($m_{b}\ll 1000\,M_{\sun}$), the $p_{\rm 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 $\Delta E$ represents the total energy, this is not directly modified by changing $\Delta{\bf p}$ (the correct thermal-kinetic breakdown is automatic).

Also note that Eq. 75 is an approximation if the gas surrounding the star particle is moving at a non-uniform velocity (with non-negligible velocities relative to the ejecta); in this limit the exact expression is given in Paper II (Appendix E). 7. 7.

Add final fluxes to the neighboring gas elements, in a fully-conservative manner:

[TABLE]

So (like our hydrodynamic update), we add conserved quantities ($m$, ${\bf p}$, $E$) and from those update primitive quantities ($Z$, ${\bf 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.

To be fully consistent with the radiative losses described above (when the cooling radius is un-resolved), we must also modify $U_{\rm internal}^{\rm new}$. Following Thornton et al. (1998), the thermal post-shock energy outside $R_{\rm cool}$ decays rapidly as $\propto(r/R_{\rm cool})^{-6.5}$; so we estimate the effective $R_{\rm cool}$ via the requirement that, at the end of the energy-conserving phase, $(1/2)\,(m_{\rm ej}+m_{\rm swept}[R_{\rm cool}])\,v_{f}^{2}=(1/2)\,m_{\rm ej}\,v_{\rm ej}^{2}$ and $p_{\rm t}=m_{\rm swept}[R_{\rm cool}]\,v_{f}$ (where $m_{\rm swept}$ is the enclose mass “swept up” by the shell), giving $R_{\rm cool}\approx 28.4\,{\rm pc}\,(n_{b}/{\rm cm^{-3}})^{-3/7}\,(E_{\rm tot,\,ej}/10^{51}\,{\rm erg})^{2/7}\,f(Z_{b})$ for $p_{\rm t}$ in Eq. 76. If $|{\bf x}_{ba}|<R_{\rm cool}$, we leave $U_{\rm internal}^{\rm new}$ un-modified. If $|{\bf x}_{ba}|>R_{\rm cool}$, we calculate the increase in internal energy from shock-heating, ignoring cooling: $\Delta U_{b}\equiv U_{\rm internal}^{\rm new}-U_{\rm internal}^{\rm old}$, and then modify it to determine the correct internal energy: $U_{\rm internal}=U_{\rm internal}^{\rm old}+\Delta U_{b}\,(|{\bf x}_{ba}|/R_{\rm 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_{\rm cool}$, it will radiate its energy rapidly, so in practice we find that if we simply leave $\Delta U_{b}$ 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), photo-ionization, and photo-electric heating. The algorithm here was first developed in Hopkins et al. (2012c), and the version used in FIRE-1, which is almost exactly identical to that here,363636The 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 FIRE-1, 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. 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.

We 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 $\sim 10^{6}\,{\rm cm^{-3}}$ and size scales $<0.1\,$pc).

In 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.

The complete algorithm is:

Determine Background Radiation (Source Luminosities): every timestep $\Delta t_{a}$ for each star particle $a$, we take the luminosity $L_{\nu}^{a}=\Psi_{\nu}^{a}\,m_{\ast,\,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_{\rm ion}$, $\lambda<912\,$Å), far-UV ($L_{\rm FUV}$, $912\,$Å$<\lambda<1550\,$Å), near-UV ($L_{\rm UV}$, $1550<\lambda<3600\,$Å), optical/near-IR ($L_{\rm Opt}$, $3600\,$Å$<\lambda<3\,\mu$), and mid/far-IR ($L_{\rm IR}$, $\lambda>3\,\mu$). 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. 2. 2.

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 $\tau_{\nu}^{a}=\kappa_{\nu}\,\Sigma_{\rm column}^{a}$, where $\kappa_{\nu}$ 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 $\Sigma_{\rm column}^{a}$ using the Sobolev approximation for the isotropic (angle-averaged) column density integrated outward from the source:

[TABLE]

Here $\rho_{a}$, $\nabla\rho_{a}$, and $h_{a}$ are the gas density, density gradient and inter-element spacing, evaluated at location ${\bf x}_{a}$ (with the same algorithm as our usual hydrodynamics). The $\rho_{a}/|\nabla\rho_{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 $\rho_{a}\,h_{a}$ term is just the column through the local cell. Since this is isotropic, the absorbed luminosity (in a narrow band) is just $L_{{\rm abs},\,\nu}^{a}=(1-\exp{(-\tau_{\nu}^{a})})\,L_{\nu}^{a}$ and the surviving un-absorbed luminosity is:

[TABLE]

Recall 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:

[TABLE]

For 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 $\dot{N}_{a}^{\rm ion}$, the rate of production of ionizing photons ($\propto L_{\rm ion}$). Now, we take all gas elements $b$ in the vicinity of $a$ and sort them by increasing distance $|{\bf x}_{ba}|$. Beginning with the closest, we test whether it is already ionized (either because $T_{b}>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: $\Delta\dot{N}_{b}=N(H)_{b}\,\beta\,n_{e,\,b}$ (where $N(H)_{b}=X_{H}\,m_{b}/\mu\,m_{p}$ is the number of H atoms in $b$, $\beta\approx 3\times 10^{-13}\,{\rm cm^{3}\,s^{-1}}$ is the recombination coefficient, and $n_{e,\,b}$ is the electron density assuming full ionization). If $\Delta\dot{N}_{b}\leq\dot{N}_{a}^{\rm ion}$, then particle $b$ is tagged as being within an HII region, and the photons are “consumed,” so $\dot{N}_{a}^{\rm ion}\rightarrow\dot{N}_{a}^{\rm ion}-\Delta\dot{N}_{b}$. We then proceed to the next particle and repeat. If we reach a particle which is not ionized but for which $\Delta\dot{N}_{b}>\dot{N}_{a}^{\rm ion}$, we determine whether or not to ionize it randomly, with probability $p=\dot{N}_{a}^{\rm ion}/\Delta\dot{N}_{b}$, 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. 85. Tests of this algorithm (both static but also dynamic tests of D-Type ionization front expansion) are also shown in Hu et al. (2017). 3. 3.

Account for Momentum of Locally-Absorbed Photons: Over a timestep $\Delta 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

[TABLE]

This momentum flux is distributed among the neighbors, directed radially away from the star particle, as described in Appendix D. 4. 4.

Transport the Locally-Extincted Radiation via an Optically-thin Network: We now have an “emergent” spectrum after local attenuation around each star, $L_{{\rm emergent},\,\nu}$. 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 ${\bf F}_{\nu}^{b}$ and photon energy density $e_{\nu}^{b}$ at a distant gas element $b$ are just

[TABLE]

This 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).373737The softening kernel, following Hopkins (2015), is given by replacing Eqs. 88-89 with ${\bf F}_{\nu,\,{\rm emergent}}=\sum_{a}\,(1/4\pi)\,L_{{\rm emergent},\,\nu}^{a}\,({\bf x}_{b}-{\bf x}_{a})\,H_{a}^{-3}\,\mathcal{F}_{s}(u_{ba})$ and $\epsilon_{\nu}^{b}=\sum_{a}\,(1/4\pi\,c)\,L_{{\rm emergent},\,\nu}^{a}\,|{\bf x}_{b}-{\bf x}_{a}|\,H_{a}^{-3}\,\mathcal{F}_{s}(u_{ba})$ where $u_{ba}\equiv|{\bf x}_{b}-{\bf x}_{a}|/H_{a}$, $H_{a}=(24/\pi)^{1/3}\,h_{a}=(24/\pi)^{1/3}\,\epsilon_{a}$ is the maximum kernel search radius, and

$\displaystyle\mathcal{F}_{s}(u)$ $\displaystyle\equiv\begin{cases}{\displaystyle\frac{32}{15}\,\left[5+3\,u^{2}\,(5\,u-6)\right]\ \ \ \ \ \ \ \ \hfill{\left(u\leq\frac{1}{2}\right)}}\\ \\ {\displaystyle\frac{32}{15}\,\left[10-\frac{45}{2}\,u+18\,u^{2}-5\,u^{3}-\frac{1}{32\,u^{3}}\right]\ \ \ \hfill{\left(\frac{1}{2}<u<1\right)}}\\ \\ {\displaystyle\frac{1}{u^{3}}\ \ \ \ \ \ \ \ \hfill{\left(u\geq 1\right)}}\end{cases}$

(90)

This becomes exactly inverse-square at $r>H$, and but prevents a $1/r^{2}$ divergence as $r\rightarrow 0$ (with flux $\sim 1/h_{a}^{2}$ instead of $1/r^{2}$). 5. 5.

Calculate incident radiative acceleration from long-range fluxes: For a gas element $b$ with effective face area $A_{b}$ and mass $m_{b}$ (hence surface density $\Sigma_{b}\equiv m_{b}/A_{b}$) seeing an incident flux ${\bf F}_{\nu}^{b}$, the exact radiative acceleration is given by

[TABLE]

In the optically thin limit ($\kappa_{\nu}\,\Sigma_{b}\ll 1$), this reduces to the common expression ${\bf a}=\kappa_{\nu}\,{\bf F}_{\nu}/c$, but in the optically thick limit ($\kappa_{\nu}\,\Sigma_{b}\gg 1$), the force saturates at $m_{b}\,{\bf a}=({\bf F}_{\nu}\,A_{b})/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 ($=m_{b}/\rho_{b}=(4\pi/3)\,h_{b}^{3}$), i.e. $\pi\,h_{b}^{2}$; using the more complicated hydrodynamic face areas introduces negligible ($\ll 10\%$) differences in the accelerations here. We adopt $\kappa_{\nu}=\langle\kappa\rangle_{\nu}$ for each band. 6. 6.

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 $\langle e^{b}_{\nu}\rangle=e^{b}_{\nu}\,\exp{(-\tau_{\nu}^{b})}$, where $e^{b}_{\nu}$ is the photon energy density given by Eq. 89, and $\tau_{\nu}^{b}$ is the optical depth estimated using the Sobolev approximation in Eq. 82, but now at the location of the absorbing gas element (instead of around the emitting star). The resulting, shielded $e_{\nu}^{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.

F.1 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 fluid-mixing 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.

Readers interested in more details of our numerical implementation of MHD should consult Hopkins & Raives (2016) and the public GIZMO source code.

F.2 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 $\boldsymbol{\Pi}$ to the momentum flux (${\bf F}_{\bf p}=\boldsymbol{\Pi}$) and energy flux (${\bf F}_{e}=\boldsymbol{\Pi}\cdot{\bf v}$), and a conductive energy flux ${\bf F}_{e}={\bf K}\cdot\nabla T$. The appropriate anisotropic tensor expressions for MHD are given by (Spitzer & Härm, 1953; Braginskii, 1965):

[TABLE]

where $\otimes$ denotes the outer product; $\hat{B}$ is the direction of the magnetic field vector; ${\bf I}$ is the identity matrix; ${\bf v}$ the velocity; “$:$” denotes the double-dot-product (${\bf A}:{\bf B}\equiv{\rm Trace}({\bf A}\cdot{\bf B})$); $\ln{\Lambda}\approx 37.8$ is the Coulomb logarithm (Sarazin, 1988); $m_{e}$, $e$, $m_{i}$, $Z_{i}\,e$ are the electron mass and charge and ion mass and charge; $f_{i}$ the ionized fraction (calculated self-consistently in our cooling routines); $k_{B}$ the Boltzmann constant; $\ell_{e}$ ($\ell_{i}$) is the electron (ion) mean-free path, and $\ell_{T}=T/|\nabla T|$ ($\ell_{v}=|{\bf v}|/||\nabla\otimes{\bf v}||$) is the temperature (velocity) gradient scale length (this term correctly accounts for saturation of $\kappa_{\rm cond}$ or $\nu_{\rm visc}$ when electrons/ions have long mean free paths, by not allowing the gradient scale length to be shorter than $\ell_{e,\,i}$). In these equations, $\kappa_{\rm cond}$ is the conductivity, and $\nu_{\rm 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 fixed-grid codes, and is able to correctly capture complicated non-linear instabilities sourced by anisotropic diffusion such as the magneto-thermal and heat-flux bouyancy instabilities.

Readers interested in more details of our numerical implementation of anisotropic diffusion should consult Hopkins (2017) and the public GIZMO source code.

F.3 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: $\partial(\rho\,Z)/\partial t=\nabla\cdot(\kappa_{\rm turb}\,\rho\,\nabla Z)$, where $Z$ is the abundance per unit mass of the scalar (i.e. the metallicity) and the “eddy diffusivity” $\kappa_{\rm turb}\sim\lambda_{\rm eddy}\,v_{\rm 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 $\lambda_{\rm eddy}$ scales with the eddy turnover time. Following the common Smagorinsky (1963) approximation, we can approximate the “eddy diffusivity” as

[TABLE]

where $C\sim 0.05-0.15$ is a constant calibrated to numerical simulations, motivated by a Kolmogorov cascade in Smagorinsky (1963), $h_{a}$ is the grid scale (for our MFM method, this is equal to the rms inter-element separation), and ${\bf S}_{a}\equiv[(\nabla\otimes{\bf v})_{a}+(\nabla\otimes{\bf v})_{a}^{T}]-{\rm Trace}(\nabla\otimes{\bf v})_{a}/3$ is the symmetric shear tensor (and $\|{\bf S}\|$ denotes the Frobenius norm). Note that we use our higher-order matrix-based gradient formalism from Hopkins (2015) to calculate ${\bf 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).

The 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. 96 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, high-resolution supersonic turbulent box simulations (with and without magnetic fields and/or shear), and verify that this prescription, with $C\approx 0.05$, is reasonable specifically in our identical MFM code with the definitions of $h$ and ${\bf 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\approx 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.

We stress that a term like Eq. 96 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 $\sim h_{a}\,\Delta v(h_{a})$ 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. 96) 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. 96. 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 small-scale 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. 96 as a purely numerical term which “restores” the desirable aspect of the numerical diffusivity present in certain numerical methods.

We caution, however, that simple diffusion prescriptions such as Eq. 96, naively applied, can substantially over-estimate the diffusivity. The critical assumption is that the resolution-scale motion $\Delta v\sim\|{\bf S}\|\,h\sim v_{\rm eddy}$ is entirely due to turbulence; if there is any bulk motion included in ${\bf S}$, this will over-estimate $\kappa_{\rm 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 $\kappa_{\rm turb}$ can over-estimate by an order of magnitude the true turbulent motion. For example, if the disk scale height $H_{\rm disk}\approx\sigma_{\rm turb}/\Omega$ is unresolved, $h_{a}\gg H_{\rm disk}$, then Eq. 96 will return $\kappa_{\rm turb}\sim h_{a}\,(h_{a}\,\Omega)$ instead of the correct maximum diffusivity for disk-scale eddies, $\kappa_{\rm turb}\sim H_{\rm disk}\,\sigma_{\rm turb}\sim H_{\rm disk}^{2}\,\Omega$. 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 finite-element 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. 96. 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. 96) 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.

Readers 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.

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 $\approx 40\%$ in previously-published work to $<10\%$ in the current code. More generally, the entire GIZMO code has been line-by-line optimized, manually unrolling or (where possible) vectorizing 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 $\sim 2$ speed improvement. 2. 2.

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 $\sim 2-3$, but extending the strong scaling to $\sim 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 $\sim 2-3$ reduction in load imbalances at the smallest timesteps. 3. 3.

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. 4. 4.

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). 5. 5.

We adopt adaptive gravitational softenings as described above. This imposes essentially no cost for gas (since $h_{i}$ 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). 6. 6.

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). 7. 7.

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 $\sim 1.3-2.5$, mostly owing to the larger neighbor number needed in P-SPH to achieve comparable accuracy. 8. 8.

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\sim 2$, because the extremely dense relic star clusters would require constant, extremely short timesteps.

Readers interested in more details should consult the public source code.

Figure 38 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.

Our optimizations allow us to extend good strong scaling, at our modest “typical FIRE-1” resolution, to $\sim 1024$ cores. Even more strikingly, the optimizations we have made allow us to maintain good weak scaling up to $\sim 4096$ cores for a simulation with $3\times 10^{8}$ particles – our “Latte” resolution – and $\sim 16,384$ cores for a zoom-in simulation with $>10^{9}$ particles (and as many as $\sim 10^{6}$ cores for large-volume simulations with $\sim 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.

For comparison, if we compare similar runs to the strong-scaling tests in Fig. 38 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 $\sim 64-128$ cores.

Appendix H Congratulations!

You made it to the end of this paper. Fig. 39 is your reward. Go Cleveland!