Feedback with Individual Stars

Simulating an Isolated Dwarf Galaxy with Multi-Channel Feedback and Chemical Yields from Individual Stars

Andrew Emerick, Greg L. Bryan, and Mordecai-Mark Mac Low
Department of Astronomy, Columbia University, New York, NY, 10027, USA
Department of Astrophysics, American Museum of Natural History, New York, NY, USA
Center for Computational Astrophysics, Flatiron Institute, 162 5th Ave, New York, NY, 10003, U.S.A

In order to better understand the relationship between feedback and galactic chemical evolution, we have developed a new model for stellar feedback at grid resolutions of only a few parsecs in global disk simulations, using the adaptive mesh refinement hydrodynamics code Enzo. For the first time in galaxy scale simulations, we simulate detailed stellar feedback from individual stars including asymptotic giant branch winds, photoelectric heating, Lyman-Werner radiation, ionizing radiation tracked through an adaptive ray-tracing radiative transfer method, and core collapse and Type Ia supernovae. We furthermore follow the star-by-star chemical yields using tracer fields for 15 metal species: C, N, O, Na, Mg, Si, S, Ca, Mn, Fe, Ni, As, Sr, Y, and Ba. We include the yields ejected in massive stellar winds, but greatly reduce the winds’ velocities due to computational constraints. We describe these methods in detail in this work and present the first results from 500 Myr of evolution of an isolated dwarf galaxy with properties similar to a Local Group, low-mass dwarf galaxy. We demonstrate that our physics and feedback model is capable of producing a dwarf galaxy whose evolution is consistent with observations in both the Kennicutt-Schmidt relationship and extended Schmidt relationship. Effective feedback drives outflows with a greater metallicity than the ISM, leading to low metal retention fractions consistent with observations. Finally, we demonstrate that these simulations yield valuable information on the variation in mixing behavior of individual metal species within the multi-phase interstellar medium.

Galaxy – Hydrodynamics – Simulation – Feedback

1 Introduction

Detailed interstellar medium (ISM) and chemical abundance properties of galaxies are sensitive tests of the underlying physical processes that govern galaxy evolution. Examining these in more detail in galaxy scale simulations is an important and exciting new discriminator between models. There is a considerable body of work studying the chemodynamical evolution of galaxies using cosmological hydrodynamics simulations (e.g. Oppenheimer & Davé, 2008; Wiersma et al., 2009; Shen et al., 2010; Simpson et al., 2013; Snaith et al., 2015; Schaye et al., 2010, 2015; Hopkins et al., 2014). These simulations, coupled with additional attention to feedback processes, have made remarkable progress in reproducing global galaxy trends such as evolution of the mass-metallicity relationship (e.g. Obreja et al., 2014; Ma et al., 2016a; Davé et al., 2017; Torrey et al., 2017) and more detailed quantities such as metallicity distribution functions (MDFs) and the evolution of individual species abundances (Marcolini et al., 2008; Revaz et al., 2009; Sawala et al., 2010; Revaz & Jablonka, 2012; Jeon et al., 2017; Hirai & Saitoh, 2017) .

However, much of this work has been done with Lagrangian smoothed particle hydrodynamics schemes, with a few recent exceptions (Few et al., 2012; Simpson et al., 2013; Few et al., 2014; Vorobyov et al., 2015; Corlies et al., 2018). In its original form, this scheme does not capture mixing between chemically inhomogeneous particles, as necessary for chemical evolution. Mixing can be modeled with sub-grid models of turbulent metal diffusion (e.g. Shen et al., 2010; Shen et al., 2013; Brook et al., 2014; Su et al., 2017b; Escala et al., 2018), but there are many possible models and each is not necessarily applicable in every regime (see Revaz et al., 2016). While mixing occurs in Eulerian codes even without sub-grid models, numerical diffusion tends to result in over-mixing in simulations without sufficiently high spatial resolution. Molecular diffusion or even turbulent mixing is certainly not resolved in any galaxy-scale simulation with either method, requiring additional sub-grid models; this can be particularly important for understanding the initial pollution of otherwise pristine gas (see Pan et al., 2013; Sarmento et al., 2017, and references therein). Moreover, metal mixing efficiencies may vary species-by-species (e.g. Cohen et al., 2013; Roederer et al., 2014; Frebel & Norris, 2015; Hirai & Saitoh, 2017; Côté et al., 2017a; Krumholz & Ting, 2018). Mixing behavior is tied critically to the feedback source (stellar winds, supernovae, and possibly more exotic sources) that inject metals into different phases of the ISM with different energies and on different timescales; the observational effect of this is poorly understood, however. The variations in how different methods handle sub-grid metal injection and metal mixing schemes can lead to uncertainties in connecting models to observations and the fundamental physics that drives galaxy evolution.

Increasing physical resolution reduces reliance on sub-grid physics for mixing. However, at high particle mass resolution ( M) standard schemes for modeling stars as simple stellar populations lose validity (as studied in detail by Revaz et al., 2016). Below 10 M, such schemes do not fully sample the initial mass function (IMF), and cannot be considered average representations of stellar populations. This is acutely problematic at low star formation rate densities with star particle masses comparable to or below the mass of the most massive individual star ( M). At high star formation rate densities, an undersampled IMF in a single low mass star particle can be compensated for by having many adjacent star particles. Various approaches exist to address this issue (e.g. Kobayashi et al., 2000; Weidner & Kroupa, 2004; Pflamm-Altenburg & Kroupa, 2006; Revaz & Jablonka, 2012; Kroupa et al., 2013; Rosdahl et al., 2015; Su et al., 2017a), but none are without caveats (Revaz et al., 2016), save for schemes which begin to track the star-by-star information within a given particle by directly sampling the IMF at formation time (e.g. Hu et al., 2017). The most straightforward solution is to remove the single stellar population formalism entirely and simply track stars as individual particles.

We introduce a new method for studying galactic chemical evolution that follows stars as individual star particles implemented in the adaptive mesh refinement code Enzo, designed for high resolution simulations of isolated galaxies. The relative simplicity of idealized, isolated galaxy evolution simulations allows for a focused, first-principles approach to studying multi-channel feedback mechanisms. We follow recent work using low mass dwarf galaxies as laboratories to study in detail how feedback governs galaxy evolution (Forbes et al., 2016; Hu et al., 2016, 2017). Our work builds upon our current understanding of feedback and galactic chemodynamics while making three notable advances: 1) direct star-by-star modeling, 2) stellar winds from both massive and asymptotic giant branch (AGB) stars, and 3) using an adaptive ray tracing method to follow stellar ionizing radiation. We also include core collapse and Type Ia supernova feedback, photoelectric heating from stellar far ultra-violet (FUV) radiation, and Lyman-Werner dissociation from stellar radiation.

Using star-by-star modeling, we capture in more detail the stellar yields from individual stars released over their lifetime. This includes yields from massive and AGB stellar winds, and supernovae (SNe). In addition to better capturing how individual metal species enrich the ISM, this allows us to chemically tag individual stars. This ability opens an exciting new channel for testing models of galaxy evolution by leveraging current and ongoing observations probing the detailed distributions of chemical abundances of stars in the Milky Way and Local Group, such as APOGEE and APOGEE2 (Majewski et al., 2010, 2016), the Gaia-ESO survey (Gilmore et al., 2012), and GALAH (De Silva et al., 2015). This paper is the first in a series examining in detail the role that individual components of multi-channel stellar feedback play in galaxy dynamical and chemical evolution. We describe each mode of our multi-channel feedback in detail in Section 2, describe their implementation in an isolated dwarf galaxy simulation in Section 3, show results from this simulation in Section 4, discuss the results in Section 5, and conclude in Section 6.

Readers who may want to only briefly skim (or skip) over the details of our included physics are advised to just read the beginning of Section 2, which contains a complete—yet brief—summary of the included physics.

2 Methods

We produce high-resolution, galaxy-scale simulations tracking stars not as single stellar populations, but as individual stars sampled from an assumed IMF. This allows us to follow star-by-star variations in feedback physics and stellar yields in detail. To properly model the ISM, we track non-equilibrium, primordial chemistry (including molecular hydrogen) using Grackle (Smith et al., 2017), with heating and approximate self-shielding from a metagalactic ultraviolet (UV) background. We assume collisional ionization equilibrium for all other elements and use updated Cloudy metal-line cooling tables consistent with our self-shielding approximation (see Appendix D). We also include an updated observationally motivated dust model for the low metallicity regimes studied here ( Z). Each star is assigned basic properties including surface gravity, effective temperature, radius, and lifetime from tabulated stellar evolution models, which inform how the stars deposit their feedback. We directly track ionizing radiation from massive stars using an adaptive ray tracing radiative transfer method that includes the effects of radiation pressure on HI. In addition, we follow the optically thin, non-ionizing UV radiation from these stars that cause photoelectric heating and Lyman-Werner dissociation of molecular hydrogen. We track the stellar wind feedback and SNe from these stars, depositing individual metal yields from both. We include AGB wind feedback and yields for lower mass stars, and track these directly as Type Ia SN progenitors. We follow yields for 15 individual metal species (C, N, O, Na, Mg, Si, S, Ca, Mn, Fe, Ni, As, Sr, Y, and Ba), chemically tagging each star as it forms with the associated local gas abundances for each species. In addition, we track a total metal density field which is the sum of all metals, including those not directly tracked. This field is used to inform the heating/cooling physics, and determines the metallicity of each star at birth. These methods are discussed in full detail below.

2.1 Hydrodynamics and Gravity

We use the adaptive mesh refinement hydrodynamics and N-body code Enzo111 to simulate the chemodynamical evolution and detailed feedback physics in a set of high resolution, isolated, low-mass dwarf galaxies. Enzo is an open-source code that is undergoing continual, active development by many researchers across several institutions. We use a substantially modified version of the current development version of Enzo (version 2.X) in this work.222This version is contained in a publicly available fork of the main repository: Specifically, simulations presented here were conducted at changeset 7001d99. We solve the equations of hydrodynamics using a direct-Eulerian piecewise parabolic method (Colella & Woodward, 1984; Bryan et al., 1995) and a two-shock approximate Riemann solver with progressive fallback to more diffusive Riemann solvers in the event that higher order methods produce negative densities or energies. We compute the total gravitational potential from gas self-gravity, stars, and a static background dark matter potential (see Section 3). Self-gravity is computed with a multigrid Poisson solver. The collisionless star particles are evolved with an adaptive particle-mesh N-body solver at an effective force resolution of , where is the local cell size.

We refine the mesh whenever the thermal Jeans length is no longer resolved by a minimum of 8 cells, continually refining a given region until this criterion is met or until the region reaches the maximum resolution (1.8 pc). At maximum resolution, the Jeans length can become under-resolved, leading to artificial numerical fragmentation. Truelove et al. (1997) showed that resolving the Jeans length by at least 4 cells is required to suppress this fragmentation.

We set the star formation density threshold to the value at which the Jeans length becomes resolved by only 4 cells in sub-200 K gas, or about 200 cm(as discussed further in Section 2.3). Forming stars from this gas will reduce the local density, ensuring the Jeans length is resolved. However, since star formation is not instantaneous, we employ a pressure floor to support gas against artificial fragmentation and collapse. A non-thermal pressure term is added to cells once their thermal Jeans length becomes unresolved. This prevents dense, self-gravitating gas from rapidly reaching densities significantly above our resolution limit. The use of a pressure floor is common in galaxy scale simulations with limited dynamic range (e.g Machacek et al., 2001; Robertson & Kravtsov, 2008).

Due to computational constraints we found it necessary to institute a temperature ceiling in low density, diffuse gas. These high temperatures, typically well above 10 K and up to 10 K, would place an onerous constraint on the limiting time step at high spatial resolution. At these temperatures, with typical velocities up to 10 km s, satisfying the Courant condition requires time steps on order of 100 yr. We institute a maximum temperature of 8 K in gas with densities between 10 g cm and 10 g cm. These densities were somewhat arbitrarily chosen, but ensure that this threshold does not impact very low density gas in the halo of our dwarf galaxy or higher density gas in supernova injection regions. This threshold decreases the required run-time by factors of a few. The value of the temperature threshold was chosen to ensure the affected hot gas remained just above the high-temperature minimum of our cooling curve (see Appendix D.)

2.2 Chemistry and Cooling Physics

We use the chemistry and cooling library Grackle333 v 3.0 to follow a nine species non-equilibrium chemistry network (H, H, He, He, He, e, H, H, and H) which includes radiative heating and cooling from these species and metals.444We use a slightly modified version of the the main Grackle repository, available at at changeset c2c0faf. Grackle is a freely available, open source, multi-code library, designed to interface with a wide variety of astrophysical codes. We outline specific model choices made in our simulations and refer the reader to Smith et al. (2017) for a detailed discussion of the code. We apply the Glover & Abel (2008) three-body rate for H formation and include a model for H formation on dust, dust heating, and dust cooling following the methods in Omukai (2000) and Omukai et al. (2005) as included in Grackle. However, we update the default dust to gas ratio scaling in Grackle to account for the steeper scaling in low metallicity regimes (), using the broken power law scalings from Rémy-Ruyer et al. (2014). For metallicities above , this is equivalent to the default behavior of Grackle, where dust content scales linearly with metallicity.

As part of the Grackle package, metal line cooling is modeled using pre-computed Cloudy (Ferland et al., 2013) 555 tables interpolated as a function of density, temperature, and redshift, using the Haardt & Madau (2012) UV metagalactic background. As discussed in more detail in Section 2.5.6, we account for approximate self-shielding of H and He against this UV background. Using this prescription with metal line cooling tables computed under an optically thin assumption can lead to an order of magnitude overestimation of the cooling rate at certain densities, as discussed in Hu et al. (2017) and Appendix D. To address this issue, we use re-computed metal line tables consistent with the self-shielding approximation. We have made these new tables public in the main Grackle repository. These are discussed in greater detail in Appendix  D. Finally, we ignore the effect the stellar radiation field (see Section 2.5.6) may have on the interpolated metal line cooling rates.

2.3 Star Formation Algorithm

In order to resolve individual star formation events on galactic scales, we implement a stochastic star formation algorithm adopted from Goldbaum et al. (2015, 2016). Each cell at the maximum refinement level is capable of forming stars if it contains gas that meets the following local criteria on number density , temperature , cell mass , and velocity : 1) , 2) , 3) , and 4) , where is a resolution dependent density threshold, is a temperature threshold, and is the local thermal Jeans mass. Our fiducial values are and  K. We limit the fraction of a cell’s gas mass that is converted into stars by requiring , where and the maximum star mass . No star formation occurs when , ensuring that a star formation episode can not produce negative densities.

We make the common ansatz that star formation occurs by converting gas into stars in a free fall time with a star formation efficiency, . At high resolution, the choice of should be irrelevant (Orr et al., 2018; Hopkins et al., 2018), as star formation is ultimately self-regulated by feedback.

The stellar mass formed during a timestep from a region with a total gas mass is


In practice, , and is smaller than the minimum star particle mass at parsec scale resolution. We therefore allow star formation to proceed stochastically, following the methods in Goldbaum et al. (2015, 2016), modified for variable stellar masses. In each cell that could potentially form stars, we compute the probability that 100 M of gas will be converted into stars in that time step, and use a random number draw to determine whether or not star formation actually occurs. If it does, we randomly sample from the adopted IMF until approximately 100 M of stars form, keeping the last sampled particle when the total stellar mass formed exceeds 100 M. The total mass of formed star particles is subtracted from the gas mass in the star forming region to ensure mass conservation. We assume a Salpeter IMF (Salpeter, 1955) with , sampling over the range between a minimum stellar mass of 1 M and an arbitrarily chosen maximum stellar mass of 100 M. Our lower limit on stellar masses ensures that we are able to both directly track all particles that contribute in some way to feedback and metal enrichment, and follow longer lived star particles, while reducing the computational expense of following a large number of low mass stars that have no dynamical impact on the galaxy evolution. By ignoring the formation of stars below 1 M, our model in effect spreads this mass into higher mass stars, changing the normalization of the IMF slightly from what would be expected for an IMF that extends below 1 M.

Formed stars are deposited with random positions within the star forming cell and assigned velocities equal to the cell bulk velocity with a 1 km s velocity dispersion. This dispersion captures some of the unresolved gas motions below the resolution limit that are smoothed out by numerical diffusion; it is comparable to, but less than, the velocity dispersion of the coldest gas in our simulations. Stars are assigned metallicities corresponding to the metallicity of the star forming zone, and are chemically tagged with the 17 individual species abundances (H, He, and the 15 metals) that we follow in our simulations.

Stars evolve during the simulation, by losing mass from stellar winds and SNe as described below, and by changing types, but persist throughout the entire simulation. For example, low mass stars are tagged as white dwarfs (WDs) at the end of their life, which may eventually explode as a Type Ia SN (discussed below), after which they persist as massless tracer particle remnants. Finally, each star is marked as a “must refine” particle, requiring that it be surrounded by a four-cell region at the highest level of refinement. This ensures that both stellar winds and SNe feedback are maximally resolved, and that any ejected yields are deposited over a consistent physical scale throughout the simulation.

2.4 Stellar Properties

Given each star’s birth mass and metallicity, we interpolate over the PARSEC grid of stellar evolution tracks (Bressan et al., 2012) to assign a lifetime and AGB phase start time (if any) to it, as well as the effective temperature and surface gravity used in computing radiation properties (see Section 2.5.5). We use the largest subset of the PARSEC models that are regularly sampled in our mass/metallicity space of interest, with 26 mass bins over M and 11 metallicity bins over . Although and evolve over time for stars, modifying stellar radiative properties, following a stellar evolution track for each of our stars is beyond the scope of this work. We instead fix these properties at their zero age main sequence values.

2.5 Stellar Feedback and Chemical Yields

2.5.1 Stellar Yields

For the first time in galaxy-scale simulations, we track galactic chemodynamical evolution using stellar yields ejected from star particles that represent individual stars. We adopt the NuGrid666 collaboration’s set of stellar yields given on a uniformly sampled grid in stellar mass and metallicity with 12 mass bins over M and five metallicity bins at metal fractions of 0.02, 0.01, 0.006, 0.001, and 10 (Pignatari et al., 2016; Ritter et al., 2018). This grid includes yields from the AGB phase of stars from 1–7 M, as well as yields from both stellar winds and core collapse SNe of massive stars from 12–25 M. We complement these tables with tables from Slemer et al. (in prep), based on the PARSEC stellar evolution tracks (Bressan et al., 2012; Tang et al., 2014), to track stellar winds for stars more massive than 25 M. We ignore SN yields from these stars (see next paragraph). We combine all stable isotope yields for a given element into a single elemental abundance for all stable elements from hydrogen to bismuth. Although we can in principle follow an arbitrary number of metal species, practical considerations of memory use prevent this in any given simulation. We refer the reader to previous uses of the NuGrid yields in one-zone galactic chemical enrichment models (Côté et al., 2016b; Côté et al., 2016a; Côté et al., 2017b) for a detailed discussion of how various model uncertainties can influence galactic chemical evolution.

Above some mass within the unsampled range of 7–12 M, stars no longer undergo AGB wind phases but end their lives instead as core collapse SNe. Where this transition occurs is uncertain, but is commonly taken to be at a mass –10 M; we take  M. In our model, stars below this mass eject their wind yields in an AGB phase only at the end of their lives, typically over a period comparable to or less than a few time steps ( kyr). Stars above this mass are assumed to eject their stellar yields via line-driven stellar winds at a constant mass loss rate throughout their lifetime (neglecting Wolf-Rayet and luminous blue variable phases), ending their lives as a core-collapse SN (see Sect. 2.5.2 for details on the wind energetics). Varying changes both the time at which yields for stars around this mass are ejected (for reference, the lifetime of an 8 M star is about 35–40 Myr), and the energy injection from these winds. Côté et al. (2017b) explores how the choice of affects galaxy abundances in a one-zone model. We neglect the effects of binary star evolution on stellar feedback, and discuss the significance of this in Sect. 5.3.

There are large uncertainties in stellar yields for stars more massive than 25 M (see Côté et al., 2016b, and references therein). Indeed, even the exact fate of these stars is uncertain (e.g. Woosley et al., 2002; Zhang et al., 2008; Ugliano et al., 2012), particularly as a function of metallicity (Fryer et al., 2012) with potentially multiple stable and unstable regimes as a function of mass (Heger et al., 2003). Due to this uncertainty, and to avoid erroneously extrapolating from our yield tables, we adopt the simplest model and assume all stars above 25 M end their life through direct collapse to a compact object with no further mass, metal, or energy ejection.

Type Ia SNe are an important additional source of galactic chemical enrichment. These iron group rich events are responsible for the 1 Gyr timescale turnover, or “knee", in vs diagrams. We use the Type Ia SN yields given in Thielemann et al. (1986), adopting a Type Ia SN model as discussed in Sect. 2.5.4. We emphasize that we only track Type Ia SNe occurring within the population of stars formed in this model, neglecting SNe from any possible pre-existing population, substantially limiting the number of Type Ia SNe occurring during the initial gigayear in our models.

2.5.2 Stellar Winds

Stellar winds are important sources of enrichment and feedback in galaxies at both early times from massive stars and late times from AGB stars. Although the energy injected via winds over the lifetime of a cluster of stars is much less than that from SNe and radiation (Shull & Saken, 1995), stellar winds are potentially important sources of pre-SN feedback. We assume complete thermalization of the wind kinetic energy, taking the total injected energy injected in timestep as , where is the thermal energy of the ejected gas mass given the star’s interpolated effective temperature . This mass and energy is injected evenly over a three-cell radius spherical volume centered on the star particle. The edges of this spherical region will only partially overlap grid cells. We use a Monte Carlo sampling method to properly compute the volume of this overlap to scale the injection in these cells appropriately. We assume constant mass loss rates for all winds as set by the yield tables over either the lifetime of the star (for massive stars) or the length of the AGB phase (for low mass stars).

Massive stellar winds have typical velocities of order 10 km s (Leitherer et al., 1992). Satisfying the Courant time step becomes prohibitively expensive following this gas, with time steps dropping to  100 yr. For this reason, we adopt the common simplification of reducing the wind velocity (e.g Offner & Arce, 2015). In our case, we fix massive stellar wind velocities to  km s for stars above 8 M. Our initial tests show that turning off energy injection from stellar winds like this does not significantly affect the global star formation rate of our galaxies. Due to the substantial additional computational expense of following stellar winds for gigayear timescales, we reserve examining the detailed importance of winds to future work. These points are discussed in more detail in Sect. 5.2.1.

Stars that only undergo an AGB phase deposit their feedback at the end of their lives, as determined by the PARSEC evolution tracks. AGB wind velocities vary dramatically over their relatively short lifetimes, but are typically on the order of 10 km s. For simplicity, we adopt a fixed wind velocity of 20 km s for all AGB stars as well.

2.5.3 Core Collapse SNe

Stars between M and 25 M end their lives as core collapse SNe, ejecting mass and metals as determined by the NuGrid stellar yield tables, along with 10 erg of thermal energy. Due to the high resolution of our simulations (1.8 pc), we generally resolve the Sedov phase of each SN explosion well (see Appendix C). We inject thermal energy alone in a three-cell radius spherical region around the star particle, which we find to be sufficient to resolve the SN explosions. We use the same Monte Carlo sampling method as used in our stellar winds to map the spherical injection region to the grid. We continue to track any remaining stellar mass after the SNe occurs as a massive remnant tracer particle. In future work these particles can be used to self-consistently account for more exotic sources of feedback and chemical enrichment such as X-ray binaries and neutron-star binary merger events which, while rare, could be important in long term galaxy evolution (e.g. Artale et al., 2015).

2.5.4 Type Ia SNe

We continue to track low mass stars ( M) after their death as WD particles, marking a subset as Type Ia SN progenitors. This is the most self-consistent model for Type Ia SNe in galaxy-scale simulations. We note however that for the low SFRs in our isolated dwarf galaxy simulation, the first Type Ia SN only appears after a few hundred megayears of simulation time. By the end of the simulation presented here (500 Myr), only 18 have gone off. At the end of their life, we assign a new mass to these particles following the initial-to-final-mass relation of Salaris et al. (2009). We follow the common assumption that progenitor stars with initial masses between 3 and 8 M form WDs that are Type Ia progenitors (see Côté et al., 2017b, and references therein).

We compute the probability that a given Type Ia progenitor will explode as a function of time using an observationally motivated delay time distribution model. The Type Ia SN rate is taken to be a power law in time, , whose slope and normalization are observables. The latter represents the number of Type Ia SNe per mass of star formation. By assuming an IMF , one can write down the fraction of stars capable of forming a Type Ia SN progenitor that will explode within a Hubble time. This is given as


where and are the lower and upper bounds of the IMF, and and are the lower and upper bounds of the range of stars that can form Type Ia candidates. The distribution slope is of order unity, with typical values between 1.0 and 1.2 (see Maoz et al., 2014, for a recent review). can be derived by taking observed values of the Type Ia SN rate and integrating over a Hubble time. Typical values for this are on order of 10 M (Maoz et al., 2014). For our fiducial values, we adopt (Maoz et al., 2010) and M (Graur & Maoz, 2013). Given our choice of IMF, and with  M,  M,  M, and  M, this gives .

Finally, we can normalize to give the probability per unit time that a Type Ia candidate will explode at a time after the formation of its main sequence progenitor. Integrating this gives the total probability at any given time as


where is the formation time of the WD and the leading term on the right hand side properly normalizes the total probability over a Hubble time to . This naturally accounts for both a prompt and delayed Type Ia SN population in our simulations. In practice, rather than drawing a random number for each candidate every timestep, we make a single random number draw, , at the formation time of the white dwarf particle. For , we interpolate its position on a pre-tabulated and inverted cumulative probability distribution function to assign a single time at which the WD particle will explode as a Type 1a supernova. We institute a minimum delay time by defining only for , such that a particle cannot be assigned an explosion time before its formation time.

2.5.5 Ionizing Radiation from Discrete Sources

Radiation feedback, including ionization, ionization heating, and radiation pressure, is an important source of feedback in galaxies. \ionH2 regions carved out by stellar radiation change the ISM structure in regions where SNe eventually explode, generally increasing their dynamical importance. However, accounting for angular effects, radiation can also allow energy from SNe to dissipate more readily by escaping out of channels carved through dense clouds. Radiation feedback effects have been included with various approximations in a wide range of simulations (e.g. Oppenheimer & Davé, 2006; Krumholz et al., 2007; Hopkins et al., 2012; Agertz et al., 2013; Renaud et al., 2013; Stinson et al., 2013; Roškar et al., 2014; Ceverino et al., 2014; Hopkins et al., 2014; Agertz & Kravtsov, 2015; Forbes et al., 2016; Hu et al., 2016, 2017; Hopkins et al., 2018), with a smaller subset using full radiation hydrodynamics (Wise et al., 2012a; Wise et al., 2012b; Wise et al., 2014; Kim et al., 2013a, b; Pawlik et al., 2013; Rosdahl et al., 2015; Aubert et al., 2015; Ocvirk et al., 2016; Baczynski et al., 2015; Pawlik et al., 2017) due to the additional computational expense of direct ray tracing. As we seek a complete accounting of stellar feedback physics, we follow HI and HeI ionizing radiation from our stars through the ray tracing methods described below.

Enzo includes an adaptive ray tracing implementation, Enzo+Moray (Wise & Abel, 2011), to solve the equations of radiative transfer coupled to the hydrodynamics of the simulation. We follow HI and HeI ionizing photons which are coupled to the Grackle primordial chemistry and heating and cooling routines to track photoionization and heating, as well as radiation pressure on hydrogen.

We determine the HI and HeI ionizing photon rates for each star using the OSTAR2002 (Lanz & Hubeny, 2003) grid of O-type stellar models, appropriate for  M at solar metallicity777The exact stellar mass range on the OSTAR2002 grid is model dependent and a function of metallicity. We use linear interpolation in stellar effective temperature, surface gravity, and metallicity to compute the ionizing photon fluxes and rates for each star. Stars less massive than about 15 M and very massive stars with sub-solar metallicity are generally not well sampled by the OSTAR2002 grid. In this case, we integrate a black body spectrum at to obtain the ionizing photon fluxes, but normalize the result to be continuous with the OSTAR2002 grid (see Appendix B).

Instead of assigning a fixed ionizing photon energy across all sources, we integrate over each star’s blackbody curve to compute the average ionizing photon energy individually for each source (see Appendix B). The average energy for HI and HeI ionizing photons changes significantly over the OSTAR2002 temperature range , ranging from 15.72 eV to 20.07 eV and 26.52 eV to 31.97 eV respectively.

We also include the effects of radiation pressure on HI. This has been shown to be important in suppressing the star formation rates of dwarf galaxies by influencing turbulence and the dense gas content of the ISM (Wise et al., 2012a; Ceverino et al., 2014). We ignore the absorption of ionizing radiation by dust and re-radiation in the infrared. This is included in other models (e.g. Rosdahl et al., 2015; Hopkins et al., 2014, 2018) as this may increase by a factor of a few to several the effective radiation pressure (Zhang & Davis, 2017). However, the importance of multiple scattering is still unclear. Other works have shown the effect to only increase the radiation pressure by a factor of order unity (Krumholz & Thompson, 2012, 2013; Reissl et al., 2017; Wibking et al., 2018). Due to these uncertainties, and given that our dwarf galaxy has a low dust content, and therefore a low infrared opacity, we ignore this effect.

2.5.6 Diffuse Heating

We include two forms of diffuse heating in our simulations, each tied directly to the non-equilibrium primordial chemistry network in Grackle: 1) the optically thin, uniform metagalactic UV background (Haardt & Madau, 2012), and 2) localized photoelectric heating from the FUV (6 eV eV) radiation from each of our star particles. The FUV flux for each star is again obtained from the OSTAR2002 grid by directly integrating over the spectral energy distributions for each gridded star. Like the ionizing radiation, we again use an adjusted black body spectrum to compute the flux for stars off of the grid (see Sect. 2.5.5 and Appendix B). Photoelectric heating can be a dominant heating mechanism in the ISM of the Milky Way (Parravano et al., 2003), and could be significant in regulating star formation in dwarf galaxies (Forbes et al., 2016). However, this conclusion warrants further research as its exact importance in dwarf galaxies relative to other feedback mechanisms is contentious (Hu et al., 2016, 2017). Generally, models for photoelectric heating and Lyman-Werner radiation in hydrodynamic simulations of galaxies adopt a constant value or a static, radial profile. Only recently has the localization and time variation of these processes been considered.

Self-shielding of gas against the metagalactic UV background is important in high-resolution simulations, particularly for low-mass, low-metallicity dwarf galaxies where the UV background is capable of gradually photoevaporating unshielded gas from the galaxy (Simpson et al., 2013). We have implemented the Rahmati et al. (2013) approximate self-shielding method in Grackle to account for HI self-shielding against the UV background (see Smith et al., 2017, for more details of this implementation). We assume HeI ionization generally follows HI. This allows us to approximate HeI self-shielding using the same form (A. Rahmati, private communication). We ignore HeII photoionization from the UV background entirely. For consistency, we additionally reduce the reaction rates for direct H ionization (15.4 eV) and H destruction (30 eV) by the same shielding factors computed for HI and HeI shielding.888Ignoring this effect leads to unrealistically high electron fractions in self-shielding gas from direct ionization, which drives significant production of H via gas-phase reactions. Accounting for self-shielding in this manner leads to an inconsistency in using tabulated, optically-thin metal line cooling rates from Cloudy (see Section 4.1.1 of Hu et al. (2017)). As mentioned previously, we have re-computed metal line cooling tables using Cloudy models of optically thick clouds to be consistent with our self-shielding prescription. This is described in more detail in Appendix  D.

We assume the galaxy is mostly optically thin to stellar FUV and use only local approximations for shielding. We calculate the stellar FUV flux in each cell as summed over the contributions from each star to parameterize the local photoelectric heating rate as (Bakes & Tielens, 1994; Wolfire et al., 2003; Bergin et al., 2004)


where is an efficiency factor that depends on , the attenuated local FUV flux


is the dust-to-gas ratio, normalized to the solar value, and is the local FUV flux normalized to the solar neighborhood (Habing, 1968). Aside from a different treatment of and the attenuation, both discussed below, this is equivalent to the method used in Hu et al. (2016, 2017).

The value of is computed consistently with our Grackle dust model, using the broken power law fit from Rémy-Ruyer et al. (2014), as described in Section 2.2. The extremely low dust-to-gas ratio in our modeled galaxies leads to a reduction in the photoelectric heating rate by approximately two orders of magnitude, as compared to a model that assumes a ratio that scales linearly with metallicity at very low metallicity. At these low metallicities, the FUV field only becomes optically thick at length scales of 100 pc for densities of  cm. Given that the ambient density of the ISM is generally 1–10 cm, we can safely assume the FUV field to be optically thin. However, we do include a localized attenuation prescription that may influence high-density or metal-enriched regions of the galaxy. We approximate given in the equation above locally, as , where is the cell width; this approximation is necessarily resolution dependent, but substantially more computationally efficient than direct ray tracing.

Properly computing in Eq. 4 requires an accurate account of the electron number density . This is non-trivial in dense, neutral regions where is dominated by contributions from carbon, dust, and PAH ionizations. Our chemical network only includes contributions from H, He, and H to n. Instead, we use a power-law fit of as a function of from the Wolfire et al. (2003) model of in the solar neighborhood (see Figure 10b of that work); we adopt .

2.5.7 Lyman-Werner Radiation

In addition to the Lyman-Werner radiation from the UV background, we account for localized Lyman-Werner flux from each of our stars to compute the total, local H dissociation rate. We compute the stellar Lyman-Werner flux again from the OSTAR2002 grid by integrating the spectral energy distributions over photon energies from 11.2 eV to 13.6 eV (see Appendix B). Given the local Lyman-Werner flux, the H dissociation rate is taken as , where is the H dissociation cross section. We account for approximate H self-shielding against these sources of Lyman-Werner flux by implementing the Sobolev-like approximation from Wolcott-Green et al. (2011) in Grackle.

3 Galaxy Initial Conditions

We apply these methods to a first test case of the evolution of an isolated, low mass dwarf galaxy. The galaxy is constructed to have initial properties similar to those observed for the Local Group dwarf galaxy Leo P (Giovanelli et al., 2013; McQuinn et al., 2013, 2015a, 2015b), although it is not intended to be a matched model to this galaxy. Leo P is gas rich, with a neutral hydrogen mass  M and stellar mass  M (McQuinn et al., 2015a) extending to an observed neutral hydrogen radius  pc. LeoP has a low metallicity, with an oxygen to hydrogen abundance ratio (O/H) of (Skillman et al., 2013), or a metallicity of (, adopting from Asplund et al. (2009)). Our dwarf galaxy model is constructed without an initial background stellar population, with a total gas mass of  M, of which  M, and , comparable to the average redshift metallicity from the stellar models computed in McQuinn et al. (2015b).

The galaxy initially consists of a smooth, exponential gas disk in hydrostatic equilibrium with a static, background dark matter potential. The gas profile follows Tonnesen & Bryan (2009) and Salem et al. (2015), with


where and are the radial and vertical gas disk scale heights, and is approximately 70% of the total gas mass. We set  pc,  pc, and  M. We adopt a Burkert (1995) dark matter potential with virial mass and radius and  kpc as defined in (Bryan & Norman, 1998) and scale radius r pc. This gives a maximum circular velocity  km s at  kpc. These parameters were adopted specifically to match the observed dynamical mass of Leo P interior to 500 pc of M, and represent virial properties within the halo mass expected for galaxies of this size (Ferrero et al., 2012; Read et al., 2017).

Following the initialization procedure of Hu et al. (2017), we use artificial SN driving to generate realistic initial densities and turbulent properties in the galaxy ISM. This prevents an otherwise uniform collapse of the gas disk at the beginning of the simulation. During this period, SNe explode at a fixed rate of  Myr, corresponding to the SFR obtained given the central HI surface density and the relation presented in (Roychowdhury et al., 2009). We stop the artificial driving 25 Myr after the first star particle forms. These artificial SNe do not drive chemical evolution of the galaxy; their metal yields correspond to the mean ISM abundances. We note this initial driving is ad hoc in that we do not include other effects from the stellar population that would have caused these SNe.

We emphasize here that our model, with no initial stellar distribution, is not intended to identically reproduce the evolution of Leo P, which has formed stars continuously over cosmological timescales. In addition, the mass fractions for the individual metals we track are set to zero so that we follow only the evolution of metals self-consistently produced in the simulations. Otherwise, the galaxy chemical properties would be dominated by the somewhat arbitrary choice of initial abundances. In some ways this is similar to the first pollution of pristine gas in the early Universe, but we note that we cannot directly make this comparison as the environmental conditions and UVB properties were different at high redshift, and we explicitly ignore Pop III and Pop II stellar evolution. Instead, this work is intended as a numerical experiment to investigate how metal enrichment from ongoing star formation proceeds in a gas / dark matter environment similar to low mass halos at both low redshift and the early Universe. The subsequent metal enrichment from the stars in our simulation can be thought of as tracking a change in abundances from arbitrary initial conditions. We will discuss how to make proper comparisons to observed stellar and gas abundance properties of dwarf galaxies in future work, where we will investigate the abundance evolution of our simulations in more detail.

4 Results

We present our initial results here, providing an overview of the morphological (Section 4.1), star formation (Section 4.2), ISM (Section 4.3), radiation field (Section 4.4), outflow (Section 4.5), and chemical (Section 4.6) properties of our dwarf galaxy during the 500 Myr after the first new star forms. Unless otherwise noted, is defined as the time at which that first star particle forms, which is 43 Myr after the actual beginning of the simulation run. The galaxy disk is defined as the fixed physical region within a cylindrical radius of 600 pc and vertical height  pc relative to the center of the galaxy. ISM properties are calculated considering only the gas contained within the disk of the galaxy. We include a resolution study in Appendix E.

Our analysis makes extensive use of the open-source yt toolkit (Turk et al., 2011). All analysis scripts used in this work can be found at at changeset dd76ad10.

4.1 Morphological Structure and Evolution

We begin by characterizing the morphological properties of our dwarf galaxy, as demonstrated in a series of face-on and edge-on images, presented in Figure 1 and Figure 2. These figures show inside-out star formation, as star formation propagates from the inner regions outward during the galaxy’s evolution. This is clear in the face-on panels, which demonstrate the growth of the stellar population from the center outward, and the declining gas densities inside-out as a result of stellar feedback driven winds. This central region quickly fills with warm and hot gas generated from radiation feedback and SNe respectively. Both the ISM and the halo gas are multi-phase, containing gas at cold, warm, and hot temperatures with a range of densities, as evident in the temperature slices in both panels. The ISM properties are quantified further in Section 4.3.

Figure 1: Edge-on views of our dwarf galaxy at four different times in its evolution, 0, 150, 300, and 500 Myr after the beginning of star formation. Shown are the density weighted projection of number density (top row), temperature slices (second row), HI column density (third row), and H column density (fourth row). Each individual main sequence star particle is shown in the number density projections as a single black dot.

The initially puffy gas distribution collapses to a thin disk, with scale heights between 10–30 pc, as shown by the blue line in Fig 3. This figure shows the scale height of all gas in the disk, averaged over 20 Myr periods centered on each given time. Stellar feedback acts to heat up this initially thin disk substantially, creating typical scale heights of 50–120 pc. Towards the end of the simulation time, the half-light radius is pc, where the uncertainty represents the 1 variation in this quantity during the final 20 Myr. Although the disk remains thin beyond the half-light radius, with a scale height of around 50 pc, it is fully resolved at all radii. By the end of the simulation, the majority of the disk has a scale height of 100 pc.

Constraining the gas scale height in ultra-faint dwarf galaxies observationally is challenging. For Leo P, located at 1.7 Mpc, HI observations that are capable of detecting the diffuse HI throughout the galaxy have a resolution of 100–200 pc, with higher resolution observations identifying only the densest HI clumps in the galaxy (e.g. Bernstein-Cooper et al., 2014). In the final column of Figure 2, the peak HI column density reaches  cm, but is located in dense regions with sizes 10pc. With a resolution of 100 pc, the peak column density in an edge-on view is  cm, and  cm in a face-on view. At 500 pc resolution, this corresponds to  cm, and  cm. These column densities are consistent with the resolution-dependent peak column densities found in the low mass dwarf galaxy sample of Teich et al. (2016), and consistent with the observed peak column density of Leo P,  cm, observed with a spatial resolution of about 33 pc.

Figure 2: Same as Figure 1, but showing face-on views.
Figure 3: The total gas scale height at various times throughout the simulation. These times match the images in Figures 1 and  2.

4.2 Star Formation Rate and Mass Evolution

We present the star formation rate (SFR) and core collapse SN rate (SNR) evolution of our dwarf galaxy as measured in 10 Myr bins in the left panel of Figure 4.999We do not show the Type Ia rate as there have only been 16 by the end of the simulation. Within the first 50 Myr of evolution the SFR rises quickly to nearly M yr, declining to  M yr until a significant drop off at about 130 Myr. The remainder of the evolution is characterized by periods of little to no star formation interspersed with periods of continual, but low star formation around  M yr. The SNR tracks the SFR with a time delay, with roughly one core collapse SNe per 100 M of star formation. Averaging over the entire simulation time, we obtain M yr. We discuss how the SFR of this galaxy compares to observed galaxies in Section 5.1.

We note that the granularity in our star formation algorithm creates a lower limit to the SFR that depends on the period over which the SFR is measured. Since we produce stars in  M sets, the smallest value for our measured SFR is . For Myr this is 10 yr. Removing the granularity requires a fundamental change in our star formation algorithm, likely at the cost of increased complexity and computational expense. Sink particles, which track pre-main sequence stellar mass accumulation, would be the most viable way to do this (see for example Krumholz et al., 2004; Federrath et al., 2010; Gong & Ostriker, 2013; Bleuler & Teyssier, 2014; Sormani et al., 2017).

At initialization, all H and He of our dwarf galaxy is neutral, with no molecular hydrogen component. By the time of first star formation ( in Figure 4), HI still dominates the mass of the galaxy, with a molecular hydrogen mass fraction of only  0.3 %. The molecular component declines rapidly as this gas is both converted into stars and is destroyed by stellar radiation feedback. For the remainder of the simulation, the H mass generally increases, with small fluctuations during periods of star formation, reaching a peak mass fraction of 5% at 500 Myr. The growth of the molecular fraction is due in part to a decline in the total gas content of our galaxy from feedback-driven galactic winds. During these outflows, the densest gas, the molecular gas, is preferentially retained over the more diffuse ISM. Examining the molecular properties of the ISM in low mass dwarf galaxies in more detail is a vital avenue of future research, as there are significant observational uncertainties in deriving H content of galaxies in this low metallicity regime (Leroy et al., 2008; McQuinn et al., 2012; Amorín et al., 2016). The molecular properties of our galaxy are discussed further in context with other works in Section 5.1.

Figure 4: Left: The SFR and core collapse SN rate in our dwarf galaxy in 10 Myr bins. Broken portions of this histogram are time periods with no star formation or supernovae. Note that the SN rate has been scaled by a factor of 100 to fit on the same vertical axis as the SFR. Right: The evolution of the total gas mass (black), HI mass (blue), H mass (orange), and stellar mass (red) in the disk of our galaxy over time.

4.3 ISM Properties

Our simulations include sufficient resolution and microphysics to capture a multi-phase medium within the ISM and halo of our simulated galaxy. We define five different gas phases following those defined in (Draine, 2011): molecular gas, cold neutral medium (CNM), warm neutral medium (WNM), warm ionized medium (WIM), and hot ionized medium (HIM). We emphasize that the molecular ISM phase is defined as all cells with M/M f, and is thus somewhat different than simply considering the total H content. By this definition, although our galaxy certainly contains molecular hydrogen, molecular gas as a phase does not exist; the peak f in any single cell remains below 30%. See Appendix A for a quantitative definition of these phases. The properties of these phases are regulated by the complex interplay between cooling, turbulence, self-gravity, and radiative and shock heating from stellar feedback throughout the galaxy’s evolution. Here we discuss the thermodynamic properties of the gas within the inner halo of our dwarf galaxy.

Figure 5 shows the temperature-density distribution of all gas within 0.25 R of the center of the galaxy, averaged over the time period 300–350 Myr. One can readily identify the two regimes containing most of the mass in the simulation: low density, warm gas produced through ionization and SN heating, and cold, high density gas that makes up most of the mass in the galaxy’s disk (see Figure 6). Several notable features of the distribution include: broad ranges of temperature even in quite dense gas, perhaps produced by photoionization and photoelectric heating, a substantial amount of extremely cold gas below 10 K, and the lack of well-defined thermal phases due to the complexity of both the heating and cooling in a turbulent medium. We note that we are likely missing important physics, such as cosmic ray heating and ionization, that would prevent the formation of the coldest gas in this diagram (below about 10 K), but we do not expect this to significantly alter our results. Our artificial temperature ceiling in diffuse gas (see Section 2.1) is seen clearly by the horizontal feature in the top left. The boxed regime in the lower right corner shows our star formation density and temperature threshold. Gas in this regime is rapidly consumed by star formation and subsequent feedback. Given the small size of our dwarf galaxy, the total amount of mass in this regime at any given instant can be small, but does appear in this time average.

Figure 5: The temperature vs. number density phase diagram of our dwarf galaxy simulation showing all gas interior to 0.25 , averaged over a 50 Myr period from t = 300 Myr to t = 350 Myr. The dashed lines are lines of constant pressure, separated by factors of 10. The region in the lower right corner indicates our density and temperature thresholds for star formation.

The mass of the ISM in our dwarf galaxy is dominated by the CNM for the entirety of the simulation, as shown in the left panel of Figure 6. The mass fraction of the remaining phases are ordered by temperature, with the WNM as the next most-significant component. The WNM is initially comparable to the CNM, but comprises a mass fraction of about 0.1 by the end of the run. The WIM and HIM fluctuate significantly, corresponding to fluctuations in the SFR and associated feedback, but are subdominant throughout the simulation. During periods of peak stellar feedback, however, the WIM can reach a mass fraction above 0.1. Although the CNM dominates the mass fraction, it is a negligible component of the ISM volume, which is WIM dominated. However, the large, anti-correlated fluctuations in the WNM and HIM make these three phases often comparable. Together, these figures better quantify the general properties observed in the panel plots in Figure 1 and Figure 2.

These results are in contrast with those found for the more massive dwarf galaxy modeled by Hu et al. (2016, 2017). They find the mass and volume fraction of the ISM are nearly entirely dominated by warm gas (defined in those works as gas with 100 KT), with cold gas having between 1 and 10% of the mass, and occupying negligible volume. Hot gas (defined as gas with  K) occupies 10% of the volume, with negligible mass, in their galaxy, while our WIM alone occupies 50% of the volume. Our lower mass, lower metallicity galaxy contains more cold gas (by mass fraction) and hot gas (by volume fraction) that seen in the more massive dwarf galaxy in these works. The driver of these differences, which are likely somewhat related to differences in the dark matter halo potential, will be investigated in future work. We have compared our cooling curves to those used in (Hu et al., 2017) and found them to be comparable; though this could contribute to the differences, it is likely not the dominant source.

Figure 6: The evolution of the mass and volume fractions for each phase of the model galaxy’s ISM. See Appendix A for definitions of each phase.

4.4 Interstellar Radiation Field

The interstellar radiation field (ISRF) of our galaxy varies dramatically in both space and time, as has been seen previously in works modeling varying radiation fields both as expected from stellar motions in our own galaxy (Parravano et al., 2003), and in models including radiation (e.g. Hu et al., 2017). This is not particularly surprising in our low SFR regime, where there can be large fluctuations over time as individual massive stars form, move about, and evolve. In Figure 7 we present azimuthally averaged radial profiles of the ISRF in various bands, time averaged over 100 Myr during the period of star formation from roughly 250–350 Myr. The top panel shows , the ISRF flux between 6–13.6 eV normalized to the value in the solar neighborhood of the Milky Way (see Sec. 2.5.6). The contribution from the metagalactic background falls below the lower vertical axis limit on this plot. The averaged profile varies between values of 0.02 and 0.1, with peaks located at radii of the few active star formation regions. At any given radius, there is over a two order magnitude variation in the ISRF during this period of time.

The bottom panel gives the HI ionizing photon flux from stellar radiation. The ionizing radiation profile follows a similar trend, yet with significantly more variation, anywhere from two to four orders of magnitude. As this radiation is followed through radiative transfer, the profile encodes information about local attenuation by dense, neutral gas. This is the main driver of the differences between the two panels. The total fluctuation in both panels is due in part to the low-level, stochastic star formation in our galaxy. A higher star formation rate would produce a more regular population of massive stars and more uniform (in time) ISRF.

Figure 7: Azimuthally averaged radial profiles of the ISRF in the mid-plane of our galaxy in two different bands, time-averaged over 50 Myr from 300 – 350 Myr. Here we define the midplane as within 2.5  of z = 0, or 4.5 pc. The top panel gives , the flux of radiation between 6–13.6 eV normalized to the value in the solar neighborhood, shaded between minimum and maximum values at each position, with the average shown as a black line. The bottom panel gives the HI ionizing stellar radiation flux. Since this radiation is tracked directly through radiative transfer, the minimum value at all radii is 0 at some point. For this reason we only shade between the first quartile and maximum values. HeI ionizing radiation is very similar to HI ionizing radiation, with a small vertical offset, and is not shown for clarity.
Figure 8: Single-snapshot 2D radial profile plots at 300 Myr of the ISRF in two flux bands, and HI ionizing radiation, illustrating the full dynamic range of radiation flux at a given radius in the galaxy. Here, we include all gas within the mid-plane of our dwarf galaxy. Since a majority of the mass of the galaxy is in the cold phase (see Figure 5), and is therefore optically thick to HI ionizing radiation, it does not show up in the HI ionizing radiation diagram. This gas readily appears in the diagram since we assume it to be optically thin, though we do apply a localized shielding approximation.

To further quantify the local variations in these radiation fields, we present the full distribution of and the HI ionizing flux in Figure 8 at a single snapshot at 300 Myr. This diagram shows how dramatic the increase in ISRF near young, massive stars is (the spikes in both diagrams), while much of the mid-plane sees an ISRF orders of magnitude lower. The striking contrast between the two diagrams is due to the shielding of the HI ionizing flux in the most massive (cold and dense) regions of the galaxy through the radiative transfer calculations; shielding of the FUV radiation is approximate and in general weaker, making these regions more prominent in the left hand figure (the pink/white clumps). From both of these diagrams, it is clear that the ISRF of a low mass dwarf galaxy varies greatly over time and space in a way that cannot be appropriately captured by an analytic profile. Although one could adopt an averaged radial profile to provide a realistic, global source of energy for thermal pressure support of the gas against collapse, it is unclear how sufficient this would be in suppressing star formation. In particular, the large increases around sites of recent star formation could be important sources of feedback to destroy molecular clouds and reduce their effective star formation efficiency. It remains to be seen which of these two modes of feedback is more important in regulating star formation.

4.5 Outflow Properties

The recent FIRE cosmological simulations of dwarf galaxies over a range of dark matter halo masses find that they exhibit large outflows, with mass loading factors () on order of 100–1000 (Muratov et al., 2015). However, comparable models of idealized dwarf galaxies with detailed feedback and physics treatments find more modest mass loading factors (Hu et al., 2016, 2017). In Figure 9 we present the mass outflow and mass loading rates for our dwarf galaxy as a function of time, computed at five different positions from the galaxy. We follow Muratov et al. (2015) in defining the mass outflow rate at any given radius to be the sum of the outflow rate in all cells in a spherical annulus of width centered at that radius,


We choose , or 2.74 kpc.

The total mass outflow rates and mass loading factors at 0.1, 0.25, 0.5, and 1.0 are shown in Figure 9. Generally, other works use gigayear timescale measurements of the SFR to compute the mass loading factor. For consistency with those works, we use the 500 Myr average SFR for computing the mass loading factor. The outflow rate at 0.1 is high, corresponding to mass loading factors between 20–100 throughout the simulation time. This declines towards larger radii, however, with substantially less outflow past the virial radius. Muratov et al. (2015) finds typical mass loading factors at 0.25 R on order of 20–40 for galaxies with km s at low redshift, consistent with our results. The fluctuations in both of these panels are directly correlated with the SFR, with increased outflow during periods of star formation, and decreased outflow during periods of quiescence.

Interestingly, the v 30 km s halos examined in Muratov et al. (2015) are more massive than the M M halo examined here by a factor of a few. Using a fit provided in Muratov et al. (2015) to extrapolate and compare at fixed halo mass, one would expect mass loading factors on order of 100 at 0.25 R for our dwarf galaxy, a factor of a few higher than what we find. These differences could be attributed to our lack of cosmological evolution in these isolated simulations, but ultimately requires a larger set of dwarf galaxy simulations to make a more robust comparison. We note, however, that our results are closer to the Muratov et al. (2015) results than those in Hu et al. (2016, 2017), which find lower mass loading factors even closer to the disk, at 0.05 , between 1 and 10 for a dwarf galaxy with M M; certainly this implies even smaller mass loading factors at 0.25 .

Figure 9: Spherical mass outflow rates (Eq. [7]) and mass loading rates over time at 4 different radii from the galaxy.

Detailed outflow properties, beyond outflow rates and mass loading factors, can help discriminate between the model dependent feedback physics included in galaxy simulations. In Figure 10 we present radial velocity distributions of all material outside our dwarf galaxy’s disk, and within the halo, broken into three gas phases. Gas with a negative velocity is moving towards the center of the halo. Roughly 25% of this mass is inflowing, mostly with modest negative velocities, and corresponds to previously ejected gas mixing and recycling throughout the halo. Half of the outflowing gas (positive velocities) is moving at velocities below 30 km s, 75% at velocities below 70 km s, and 95% at velocities below 100 km s. Although the mass contained in the tails of these distributions is a sub-dominant fraction of the total, there is still a non-negligible amount of gas moving at velocities of a few hundreds of km/s, with a peak velocity of over 700 km s. The WNM and WIM together dominate the mass of both the inflowing and outflowing gas, with the WIM and HIM dominating at velocities above 200 km s. The dominant launching mechanism in this simulation is SN feedback, which generates a rapidly moving and volume-filling WIM and HIM, consistent with the results in Hu et al. (2016, 2017). However, as shown, the HIM, which is mostly the SN ejecta itself, comprises very little of the outflow by mass. Most of the outflowing gas (by mass) comes from the warm phase, pushed out by the high pressure, fast moving HIM. Some of this warm gas certainly originates from adiabatically and radiatively cooled HIM, however. The amount of transfer between phases in the halo of our galaxy will be investigated in future work.

Figure 10: The time averaged radial velocity distribution of outflowing material external to the disk and within the virial radius of our dwarf galaxy. This is averaged over the same time interval as Figure 7. The outflowing material is multiphase, broken into WNM, WIM, and HIM. See Section 4.3 for definitions of these regimes. We note that WNM is often labeled simply as “cold”, or gas below  K. There is little to no outflowing mass in the CNM.

4.6 Metal Evolution

4.6.1 Metal Enriched Outflows

Dwarf galaxies efficiently, and preferentially, eject metals released in stellar feedback from their shallow potential wells (Mac Low & Ferrara, 1999; Ferrara & Tolstoy, 2000). This has been better quantified recently both observationally (e.g. Kirby et al., 2011; Zahid et al., 2012; Peeples et al., 2014; McQuinn et al., 2015b) and with more detailed cosmological simulations (Simpson et al., 2013; Anglés-Alcázar et al., 2017; Muratov et al., 2017). In the top panel of Fig. 11, we give the metal mass loading factor for our galaxy over time, at the same radii as in Fig. 9. The parameter used to quantify metal outflow efficiencies varies among works. Here, we define the metal mass loading factor as the metal outflow rate divided by the metal production rate, or


where is the metal mass outflow rate, is the total mass in metals produced, and is the total mass in stars. These metal loading factors fluctuate significantly with the SFR, just as was shown in Figure 9, reaching a minimum of about 0.05, but peaking at around 5. On average, over the simulation time, is below unity (around 0.5). Recent simulations of outflows from a Milky Way type disk indicate typical comparable to our results, usually between 0.5 and 1 (Li et al., 2017; Fielding et al., 2017). Muratov et al. (2017) computes a slightly different quantity for their galaxies, the normalized metal outflow rate , finding values of about 0.02 at 0.25  regardless of galaxy circular velocity. Our galaxy is consistent with this value, with an average , fluctuating between 0.007 and 0.02.

Figure 11: Top: Metal mass loading factor (see Eq.[ 8]) at the same radii as in Figure 9. This is the ratio between the metal outflow rate and the metal production rate. Bottom: The fraction of metals contained in the disk, CGM, and outside the halo of our dwarf galaxy over time. In both panels, we consider the total mass of all individually tracked metal species, which is zero at initialization, not the aggregate total metallicity field, which is non-zero at initialization.

These large metal mass loading factors indicate that a majority of the metals produced in our dwarf galaxy are ejected from the disk. This is quantified in the lower panel of Fig. 11, where we show the mass fraction of metals in the disk, circumgalactic medium (CGM), and outside the virial radius of our galaxy over time. After the first 20 Myr, SN-driven winds rapidly drive out large quantities of metals from the disk and into the galaxy’s halo. This continues throughout the simulation, with only 4% of produced metals residing in the disk of the galaxy. It only takes 150 Myr for some metals to reach the virial radius of the halo, with a steadily increasing fraction continually leaving the virial radius until about 350 Myr where the fraction levels off to just above 50%. Likewise, the CGM metal content continually decreases until the end of the simulation from loss through the virial radius of the halo. See Section 5.1.3 for further discussion.

4.6.2 Differential Evolution of Elements Within the ISM

It is important to understand how metals from each source of stellar yields enrich the ISM. Observations of more massive dwarf galaxies than those simulated here indicate fairly uniform radial gas-phase metallicity profiles, even beyond the stellar radius (e.g. Werk et al., 2011; Belfiore et al., 2017). This requires that metal mixing and transport occur on hundred megayear timescales, much more rapidly than the gigayear timescale expected from assuming transport at the cold gas sound speed. Therefore, either metals are transported first through a hot phase with high sound speed, or through efficient turbulent mixing within the ISM (e.g. de Avillez & Mac Low, 2002; Tassis et al., 2008; Yang & Krumholz, 2012). It remains uncertain how metal abundances vary in detail within these galaxies, beyond one-dimensional radial profiles, and whether or not abundance distributions depend on the metal species. It is even more unclear how metals are transported and distributed within low-mass dwarf galaxies, which generally host too few \ionH2 regions for a detailed examination.

We demonstrate the power of our simulations, which capture a realistic ISM at high resolution with multiple feedback sources, by addressing these questions in Fig 12. The left panel gives the abundance ratio of N to O throughout the ISM. The right two panels give the slices of number density (top right) and temperature (bottom right) in the mid-plane of a portion of our dwarf galaxy. These show regions with dense, cold gas clouds (, K) connected by cold filaments, warm, diffuse gas ( cm,  K), and hot gas from a recent SN explosion ( K).

Figure 12: Three slices in the mid-plane of our dwarf galaxy at 300 Myr after the start of star formation showing the variation in gas phase metal abundances. The left slice gives the ratio of the abundance between N and O, normalized to the solar abundance, while the number density and temperature are shown on the right. In each, we mark massive stars with active stellar winds as white points and SNe and AGB-phase enrichment events that occurred in the preceding 5 Myr as black stars and orange diamonds respectively.

As shown in the left panel, [N/O] varies significantly over this section of the ISM with notable differences across the various phases and ISM structures. The hottest gas, dominated by recent supernova explosions is overabundant in oxygen, relative to solar (purple). However, the relative abundance of nitrogen increases in the WIM and WNM, being overabundant relative to solar (green).

At our adopted metallicity, nitrogen is predominantly produced in AGB star winds, with very little production in core collapse SNe and winds from more massive stars. Therefore, nitrogen is injected into the ISM with significantly less energy ( km s) than elements produced in SNe, like oxygen, ( km s). Given the variations in Figure 12, the energetic differences between injection sources can drive abundance variations within the ISM of our dwarf galaxy. The regions most rich in N are sites of recent AGB winds that have yet to mix with the rest of the ISM. This suggests that metal mixing within the ISM (and also metal ejection from the ISM) is species dependent. A more detailed analysis is beyond the scope of this work, but we investigate this in detail in Emerick et. al. (in prep.).

5 Discussion

5.1 Comparison to Observed Low Mass Dwarf Galaxies

As noted in Section 3, our galaxy model is not intended to directly reproduce the observed properties of Local Group ultra-faint dwarfs. Notably, our initial conditions neglect a pre-existing stellar population and are only followed for 500 Myr, a fraction of the age of dwarf galaxies. However, we can still place our model in context with observations using simple comparisons to the star formation rate (Section 5.1.1), molecular gas (Section 5.1.2), and metal retention fraction (Section 5.1.3) properties of observed dwarf galaxies. We show that these properties are broadly consistent with observations.

5.1.1 Gas and Star Formation

The observational sample of isolated, gaseous, low mass dwarf galaxies is limited compared to more massive galaxies, but has improved substantially with recent blind and targeted HI surveys (e.g. Giovanelli et al., 2005; Geha et al., 2006; Geha et al., 2012; Walter et al., 2008; Cannon et al., 2011; Haynes et al., 2011; Hunter et al., 2012; Bradford et al., 2015; James et al., 2015; Tollerud et al., 2015; Sand et al., 2015; Wang et al., 2017). However, the sample of isolated, gaseous dwarf galaxies with  M remains small. In Figure  13 we show where our galaxy lies relative to the observed Kennicutt-Schmidt relation and extended Schmidt law for low mass galaxies. In both diagrams, our simulations are given by the colored points, sampled every megayear throughout the entire simulation.

Although simple to measure in simulations, these quantities are challenging to directly compare to observations. We have attempted to make a reasonable analog to how and are measured observationally for low mass dwarfs Roychowdhury et al. (see 2014). We define , where is the SFR measured over the preceding 10 Myr, and is the area of the disk within the radius of the outermost star formed within the previous 10 Myr. Likewise, , where M is the total gas mass within this defined disk. However, the total gas content cannot be determined observationally. To match this limitation, we follow Roychowdhury et al. (2014) and take , where the factor 1.34 attempts to account for He. We note that there is generally no correction made for any possible H or HII content. As shown in Section 5.1.2, these components may be significant.

Figure 13: The Kennicutt-Schmidt law (left) and extended Schmidt law (right) relationships for our galaxy as measured every megayear, plotted as points colored by time, with dark / purple early and light / green late. See text for details of the calculation. Recent observations from the SHIELD sample (Teich et al., 2016) are plotted as black points with error bars. On the left, we also give the best fit line to galaxies from the FIGGS sample from Roychowdhury et al. (2014), and on the right we also show the best fit line and 1 errors from Shi et al. (2011). There is no clear correlation with time in this diagram.

We include recent observational constraints on these relationships in Figure  13. Our galaxy fluctuates significantly about both relationships with no clear trends in time. However, in both cases, it is consistent with the available observational sample. At times our galaxy exhibits gas surface densities below the observational constraints. The trend is still consistent with higher densities at this point, but with a larger scatter towards lower star formation rate densities and efficiencies.

In constructing our galaxy model, we employed no tuning of the underlying physics, adopting only canonical values for any available free parameters. It is thus non-trivial that our galaxy should oscillate about the median relationships in Fig. 13, and signifies a proper accounting of the relevant physics governing gas and star forming properties in our galaxy. This result is consistent with galactic evolution simulations run at high resolution with a detailed accounting of stellar feedback physics (see Naab & Ostriker, 2017, and references therein) and with the demonstration by Li et al. (2005) that the Kennicutt-Schmidt law can be reproduced by gravitation acting on an isothermal disk. The net result of the complex interactions of heating, cooling , chemistry, and feedback physics on star formation is to offset to a level not too dissimilar to more simple simulations considering gravity alone.

5.1.2 Molecular Gas Content

The molecular gas content of low-mass dwarf galaxies is generally assumed to be small, but is not well constrained by either theory or observations. Assuming the relationship in Leroy et al. (2013) and Momose et al. (2013), Roychowdhury et al. (2014) finds typical molecular gas mass fractions can be anywhere from to ; a significant range in values. As discussed in Section 2.4, our galaxy has of 0.001 – 0.05, just overlapping with this range. H formation in our model is possible through formation on dust, the three body interaction, or the gas-phase reaction H + H  H + e. The gas-phase reaction dominates in our low-metallicity galaxy over the other two channels by several orders of magnitude. In our model, H is produced solely through the reaction H + e H + . Thus the presence of some ionizing background is required to generate the molecular fractions we find in our simulations, as confirmed in test simulations run only to the onset of first star formation.

In contrast, Hu et al. (2016) and Hu et al. (2017) find low molecular fractions () even in their simulations without any feedback (). However, although these works do contain a non-equilibrium chemical model, they do not include either H or a background radiation field. Our model suggests that both H and the UV background are critical components in H formation in small dwarf galaxies. Indeed it is not surprising that the gas phase reactions dominate over grain catalysis at low metallicity (Glover, 2003), though it is worth noting that the rate coefficients associated with gas-phase H formation are still uncertain by an order of magnitude (Glover et al., 2006; Glover & Mac Low, 2007). Our model does lack additional chemistry that may be important to the formation and destruction of H, including HD chemistry, C and O chemistry, a detailed dust model, and cosmic rays. It is unclear how the combination of all of these effects will behave, especially considering that even at this resolution, we are unable to resolve the high density turbulent density perturbations in which H forms most efficiently (Glover & Mac Low, 2007). This uncertainty certainly warrants further study of molecular gas in low metallicity dwarf galaxies.

5.1.3 Metal Retention

The simulations presented here have not yet been run for the gigayear timescales required to begin to make direct comparisons to the observed stellar and gas phase metallicities of comparable dwarf galaxies at . However, we can compare to a key observable: the retention fraction of metals within stars and the galaxy’s ISM compared to what would be expected from closed box stellar evolution models given the galaxy’s star formation history. This can be done readily with Milky Way dSph’s. Their stars seem to retain very little of the expected metal production: on order of a few percent or less depending on the galaxy and the species (Kirby et al., 2011). However, environmental effects, namely ram pressure and tidal stripping, complicate the understanding of how these metals were removed from the galaxy.

Leo P, the dwarf galaxy we approximate in our initial conditions, is extremely valuable as a gas-rich, star forming, low-mass dwarf galaxy, with an observable HII region, necessary for determining gas phase abundances, that is close enough to the Milky Way to conduct this experiment. Leo P retains 5% 2% of its metals, 1% in stars and the rest in ISM gas (McQuinn et al., 2015b). As discussed in Section 4.6, more than 90% of the tracked metals produced during our simulation no longer reside within the galaxy’s disk, agreeing with observations. However, this is an evolving quantity that also depends on how much (if any) subsequent re-accretion of these metals occurs. Although more than half of these metals are expected to eventually re-accrete (Anglés-Alcázar et al., 2017), it is still possible that most of the metals that have been produced at these early times will remain outside the galaxy disk.

It is interesting to consider whether ejected metals in our model reside in the CGM or have been ejected into the intergalactic medium. While this cannot be observationally determined for dwarf galaxies at this mass, cosmological simulations of an M galaxy show that, by redshift zero 40% is ejected from the galaxy’s halo (Anglés-Alcázar et al., 2017). We find 5.3% of all metals lie within 1 kpc of the center of our galaxy, 24.5% within 0.25 (or 6.85 kpc), and 51% outside . This is consistent with previous works, but this is again an evolving quantity. In addition, the amount of gas that escapes the virial radius is certainly sensitive to the details of gas accretion from the IGM on this galaxy, which we cannot capture in this model.

The re-accretion or final ejection of this gas is directly relevant to the chemical evolution of low mass dwarf galaxies. Recycling of metal enriched gas could be a significant driver of long-term chemical evolution in low mass galaxies, particularly if a majority of metals ejected from the disk (itself nearly all the metals produced by the galaxy) return. In addition, the accretion of pristine gas from the intergalactic medium could significantly affect the gas flows around the galaxy, possibly promoting the retention of ejected metals. This effect is not included in our isolated galaxy simulations, and its role is beyond the scope of this work.

5.2 Missing Physics

Although we include many detailed physical models in our simulations, there remain additional physical processes that may be relevant, which we now discuss.

5.2.1 Massive Stellar Wind Energy

Our massive stellar wind model drastically reduces the injected wind velocity from 1000 km s to 20 km s. Although our algorithm is entirely capable of generating realistic stellar winds with velocities comparable to those observed, such fast winds place a near constant and severe constraint on the Courant time step that renders 100 Myr simulations impractical. When considered in isolation, stellar winds are an important source of pre-SN feedback and can dramatically influence dynamical evolution on molecular cloud and galaxy scales (Dale & Bonnell, 2008; Peters et al., 2017; Gatto et al., 2017). However, when considered together with ionizing radiation, stellar winds contain less total energy (Agertz et al., 2013) and do not seem to have a significant dynamical influence in either idealized simulations (Geen et al., 2015) or individual giant molecular clouds (Dale et al., 2014), unless densities near the ionizing source are high enough to trap the HII region in the source cell. In that case, they can clear out a cavity to allow initial establishment of the HII region.

They are even less relevant in the low-metallicity regime studied here, as stellar winds become weaker with decreasing metallicity (Puls et al., 2000; Vink & de Koter, 2005). Although they likely have minimal dynamical importance at resolutions where peak densities are not anyway high enough to trap ionization fronts, a full model of stellar winds may affect detailed ISM properties and metal mixing, warranting closer examination in future work.

5.2.2 Cosmic Rays

Recent work has explored the importance of cosmic ray feedback in regulating the ISM and wind properties in galactic disks (Hanasz et al., 2013; Girichidis et al., 2016; Simpson et al., 2016; Farber et al., 2017), isolated galaxies (Salem et al., 2016; Salem et al., 2015; Pakmor et al., 2016; Ruszkowski et al., 2017), and galaxies in cosmological context (Salem et al., 2014). These relativistic charged particles act as a source of non-thermal pressure support in the galaxy’s ISM, capable of driving outflows at different velocities and containing different thermal phases than those driven through thermal feedback alone (Salem et al., 2016). Modeling cosmic rays is challenging, however, as they encompass a wide range of energies, and there are significant uncertainties in how they propagate through the ISM (e.g. Wiener et al., 2017). Their propagation is often modeled as a diffusive process, but in truth this diffusion should vary depending on cosmic ray energy. In addition, cosmic rays couple effectively to the magnetic fields of galaxies, diffusing preferentially along structured magnetic field lines within the ISM. Modeling cosmic ray feedback completely thus requires both an accurate cosmic ray model and magnetohydrodynamics in order to capture the interplay of these two physical phenomena. Finally, including MHD presents additional difficulties in untangling the effects of each individual feedback mechanism on galaxy chemodynamics.

We do note that an isotropic, two-fluid model for cosmic ray feedback exists in Enzo (Salem & Bryan, 2014; Salem et al., 2015) and has been well tested. Mechanically, including this relatively simple treatment of cosmic ray feedback in our model is trivial. However, the cosmic ray population, their diffusion coefficient, and the magnetic field structure of the lowest mass dwarf galaxies each have significant enough uncertainties to warrant reserving their full inclusion into our model to later work.

5.3 Detailed Stellar Evolution and Binary Stars

Roughly half of massive stars live in binary pairs (Sana et al., 2013). Their interactions, primarily through mass transfer, can significantly alter their radiation properties and lifetimes. This can change both how much and how long these stars emit ionizing radiation, an important source of stellar feedback, and where and when these stars explode as SNe. This effect could be significant, but is rarely accounted for in galaxy evolution models, which are commonly based on calculations of single star evolution (e.g. STARBURST99). For example, Zapartas et al. (2017) finds that binarity extends the timescales over which core collapse SNe occur from a given star formation event, from a maximum time of 50 Myr to  200 Myr. Although they find only of core collapse SNe explode after 50 Myr, this could still be an important effect. Properly accounting for the delay times due to variations in individual star lifetimes has already been shown to change the significance of feedback and influence galaxy metallicity properties (Kimm et al., 2015). Extending the lifetimes of these stars could increase the role of radiation feedback, however this may be less important as these additional photons are more likely to escape the galaxy (e.g. Ma et al., 2016b).

Since we model stars on a star-by-star basis, both of these effects could be reasonably accounted for by stochastically assigning binary star properties to some subset of our individual stars. This is beyond the scope of this project, but will be investigated in future work.

6 Conclusion

We have developed a new method for simulating galaxy evolution with detailed feedback and chemical enrichment. For the first time on galaxy scales, we simultaneously model multi-channel stellar feedback in detail, using individual star particles to model core collapse and Type Ia SNe, ionizing radiation followed through radiative transfer, photoelectric heating, Lyman-Werner radiation and pollution from AGB and massive stellar winds. This treatment of feedback, coupled with the detailed chemistry and heating/cooling physics followed with Grackle, allows us to capture realistic galaxy evolution in detail. In this work, we apply these methods to simulate the evolution of an isolated, low-mass, dwarf galaxy modeled after the properties of the Leo P dwarf galaxy. We present an overview of the properties of this simulation in this work.

For our simulated dwarf galaxy, we find:

  1. Multi-channel feedback is effective in regulating star formation to a rate consistent with the Kennicutt-Schmidt relationship and the extended Schmidt law in observed galaxies. (See Figures 4 and 13).

  2. This feedback drives large outflows having mass loading factors of at 0.25 R, falling to at R, and metal mass loading factors near unity. By mass, nearly all of this outflow is moving with velocities below 100 km s, but there is a significant tail towards velocities up to 1000 km s (See Figure 10).

  3. Only 4% of metals are retained in the disk of our simulated galaxy, consistent with the observed metal retention fractions of low-mass dwarfs. By the end of the simulation 45% of the remaining metals stay within the virial radius (but outside the galaxy), while 50% have been ejected beyond the virial radius. (See Figures 9, 10, and 11.)

  4. Beyond the stellar radius, the gas scale height is thin ( pc), yet resolved, with larger scale heights (–200 pc) driven by feedback interior to the stellar radius. This is comparable to the resolution limit of the diffuse HI in observed, gaseous low-mass dwarfs. At a spatial resolution of 100 pc, our galaxy has a peak HI column density  cm, depending on inclination (See Figure 3).

  5. The ISRF of our galaxy varies strongly in both space and time by orders of magnitude. It is unclear how important these fluctuations are as a source of feedback, or if the affect can be approximated with a time-averaged radial profile. The importance of radiation feedback in our model will be investigated in a future work. (See Figures 7 and 8.)

  6. We find H fractions below 5% in our dwarf galaxy, consistent with the poor constraints on molecular gas formation in low metallicity dwarf galaxies. This H forms entirely through gas-phase reactions facilitated by H in self-shielding regions; H formation on dust grains and in the three body reaction are both insignificant. Cold, neutral hydrogen dominates the mass of our galaxy. While warm, neutral hydrogen is present, it does not dominate the mass fraction (See Figure 4 and Figure 6.)

  7. Finally, we present gas-phase oxygen and nitrogen distributions as examples to briefly demonstrate that there are marked differences in how individual metal species are distributed within the ISM of our galaxy. These variations could be tied to differences in elemental yields among different sources (for example, AGB winds vs. SNe), as suggested by Krumholz & Ting (2018). This will be explored in more detail in future work (Emerick et. al. in prep.).


We would like to thank the following for their advice and valuable discussions, without which this work would not have been possible: B. Côté, S. Glover, K. Hawkins, C. Hu, K. Johnston, B. O’Shea, M. Putman, B. Smith, J. Wall, and J. Wise. We would additionally like to thank the referee for their careful report. A.E. is funded by the NSF Graduate Research Fellowship DGE 16-44869. G.L.B. is funded by NSF AST-1312888, NASA NNX15AB20G, and NSF AST-1615955. M.-M.M.L. was partly funded by NASA grant NNX14AP27G. We gratefully recognize computational resources provided by NSF XSEDE through grant number TGMCA99S024, the NASA High-End Computing Program through the NASA Advanced Supercomputing Division at Ames Research Center, Columbia University, and the Flatiron Institute. This work made significant use of many open source software packages, including yt, Enzo, Grackle, Python, IPython, NumPy, SciPy, Matplotlib, HDF5, h5py, Astropy, Cloudy and deepdish. These are products of collaborative effort by many independent developers from numerous institutions around the world. Their commitment to open science has helped make this work possible.


  • Agertz & Kravtsov (2015) Agertz O., Kravtsov A. V., 2015, ApJ, 804, 18
  • Agertz et al. (2013) Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, ApJ, 770, 25
  • Amorín et al. (2016) Amorín R., Muñoz-Tuñón C., Aguerri J. A. L., Planesas P., 2016, A&A, 588, A23
  • Anglés-Alcázar et al. (2017) Anglés-Alcázar D., Faucher-Giguère C.-A., Kereš D., Hopkins P. F., Quataert E., Murray N., 2017, MNRAS, 470, 4698
  • Artale et al. (2015) Artale M. C., Tissera P. B., Pellizza L. J., 2015, MNRAS, 448, 3071
  • Asplund et al. (2009) Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481
  • Aubert et al. (2015) Aubert D., Deparis N., Ocvirk P., 2015, MNRAS, 454, 1012
  • Baczynski et al. (2015) Baczynski C., Glover S. C. O., Klessen R. S., 2015, MNRAS, 454, 380
  • Bakes & Tielens (1994) Bakes E. L. O., Tielens A. G. G. M., 1994, ApJ, 427, 822
  • Belfiore et al. (2017) Belfiore F., et al., 2017, MNRAS, 469, 151
  • Bergin et al. (2004) Bergin E. A., Hartmann L. W., Raymond J. C., Ballesteros-Paredes J., 2004, ApJ, 612, 921
  • Bernstein-Cooper et al. (2014) Bernstein-Cooper E. Z., et al., 2014, AJ, 148, 35
  • Bleuler & Teyssier (2014) Bleuler A., Teyssier R., 2014, MNRAS, 445, 4015
  • Bradford et al. (2015) Bradford J. D., Geha M. C., Blanton M. R., 2015, ApJ, 809, 146
  • Bressan et al. (2012) Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127
  • Brook et al. (2014) Brook C. B., Stinson G., Gibson B. K., Shen S., Macciò A. V., Obreja A., Wadsley J., Quinn T., 2014, MNRAS, 443, 3809
  • Bryan & Norman (1998) Bryan G. L., Norman M. L., 1998, ApJ, 495, 80
  • Bryan et al. (1995) Bryan G. L., Norman M. L., Stone J. M., Cen R., Ostriker J. P., 1995, Computer Physics Communications, 89, 149
  • Burkert (1995) Burkert A., 1995, ApJ, 447, L25
  • Cannon et al. (2011) Cannon J. M., et al., 2011, ApJ, 739, L22
  • Ceverino et al. (2014) Ceverino D., Klypin A., Klimek E. S., Trujillo-Gomez S., Churchill C. W., Primack J., Dekel A., 2014, MNRAS, 442, 1545
  • Cohen et al. (2013) Cohen J. G., Christlieb N., Thompson I., McWilliam A., Shectman S., Reimers D., Wisotzki L., Kirby E., 2013, ApJ, 778, 56
  • Colella & Woodward (1984) Colella P., Woodward P. R., 1984, Journal of Computational Physics, 54, 174
  • Corlies et al. (2018) Corlies L., Johnston K. V., Wise J. H., 2018, MNRAS, 475, 4868
  • Côté et al. (2016a) Côté B., West C., Heger A., Ritter C., O’Shea B. W., Herwig F., Travaglio C., Bisterzo S., 2016a, MNRAS, 463, 3755
  • Côté et al. (2016b) Côté B., Ritter C., O’Shea B. W., Herwig F., Pignatari M., Jones S., Fryer C. L., 2016b, ApJ, 824, 82
  • Côté et al. (2017a) Côté B., Silvia D., O’Shea B. W., Smith B., Wise J. H., 2017a, ApJ, submitted,
  • Côté et al. (2017b) Côté B., O’Shea B. W., Ritter C., Herwig F., Venn K. A., 2017b, ApJ, 835, 128
  • Dale & Bonnell (2008) Dale J. E., Bonnell I. A., 2008, MNRAS, 391, 2
  • Dale et al. (2014) Dale J. E., Ngoumou J., Ercolano B., Bonnell I. A., 2014, MNRAS, 442, 694
  • Davé et al. (2017) Davé R., Rafieferantsoa M. H., Thompson R. J., Hopkins P. F., 2017, MNRAS, 467, 115
  • De Silva et al. (2015) De Silva G. M., et al., 2015, MNRAS, 449, 2604
  • Draine (2011) Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium
  • Escala et al. (2018) Escala I., et al., 2018, MNRAS, 474, 2194
  • Farber et al. (2017) Farber R., Ruszkowski M., Yang H.-Y. K., Zweibel E. G., 2017, ApJ, in press,
  • Federrath et al. (2010) Federrath C., Banerjee R., Clark P. C., Klessen R. S., 2010, ApJ, 713, 269
  • Ferland et al. (2013) Ferland G. J., et al., 2013, Rev. Mex. Astron. Astrofis., 49, 137
  • Ferrara & Tolstoy (2000) Ferrara A., Tolstoy E., 2000, MNRAS, 313, 291
  • Ferrero et al. (2012) Ferrero I., Abadi M. G., Navarro J. F., Sales L. V., Gurovich S., 2012, MNRAS, 425, 2817
  • Few et al. (2012) Few C. G., Courty S., Gibson B. K., Kawata D., Calura F., Teyssier R., 2012, MNRAS, 424, L11
  • Few et al. (2014) Few C. G., Courty S., Gibson B. K., Michel-Dansac L., Calura F., 2014, MNRAS, 444, 3845
  • Fielding et al. (2017) Fielding D., Quataert E., Martizzi D., Faucher-Giguère C.-A., 2017, MNRAS, 470, L39
  • Forbes et al. (2016) Forbes J. C., Krumholz M. R., Goldbaum N. J., Dekel A., 2016, Nature, 535, 523
  • Frebel & Norris (2015) Frebel A., Norris J. E., 2015, ARA&A, 53, 631
  • Fryer et al. (2012) Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, ApJ, 749, 91
  • Gatto et al. (2017) Gatto A., et al., 2017, MNRAS, 466, 1903
  • Geen et al. (2015) Geen S., Rosdahl J., Blaizot J., Devriendt J., Slyz A., 2015, MNRAS, 448, 3248
  • Geha et al. (2006) Geha M., Blanton M. R., Masjedi M., West A. A., 2006, ApJ, 653, 240
  • Geha et al. (2012) Geha M., Blanton M. R., Yan R., Tinker J. L., 2012, ApJ, 757, 85
  • Gilmore et al. (2012) Gilmore G., et al., 2012, The Messenger, 147, 25
  • Giovanelli et al. (2005) Giovanelli R., et al., 2005, AJ, 130, 2598
  • Giovanelli et al. (2013) Giovanelli R., et al., 2013, AJ, 146, 15
  • Girichidis et al. (2016) Girichidis P., et al., 2016, ApJ, 816, L19
  • Glover (2003) Glover S. C. O., 2003, ApJ, 584, 331
  • Glover & Abel (2008) Glover S. C. O., Abel T., 2008, MNRAS, 388, 1627
  • Glover & Mac Low (2007) Glover S. C. O., Mac Low M.-M., 2007, ApJ, 659, 1317
  • Glover et al. (2006) Glover S. C., Savin D. W., Jappsen A.-K., 2006, ApJ, 640, 553
  • Goldbaum et al. (2015) Goldbaum N. J., Krumholz M. R., Forbes J. C., 2015, ApJ, 814, 131
  • Goldbaum et al. (2016) Goldbaum N. J., Krumholz M. R., Forbes J. C., 2016, ApJ, 827, 28
  • Gong & Ostriker (2013) Gong H., Ostriker E. C., 2013, ApJS, 204, 8
  • Graur & Maoz (2013) Graur O., Maoz D., 2013, MNRAS, 430, 1746
  • Haardt & Madau (2012) Haardt F., Madau P., 2012, ApJ, 746, 125
  • Habing (1968) Habing H. J., 1968, Bull. Astron. Inst. Netherlands, 19, 421
  • Hanasz et al. (2013) Hanasz M., Lesch H., Naab T., Gawryszczak A., Kowalik K., Wóltański D., 2013, ApJ, 777, L38
  • Haynes et al. (2011) Haynes M. P., et al., 2011, AJ, 142, 170
  • Heger et al. (2003) Heger A., Fryer C. L., Woosley S. E., Langer N., Hartmann D. H., 2003, ApJ, 591, 288
  • Hirai & Saitoh (2017) Hirai Y., Saitoh T. R., 2017, ApJ, 838, L23
  • Hopkins et al. (2012) Hopkins P. F., Quataert E., Murray N., 2012, MNRAS, 421, 3488
  • Hopkins et al. (2014) Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS, 445, 581
  • Hopkins et al. (2018) Hopkins P. F., et al., 2018, MNRAS,
  • Hu et al. (2016) Hu C.-Y., Naab T., Walch S., Glover S. C. O., Clark P. C., 2016, MNRAS, 458, 3528
  • Hu et al. (2017) Hu C.-Y., Naab T., Glover S. C. O., Walch S., Clark P. C., 2017, MNRAS, 471, 2151
  • Hunter et al. (2012) Hunter D. A., et al., 2012, AJ, 144, 134
  • James et al. (2015) James B. L., Koposov S., Stark D. P., Belokurov V., Pettini M., Olszewski E. W., 2015, MNRAS, 448, 2687
  • Jeon et al. (2017) Jeon M., Besla G., Bromm V., 2017, ApJ, 848, 85
  • Kim et al. (2013a) Kim J.-h., Krumholz M. R., Wise J. H., Turk M. J., Goldbaum N. J., Abel T., 2013a, ApJ, 775, 109
  • Kim et al. (2013b) Kim J.-h., Krumholz M. R., Wise J. H., Turk M. J., Goldbaum N. J., Abel T., 2013b, ApJ, 779, 8
  • Kimm et al. (2015) Kimm T., Cen R., Devriendt J., Dubois Y., Slyz A., 2015, MNRAS, 451, 2900
  • Kirby et al. (2011) Kirby E. N., Martin C. L., Finlator K., 2011, ApJ, 742, L25
  • Kobayashi et al. (2000) Kobayashi C., Tsujimoto T., Nomoto K., 2000, ApJ, 539, 26
  • Kroupa et al. (2013) Kroupa P., Weidner C., Pflamm-Altenburg J., Thies I., Dabringhausen J., Marks M., Maschberger T., 2013, The Stellar and Sub-Stellar Initial Mass Function of Simple and Composite Populations. p. 115, doi:10.1007/978-94-007-5612-0_4
  • Krumholz & Thompson (2012) Krumholz M. R., Thompson T. A., 2012, ApJ, 760, 155
  • Krumholz & Thompson (2013) Krumholz M. R., Thompson T. A., 2013, MNRAS, 434, 2329
  • Krumholz & Ting (2018) Krumholz M. R., Ting Y.-S., 2018, MNRAS, 475, 2236
  • Krumholz et al. (2004) Krumholz M. R., McKee C. F., Klein R. I., 2004, ApJ, 611, 399
  • Krumholz et al. (2007) Krumholz M. R., Klein R. I., McKee C. F., 2007, ApJ, 656, 959
  • Lanz & Hubeny (2003) Lanz T., Hubeny I., 2003, ApJS, 146, 417
  • Leitherer et al. (1992) Leitherer C., Robert C., Drissen L., 1992, ApJ, 401, 596
  • Leroy et al. (2008) Leroy A. K., Walter F., Brinks E., Bigiel F., de Blok W. J. G., Madore B., Thornley M. D., 2008, AJ, 136, 2782
  • Leroy et al. (2013) Leroy A. K., et al., 2013, AJ, 146, 19
  • Li et al. (2005) Li Y., Mac Low M.-M., Klessen R. S., 2005, ApJ, 620, L19
  • Li et al. (2017) Li M., Bryan G. L., Ostriker J. P., 2017, ApJ, 841, 101
  • Ma et al. (2016a) Ma X., Hopkins P. F., Faucher-Giguère C.-A., Zolman N., Muratov A. L., Kereš D., Quataert E., 2016a, MNRAS, 456, 2140
  • Ma et al. (2016b) Ma X., Hopkins P. F., Kasen D., Quataert E., Faucher-Giguère C.-A., Kereš D., Murray N., Strom A., 2016b, MNRAS, 459, 3614
  • Mac Low & Ferrara (1999) Mac Low M.-M., Ferrara A., 1999, ApJ, 513, 142
  • Machacek et al. (2001) Machacek M. E., Bryan G. L., Abel T., 2001, ApJ, 548, 509
  • Majewski et al. (2010) Majewski S. R., Wilson J. C., Hearty F., Schiavon R. R., Skrutskie M. F., 2010, in Cunha K., Spite M., Barbuy B., eds, IAU Symposium Vol. 265, Chemical Abundances in the Universe: Connecting First Stars to Planets. pp 480–481, doi:10.1017/S1743921310001298
  • Majewski et al. (2016) Majewski S. R., APOGEE Team APOGEE-2 Team 2016, Astronomische Nachrichten, 337, 863
  • Maoz et al. (2010) Maoz D., Sharon K., Gal-Yam A., 2010, ApJ, 722, 1879
  • Maoz et al. (2014) Maoz D., Mannucci F., Nelemans G., 2014, ARA&A, 52, 107
  • Marcolini et al. (2008) Marcolini A., D’Ercole A., Battaglia G., Gibson B. K., 2008, MNRAS, 386, 2173
  • McQuinn et al. (2012) McQuinn K. B. W., Skillman E. D., Dalcanton J. J., Dolphin A. E., Cannon J. M., Holtzman J., Weisz D. R., Williams B. F., 2012, ApJ, 751, 127
  • McQuinn et al. (2013) McQuinn K. B. W., et al., 2013, AJ, 146, 145
  • McQuinn et al. (2015a) McQuinn K. B. W., et al., 2015a, ApJ, 812, 158
  • McQuinn et al. (2015b) McQuinn K. B. W., et al., 2015b, ApJ, 815, L17
  • Momose et al. (2013) Momose R., et al., 2013, ApJ, 772, L13
  • Muratov et al. (2015) Muratov A. L., Kereš D., Faucher-Giguère C.-A., Hopkins P. F., Quataert E., Murray N., 2015, MNRAS, 454, 2691
  • Muratov et al. (2017) Muratov A. L., et al., 2017, MNRAS, 468, 4170
  • Naab & Ostriker (2017) Naab T., Ostriker J. P., 2017, ARA&A, 55, 59
  • Obreja et al. (2014) Obreja A., Brook C. B., Stinson G., Domínguez-Tenreiro R., Gibson B. K., Silva L., Granato G. L., 2014, MNRAS, 442, 1794
  • Ocvirk et al. (2016) Ocvirk P., et al., 2016, MNRAS, 463, 1462
  • Offner & Arce (2015) Offner S. S. R., Arce H. G., 2015, ApJ, 811, 146
  • Omukai (2000) Omukai K., 2000, ApJ, 534, 809
  • Omukai et al. (2005) Omukai K., Tsuribe T., Schneider R., Ferrara A., 2005, ApJ, 626, 627
  • Oppenheimer & Davé (2006) Oppenheimer B. D., Davé R., 2006, MNRAS, 373, 1265
  • Oppenheimer & Davé (2008) Oppenheimer B. D., Davé R., 2008, MNRAS, 387, 577
  • Orr et al. (2018) Orr M. E., et al., 2018, MNRAS, 478, 3653
  • Pakmor et al. (2016) Pakmor R., Pfrommer C., Simpson C. M., Springel V., 2016, ApJ, 824, L30
  • Pan et al. (2013) Pan L., Scannapieco E., Scalo J., 2013, ApJ, 775, 111
  • Parravano et al. (2003) Parravano A., Hollenbach D. J., McKee C. F., 2003, ApJ, 584, 797
  • Pawlik et al. (2013) Pawlik A. H., Milosavljević M., Bromm V., 2013, ApJ, 767, 59
  • Pawlik et al. (2017) Pawlik A. H., Rahmati A., Schaye J., Jeon M., Dalla Vecchia C., 2017, MNRAS, 466, 960
  • Peeples et al. (2014) Peeples M. S., Werk J. K., Tumlinson J., Oppenheimer B. D., Prochaska J. X., Katz N., Weinberg D. H., 2014, ApJ, 786, 54
  • Peters et al. (2017) Peters T., et al., 2017, MNRAS, 466, 3293
  • Pflamm-Altenburg & Kroupa (2006) Pflamm-Altenburg J., Kroupa P., 2006, MNRAS, 373, 295
  • Pignatari et al. (2016) Pignatari M., et al., 2016, ApJS, 225, 24
  • Puls et al. (2000) Puls J., Springmann U., Lennon M., 2000, A&AS, 141, 23
  • Rahmati et al. (2013) Rahmati A., Pawlik A. H., Raičevic̀ M., Schaye J., 2013, MNRAS, 430, 2427
  • Read et al. (2017) Read J. I., Iorio G., Agertz O., Fraternali F., 2017, MNRAS, 467, 2019
  • Reissl et al. (2017) Reissl S., Klessen R. S., Mac Low M.-M., Pellegrini E. W., 2017, A&A, p. submitted
  • Rémy-Ruyer et al. (2014) Rémy-Ruyer A., et al., 2014, A&A, 563, A31
  • Renaud et al. (2013) Renaud F., et al., 2013, MNRAS, 436, 1836
  • Revaz & Jablonka (2012) Revaz Y., Jablonka P., 2012, A&A, 538, A82
  • Revaz et al. (2009) Revaz Y., et al., 2009, A&A, 501, 189
  • Revaz et al. (2016) Revaz Y., Arnaudon A., Nichols M., Bonvin V., Jablonka P., 2016, A&A, 588, A21
  • Ritter et al. (2018) Ritter C., Herwig F., Jones S., Pignatari M., Fryer C., Hirschi R., 2018, MNRAS,
  • Robertson & Kravtsov (2008) Robertson B. E., Kravtsov A. V., 2008, ApJ, 680, 1083
  • Roederer et al. (2014) Roederer I. U., Cowan J. J., Preston G. W., Shectman S. A., Sneden C., Thompson I. B., 2014, MNRAS, 445, 2970
  • Rosdahl et al. (2015) Rosdahl J., Schaye J., Teyssier R., Agertz O., 2015, MNRAS, 451, 34
  • Roškar et al. (2014) Roškar R., Teyssier R., Agertz O., Wetzstein M., Moore B., 2014, MNRAS, 444, 2837
  • Roychowdhury et al. (2009) Roychowdhury S., Chengalur J. N., Begum A., Karachentsev I. D., 2009, MNRAS, 397, 1435
  • Roychowdhury et al. (2014) Roychowdhury S., Chengalur J. N., Kaisin S. S., Karachentsev I. D., 2014, MNRAS, 445, 1392
  • Ruszkowski et al. (2017) Ruszkowski M., Yang H.-Y. K., Zweibel E., 2017, ApJ, 834, 208
  • Salaris et al. (2009) Salaris M., Serenelli A., Weiss A., Miller Bertolami M., 2009, ApJ, 692, 1013
  • Salem & Bryan (2014) Salem M., Bryan G. L., 2014, MNRAS, 437, 3312
  • Salem et al. (2014) Salem M., Bryan G. L., Hummels C., 2014, ApJ, 797, L18
  • Salem et al. (2015) Salem M., Besla G., Bryan G., Putman M., van der Marel R. P., Tonnesen S., 2015, ApJ, 815, 77
  • Salem et al. (2016) Salem M., Bryan G. L., Corlies L., 2016, MNRAS, 456, 582
  • Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
  • Sana et al. (2013) Sana H., et al., 2013, A&A, 550, A107
  • Sand et al. (2015) Sand D. J., et al., 2015, ApJ, 806, 95
  • Sarmento et al. (2017) Sarmento R., Scannapieco E., Pan L., 2017, ApJ, 834, 23
  • Sawala et al. (2010) Sawala T., Scannapieco C., Maio U., White S., 2010, MNRAS, 402, 1599
  • Schaye et al. (2010) Schaye J., et al., 2010, MNRAS, 402, 1536
  • Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  • Shen et al. (2010) Shen S., Wadsley J., Stinson G., 2010, MNRAS, 407, 1581
  • Shen et al. (2013) Shen S., Madau P., Guedes J., Mayer L., Prochaska J. X., Wadsley J., 2013, ApJ, 765, 89
  • Shi et al. (2011) Shi Y., Helou G., Yan L., Armus L., Wu Y., Papovich C., Stierwalt S., 2011, ApJ, 733, 87
  • Shull & Saken (1995) Shull J. M., Saken J. M., 1995, ApJ, 444, 663
  • Simpson et al. (2013) Simpson C. M., Bryan G. L., Johnston K. V., Smith B. D., Mac Low M.-M., Sharma S., Tumlinson J., 2013, MNRAS, 432, 1989
  • Simpson et al. (2016) Simpson C. M., Pakmor R., Marinacci F., Pfrommer C., Springel V., Glover S. C. O., Clark P. C., Smith R. J., 2016, ApJ, 827, L29
  • Skillman et al. (2013) Skillman E. D., et al., 2013, AJ, 146, 3
  • Smith et al. (2008) Smith B., Sigurdsson S., Abel T., 2008, MNRAS, 385, 1443
  • Smith et al. (2017) Smith B. D., et al., 2017, MNRAS, 466, 2217
  • Snaith et al. (2015) Snaith O., Haywood M., Di Matteo P., Lehnert M. D., Combes F., Katz D., Gómez A., 2015, A&A, 578, A87
  • Sormani et al. (2017) Sormani M. C., Treß R. G., Klessen R. S., Glover S. C. O., 2017, MNRAS, 466, 407
  • Stinson et al. (2013) Stinson G. S., Brook C., Macciò A. V., Wadsley J., Quinn T. R., Couchman H. M. P., 2013, MNRAS, 428, 129
  • Su et al. (2017a) Su K.-Y., et al., 2017a, preprint, (arXiv:1712.02795)
  • Su et al. (2017b) Su K.-Y., Hopkins P. F., Hayward C. C., Faucher-Giguère C.-A., Kereš D., Ma X., Robles V. H., 2017b, MNRAS, 471, 144
  • Tang et al. (2014) Tang J., Bressan A., Rosenfield P., Slemer A., Marigo P., Girardi L., Bianchi L., 2014, MNRAS, 445, 4287
  • Tassis et al. (2008) Tassis K., Kravtsov A. V., Gnedin N. Y., 2008, ApJ, 672, 888
  • Teich et al. (2016) Teich Y. G., et al., 2016, ApJ, 832, 85
  • Thielemann et al. (1986) Thielemann F.-K., Nomoto K., Yokoi K., 1986, A&A, 158, 17
  • Tollerud et al. (2015) Tollerud E. J., Geha M. C., Grcevich J., Putman M. E., Stern D., 2015, ApJ, 798, L21
  • Tonnesen & Bryan (2009) Tonnesen S., Bryan G. L., 2009, ApJ, 694, 789
  • Torrey et al. (2017) Torrey P., et al., 2017, ApJ, submitted,
  • Truelove et al. (1997) Truelove J. K., Klein R. I., McKee C. F., Holliman II J. H., Howell L. H., Greenough J. A., 1997, ApJ, 489, L179
  • Turk et al. (2011) Turk M. J., Smith B. D., Oishi J. S., Skory S., Skillman S. W., Abel T., Norman M. L., 2011, ApJS, 192, 9
  • Ugliano et al. (2012) Ugliano M., Janka H.-T., Marek A., Arcones A., 2012, ApJ, 757, 69
  • Verner et al. (1996) Verner D. A., Ferland G. J., Korista K. T., Yakovlev D. G., 1996, ApJ, 465, 487
  • Vink & de Koter (2005) Vink J. S., de Koter A., 2005, A&A, 442, 587
  • Vorobyov et al. (2015) Vorobyov E. I., Recchi S., Hensler G., 2015, A&A, 579, A9
  • Walter et al. (2008) Walter F., Brinks E., de Blok W. J. G., Bigiel F., Kennicutt Jr. R. C., Thornley M. D., Leroy A., 2008, AJ, 136, 2563
  • Wang et al. (2017) Wang J., et al., 2017, MNRAS, 472, 3029
  • Weidner & Kroupa (2004) Weidner C., Kroupa P., 2004, MNRAS, 348, 187
  • Werk et al. (2011) Werk J. K., Putman M. E., Meurer G. R., Santiago-Figueroa N., 2011, ApJ, 735, 71
  • Wibking et al. (2018) Wibking B. D., Thompson T. A., Krumholz M. R., 2018, MNRAS, submitted,
  • Wiener et al. (2017) Wiener J., Pfrommer C., Peng Oh S., 2017, MNRAS, 467, 906
  • Wiersma et al. (2009) Wiersma R. P. C., Schaye J., Theuns T., Dalla Vecchia C., Tornatore L., 2009, MNRAS, 399, 574
  • Wise & Abel (2011) Wise J. H., Abel T., 2011, MNRAS, 414, 3458
  • Wise et al. (2012a) Wise J. H., Abel T., Turk M. J., Norman M. L., Smith B. D., 2012a, MNRAS, 427, 311
  • Wise et al. (2012b) Wise J. H., Turk M. J., Norman M. L., Abel T., 2012b, ApJ, 745, 50
  • Wise et al. (2014) Wise J. H., Demchenko V. G., Halicek M. T., Norman M. L., Turk M. J., Abel T., Smith B. D., 2014, MNRAS, 442, 2560
  • Wolcott-Green et al. (2011) Wolcott-Green J., Haiman Z., Bryan G. L., 2011, MNRAS, 418, 838
  • Wolfire et al. (2003) Wolfire M. G., McKee C. F., Hollenbach D., Tielens A. G. G. M., 2003, ApJ, 587, 278
  • Woosley et al. (2002) Woosley S. E., Heger A., Weaver T. A., 2002, Reviews of Modern Physics, 74, 1015
  • Yang & Krumholz (2012) Yang C.-C., Krumholz M., 2012, ApJ, 758, 48
  • Zahid et al. (2012) Zahid H. J., Dima G. I., Kewley L. J., Erb D. K., Davé R., 2012, ApJ, 757, 54
  • Zapartas et al. (2017) Zapartas E., et al., 2017, A&A, 601, A29
  • Zhang & Davis (2017) Zhang D., Davis S. W., 2017, ApJ, 839, 54
  • Zhang et al. (2008) Zhang W., Woosley S. E., Heger A., 2008, ApJ, 679, 639
  • de Avillez & Mac Low (2002) de Avillez M. A., Mac Low M.-M., 2002, ApJ, 581, 1047

Appendix A Gas Phases of the ISM

To aid in comparison to other works, we define our ISM phases here, as adopted and modified from Draine (2011), Table 1.3, and used consistently throughout this analysis. By construction, these phases are mutually exclusive. We take to be the molecular hydrogen fraction of the total gas mass.

  1. Hot Ionized Medium (HIM): K

  2. Warm Ionized Medium (WIM):

  3. Warm Neutral Medium (WNM):

  4. Cold Neutral Medium (CNM): K,

  5. Molecular:  K,

Appendix B Stellar Radiation Properties

We determine the radiation properties of our model stars as a function of stellar mass and metallicity from the OSTAR2002 grid (Lanz & Hubeny, 2003) in the regime that it covers, or integrated from an adjusted black body curve for other stars, given the ZAMS stellar radius obtained from the PARSEC (Bressan et al., 2012; Tang et al., 2014) stellar evolution data set. In Figure B.1 we plot these properties. We use a constant factor across metallicities for stars not sampled on the OSTAR2002 grid to shift the black body radiation fluxes to be roughly continuous with the ionizing photon rates and luminosities as a function of stellar mass. This requires two factors for each radiation type, one for low mass and one for high mass stars. We use the following multiplicative factors to adjust the black body spectrum for HI and HeI ionizing radiation respectively: [0.1, 3.2] and [0.0001, 4.0]. We do not find this adjustment to be necessary for the FUV and Lyman-Werner radiation bands.

The ionizing radiation photon energies are taken as the average ionizing photon energy for a black body of the star’s given T obtained from the PARSEC grid. A more accurate approach to compute this energy would convolve the full stellar spectrum and the frequency-dependent absorption cross section. We tested this approach using the frequency-dependent photoionization cross sections from Verner et al. (1996).101010Source code containing the analytic fits given in Verner et al. (1996) was obtained from The blackbody approximation is accurate to within 5%, yet substantially easier to compute on the fly, as integrals over the blackbody spectrum can be expressed as infinite series that rapidly converge to high precision. Unlike for the ionizing radiation, we assume constant FUV and Lyman-Werner band energies for each star, at 9.8 eV and 12.8 eV respectively.

Figure B.1: Radiation properties for our model stars, showing the ionizing photon luminosities for HI (top left) and HeI (top right) for each star, the FUV (middle left) and Lyman-Werner (middle right) luminosities, and finally the average ionizing photon energy for HI (bottom left) and HeI (bottom right). Note, we only track radiation from stars above 8.0 M, which dominate over less massive stars, even when accounting for IMF weighting.

Appendix C Typical Gas Densities in Supernova Injection Regions

Modeling SN feedback with the injection of thermal energy alone can lead to rapid overcooling of the injected energy, and a significant underestimate of the effects of SN feedback. Often, ad hoc solutions to this problem are used in large simulations with coarse resolution, such as momentarily turning off cooling in feedback affected regions or decoupling affected regions from the hydrodynamics for some time. However, physically consistent solutions have been developed (e.g Simpson et al., 2016) that inject a mixture of kinetic and thermal energy with a ratio that depends on resolution and local gas density. Overcooling becomes less important with higher resolution and lower ISM densities, until eventually pure thermal energy injection is sufficient to resolve each SN. We take advantage of our high resolution and employ a simple thermal energy injection model for our SN feedback. We demonstrate how well this resolves our SNe in the comparison simulation presented in this work.

The left panel of Figure C.1 gives the distributions of the peak and average ISM number densities in the SN injection regions for each SN in our fiducial simulation. As shown, a majority explode in regions at substantially lower densities than the star formation threshold of 200 cm. For most SNe,  cm, which is due to the substantial pre-SN feedback included in our simulations. The right panel of Figure  C.1 gives the fractional distribution of the calculated radius of the pressure-driven snowplow () phase of the Sedov-Taylor expansion, adopting the definitions from Simpson et al. (2016)


where is the injection energy in units of 10 erg, is the number density of the medium, and is the metallicity in units of . Simpson et al. (2016) found that is needed to resolve the SN explosion with thermal energy injection alone, assuming uniform density in the injection region. In practice, the injection region is never uniform. We give assuming uniform density at both the average and maximum number density in the injection region of each SN. As shown, a majority of SNe are resolved (to the right of the 4.5 line), but up to 8.3% are not well resolved and 2.0% are completely unresolved (below a single cell width). In general, these are SN explosions from the most massive (i.e. most prompt) stars. We do not expect that resolving these SNe will dramatically alter our results. However, we note that using our feedback model at resolutions much less than a few parsecs will require implementing an different injection mechanism.

Figure C.1: Left: Distribution of peak (orange) and average (black) gas densities in the injection region of each of the SNe from the simulation. Right: Distribution of the radius of the pressure-driven snowplow phase for each SN assuming a uniform medium at either the peak or average density in the injection region. The vertical lines show one and 4.5 cell widths. This shows that a fraction of these SNe are certainly unresolved, with R less than a single cell size, and somewhat more don’t satisfy the 4.5 cell criterion, but the majority are well resolved.

Appendix D Cooling and Heating Rates

Figure D.1: The absolute value of the net cooling or heating rate in two models that account for self-shielding of primordial gas against a metagalactic UV background. Our model (black) uses self-consistent metal line cooling rates, and is compared to an incorrect model (orange) that adopts the uncorrected (i.e. optically thin) metal line cooling rates. The regimes where heating dominates are plotted with dashed lines for clarity.

The cooling curves used in this work include a correction to a significant conceptual inconsistency one may encounter when using tabulated metal line cooling rates in combination with a self-shielding approximations. Hu et al. (2017) details this inconsistency: metal cooling rates computed under the assumption of an optically thin UV background overestimate the electron fraction in regimes where self-shielding is important (). This results in an artificially enhanced metal line cooling rate in regions of self-shielding. This can be a significant effect, particularly at higher metallicities. We demonstrate this in Figure D.1, which gives the absolute value of the net cooling rate for our full model at  Z (discussed below) as compared to a model using the optically thin metal line cooling rates at a range of densities. The effect is most significant at moderate densities, where the net cooling rate can be an order of magnitude higher. In low temperature regions, this can significantly shift the inferred equilibrium temperature and reduce the cooling time of dense gas. Figure D.2 shows the net cooling rate and the individual heating and cooling rates from our generated tables at low metallicity,  Z, comparable to that used in our dwarf galaxy. We have made these tables publicly available in the main distribution of Grackle and discuss how they were generated below.

For densities where self-shielding is important, we computed Cloudy models of one-dimensional clouds at each temperature and density pair with a physical size corresponding to the Jeans length. In cases where the Jeans length is very large, we limit the size to 100 pc, a reasonable approximation for the maximum size of a self-gravitating cloud. We then adopt the metal line cooling rates obtained from the center of these clouds, whose outside is exposed to the UV background. These models were computed with a modified version of the CIAOLoop111111 (our version: code used in Smith et al. (2008). In computing these tables we had to include the cosmic ray ionization that dominates in optically thick clouds. Without this effect, Cloudy becomes unstable in entirely neutral regions. We adopted a cosmic ray ionization rate of 10% the Milky Way value, though we note that varying this from 1% to 100% the Milky Way value did not have an effect on the extracted metal line cooling rates.

Finally, we note that the cooling and heating rates in our simulation will deviate from these curves, as they depend upon local conditions that cannot be accounted for here. The primordial cooling and heating rates are computed consistently with the non-equilibrium chemistry solver, leading to deviations from the tabulated primordial rates, and the heating rates depend upon the local stellar radiation field, which is not included in these diagrams.

Figure D.2: The total heating (dotted) and cooling (dashed) rates extracted from the core of a Jeans-length sized cloud irradiated by a Haardt & Madau (2012) metagalactic UV background as modeled in Cloudy. The absolute value of the net cooling or heating rate is shown with solid lines.

Appendix E Resolution Study

We demonstrate the effect of varying resolution in the evolution of our galaxy model in Figure E.1. We compare the star formation rate, mass evolution, and metal ejection / retention fractions between our fiducial high-resolution simulation (solid lines) with the same simulation at two lower maximum resolutions, 3.6 pc (dashed lines), and 7.2 pc (dash-dotted lines). The physics in each simulation remains the same with the exception of star formation, which employs a resolution-dependent density threshold. A factor of two decrease in resolution translates to a factor of four decrease in the star formation density threshold, from 200 cm in our fiducial runs, to 50 cm and 12.5 cm in our 3.6 pc and 7.2 pc resolution runs respectively. As a result, the time at which star formation first occurs varies between the simulations, but we again define t = 0 as the time of first star formation in each simulation.

Figure E.1: A comparison of the SFR (left), disk gas mass (middle), and metal retention (right) evolution for our fiducial run (solid), 3.6 pc run (dotted), 3.6 pc run with a doubled initial supernova driving rate (dashed), and 7.2 pc run (dash-dotted). The line colors in the left panel are for better clarity between the SFRs.

Figure E.1 shows that, while decreasing the resolution certainly changes the evolution, the simulations are similar to within a factor of a few. We argue that these differences are within the stochastic variation one would expect from simulations with stochastic star formation. The star formation rate and total stellar mass is correlated with resolution, with our fiducial run having the highest SFR and stellar mass among the three runs. We note that the self-consistent SNe are still resolved in the 3.6 pc runs since these still explode at the low densities seen in Figure C.1. However, although feedback in general is still effective in the 7.2 pc run, SNe are not well resolved at this resolution.

Generally ineffective feedback leads to increased, not decreased, star formation rates. However, the lower resolution run has the lowest SFR. The source of this difference is likely related to the ability to resolve dense, cold gas in the lower resolution simulations. In particular, the low resolution simulation clearly is unable to resolve the densities at which gas-phase H formation becomes efficient, leading to significantly lower H fractions (middle panel) and lower cooling rates. In addition, whatever H does form at these lower resolutions exists at lower gas densities than in our fiducial run, and is less well protected by self-shielding effects. This is shown clearly in the significant plummet in the H mass during the phase of first star formation that is less severe in our fiducial run. In this low metallicity regime, where H is an important coolant, the ability to resolve the densities responsible for H formation and self-shielding is important.

Theses differences in densities are shown clearly by comparing Figure 5 with Figure E.2. Although the lowest temperature reached in each simulation is comparable, there is less gas in the coldest, densest portion of the phase diagram in each of the lower resolution simulations. In the 7.2 pc run, the gas phases are notably less distinct, with gas more evenly smoothed out in the multi-phase region between cold, dense gas and warm, ionized gas. We note the 7.2 pc simulation, which ran quickly, was conducted without the temperature ceiling used in the 3.6 pc and fiducial runs; we do not expect this to have a significant effect on the comparison.

Figure E.2: Temperature-density phase diagrams for our 3.6 pc run (left) and 7.2 pc run (right), as presented for our fiducial run in Figure 5.

Finally, the metal retention fraction (right panel) between the fiducial and 3.6 pc runs are quite similar, but is significantly higher for the lowest resolution simulation. It is likely, however, that this is a result of the lower SFR in the 7.2 pc run and is not directly a resolution effect. This also affects the fraction of metals ejected outside the virial radius, with less metals being ejected in the lower resolution runs.

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description