Modeling for Stellar Feedback in Galaxy Formation Simulations
Alejandro Núñez, Jeremiah P. Ostriker, Thorsten Naab, Ludwig Oser, Chia-Yu Hu, Ena Choi
DDraft version January 5, 2017
Preprint typeset using L A TEX style emulateapj v. 5/2/11
MODELING FOR STELLAR FEEDBACK IN GALAXY FORMATION SIMULATIONS
Alejandro Núñez , Jeremiah P. Ostriker , Thorsten Naab , Ludwig Oser , Chia-Yu Hu and Ena Choi Draft version January 5, 2017
ABSTRACTVarious heuristic approaches to model unresolved supernova (SN) feedback in galaxy formation sim-ulations exist to reproduce the formation of spiral galaxies and the overall inefficient conversion of gasinto stars. Some models, however, require resolution dependent scalings. We present a sub-resolutionmodel representing the three major phases of supernova blast wave evolution —free expansion, energyconserving Sedov-Taylor, and momentum conserving snowplow— with energy scalings adopted fromhigh-resolution interstellar-medium simulations in both uniform and multiphase media. We allow forthe effects of significantly enhanced SN remnant propagation in a multiphase medium with the cool-ing radius scaling with the hot volume fraction, f hot , as (1 − f hot ) − / . We also include winds fromyoung massive stars and AGB stars, Strömgren sphere gas heating by massive stars, and a gas coolinglimiting mechanism driven by radiative recombination of dense H ii regions. We present initial testsfor isolated Milky-Way like systems simulated with the Gadget based code SPHgal with improvedSPH prescription. Compared to pure thermal SN input, the model significantly suppresses star for-mation at early epochs, with star formation extended both in time and space in better accord withobservations. Compared to models with pure thermal SN feedback, the age at which half the stellarmass is assembled increases by a factor of 2.4, and the mass loading parameter and gas outflow ratefrom the galactic disk increase by a factor of 2. Simulation results are converged for a two order ofmagnitude variation in particle mass in the range (1.3–130) × solar masses. Subject headings: methods: numerical – galaxies: formation – galaxies: evolution – galaxies: stellarcontent – galaxies: abundances INTRODUCTION
Current cosmological simulations of baryonic and darkmatter (DM) have made great advances in recreat-ing the epochs and spatial distribution of galaxies (seeBertschinger 1998; Springel 2010; Somerville & Davé2015, for reviews). A major challenge is the reproduc-tion of the overall inefficiency of galaxy formation withphysically self-consistent models. Typically, there is atendency for gas in simulated galaxies to transform intostars too efficiently and too early, resulting in stellar massfractions larger than observed and stellar populationsolder than observed (e.g., White & Frenk 1991; Kerešet al. 2009; Guo et al. 2010; Moster et al. 2013).In the case of low mass galaxies, early work has pointedat efficient feedback from supernova (SN) explosions asthe appropriate solution, as it would suppress stellar for-mation at early epochs by heating up and diffusing gas,resulting in systems with lower stellar mass fractionsand younger stellar ages (White & Frenk 1991). Themost widely used feedback is SN-driven outflows (Dekel& Silk 1986; Larson 1974), and such flows are detected inmany observations at low and high redshifts (e.g., Mar-tin 1999; Strickland et al. 2000; Rupke et al. 2005; Steidelet al. 2010; Newman et al. 2012). Considering the largevariations in numerical implementations to create galac-tic outflows (e.g., Springel & Hernquist 2003; Dubois & Department of Astronomy, Columbia University, 550 West120th Street, New York, NY 10027, USA Princeton University Observatory, Princeton, NJ 08544,USA Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Strasse 1, D-85741 Garching, Germany Department of Physics and Astronomy, Rutgers, The StateUniversity of New Jersey, Piscataway, NY 08854, USA
Teyssier 2008a; Hopkins et al. 2014; Davé et al. 2013;Teyssier et al. 2013; Vogelsberger et al. 2014; Kimm et al.2015; Schaye et al. 2015; Simpson et al. 2015) in galacticand cosmological simulations, it appears that it is notyet clear how best to implement SN feedback into sim-ulations with limited resolution, which cannot fully re-solve individual SN explosions and turbulent interstellarmedium (ISM).The simplest implementation of SN feedback is to in-ject thermal energy locally to the surrounding medium.This energy, however, is quickly radiated away in denseregions and the blast wave evolution is followed incor-rectly (Navarro & White 1993) (see Stinson et al. 2006,for attempts to follow the blast wave evolution by turningoff “feedback” or “radiative losses”). This situation is anartifact of low resolution simulations, preventing previ-ous SN explosions from vacating their surroundings, thuskeeping subsequent SN explosions in dense areas (see e.g.,Dalla Vecchia & Schaye 2012; Kimm et al. 2015). Somestudies attempt to avoid this by applying a stochasticstellar feedback that is, first, calibrated to match theobserved galaxy stellar mass function and, second, dis-tributed over less mass than heating all neighbors, whichresults in longer cooling times (see studies with the Evo-lution and Assembly of Galaxies and their Environments—
Eagle — simulation, e.g., Schaye et al. 2015; Furlonget al. 2015).A second order approach consists of transferring partor all of SN energy as kinetic energy to the surroundingmedium by conserving momentum from the initial explo-sion (see e.g., Dubois & Teyssier 2008b). In this scenario,however, most of the energy transferred might still endup in thermal form, since the ratio of transferred kinetic a r X i v : . [ a s t r o - ph . GA ] J a n A. Núñez et al.energy to total SN energy is proportional to the ratio oftransferred SN ejecta mass to the mass of the surround-ing medium when momentum conservation is imposed(see e.g., Hu et al. 2014; Kimm et al. 2015). Thus, whilethis approach explicitly includes the momentum of theejecta, it effectively evolves to the usual “thermal feed-back” algorithm. Therefore, overcooling still plays a partin hindering SN energy feedback from noticeably affect-ing the structure and evolution of simulated galaxies.To prevent overcooling from nullifying SN energy feed-back, recent studies adopted approximate approaches,such as decoupling SN-affected gas from hydrodynamicinteractions (e.g., Springel et al. 2005; Oppenheimer &Davé 2008; Vogelsberger et al. 2013, 2014), implement-ing decoupled multi-phase ISM (Scannapieco et al. 2006;Aumer et al. 2013), suppressing radiative cooling effectson SN-heated gas for several tens of Myr (Stinson et al.2006; Guedes et al. 2011), constructing feedback mech-anisms based on the evolution of superbubbles (Kelleret al. 2014, 2015), and stochastically enhancing heatingof SN-affected gas (Dalla Vecchia & Schaye 2012; Schayeet al. 2015). As encouraging as the results are at repli-cating SN-driven winds and star formation rates (SFR),the methods have limited predictive power as they arenot based on ab initio physical modeling of the stellarfeedback process.In this study we attempt to develop a numerical stel-lar feedback prescription based as closely as possible onstandard physical principles. It includes winds and heat-ing from young massive stars as well as feedback fromSN explosions, the latter following the known physicalprocesses that govern the different phases of SN rem-nants. We use the smoothed particle hydrodynamics(SPH) code SPHGal (Hu et al. 2014), an implementationthat incorporates several improvements into the
Gad-get code (Springel 2005) to address some known SPHproblems. Our approach is most similar to the stellarfeedback in Agertz et al. (2013) and Agertz & Kravtsov(2016), as well as the Feedback in Realistic Environ-ment (
Fire ) project, presented in Hopkins et al. (2014)and used by Hopkins et al. (2012, 2013), Muratov et al.(2015), and Su et al. (2016). We apply our stellar feed-back to an isolated Milky Way-like galaxy within an iso-lated dark halo. The tools developed here are in use bythe cosmological code in Choi et al. (2016).We break our approach into three parts: feedback fromyoung, pre-SN stars, feedback from dying low-mass stars,and feedback from SN events. (1) We include the effectsof stellar winds from young stars by transferring momen-tum to neighboring gas particles for the first few Myr.We also heat gas particles within the Strömgren sphereof young star particles, mimicking the effects of ioniz-ing radiation from massive stars (see e.g., Renaud et al.2013). (2) We transfer momentum from old, low-massstars to the surrounding gas throughout several Gyr af-ter a type II SN has occurred, to imitate the slow windsfrom asymptotic giant branch (AGB) stars.(3) In our SN feedback model, we transfer SN energy tothe closest gas particles to the SN event, each one receiv-ing mass and energy according to a computed high res-olution evolution of SN remnants (SNR). Detailed, veryhigh resolution hydrodynamic simulations of SNR evolu-tion (see Li et al. 2015, for details and a review) showthat they pass through three successive phases: free ex- pansion (FE), energy-conserving Sedov-Taylor (ST), andmomentum-conserving snowplow (SP) (see also Kim &Ostriker 2015; Martizzi et al. 2015; Walch & Naab 2015).In our model, which we call the SP-model, each neigh-boring gas particle receives energy following the physicsgoverning the SNR phase in which the gas particles lie.If the gas particle lies within the radius of the ejecta-dominated FE phase, SN energy is transferred by con-serving the ejecta momentum. If the gas particle liesoutside the FE radius but within the radius of the adia-batic blast wave ST phase, SN energy is transferred 30%as kinetic energy and 70% as thermal energy. Finally,if the gas particle is beyond the ST radius and into therealm of the pressure driven and momentum conservingSP phase (Chevalier 1974), then part of the original SNenergy is dissipated and only a fraction is transferredto the gas particle, with the thermal portion dissipatingmore quickly than the kinetic portion following a pre-scription based on the detailed, high resolution, hydro-dynamic treatments of SNR evolution in Li et al. (2015).Our treatment of feedback in the SP phase allows for amulti-phase medium by allowing for the volume fractionthat is hot and tenuous.The paper is organized as follows: in Sec. 2 we de-scribe our numerical code simulation setup and feedbackimplementation, in Sec. 3 we present our results, and inSec. 4 we discuss these results and possible further de-velopments. SIMULATIONS
Numerical Code
To test our SN feedback treatment, we used a recentlymodified version of
Gadget (Springel 2005). The im-proved version —SPHGal— is described in detail byHu et al. (2014). SPHGal includes an artificial viscos-ity implementation modified from the version presentedin Cullen & Dehnen (2010) that uses a shock indicatorbased on the time derivative of the velocity divergenceso as to improve the shock treatment. It also includes anenergy diffusion implementation based on Read & Hay-field (2012) to correct for mixing instabilities. Finally, itincludes a switch from a density–entropy to a pressure–entropy formulation to more accurately compute contactdiscontinuities as described in Hopkins (2013), Ritchie& Thomas (2001), and Saitoh & Makino (2013). As inthe original
Gadget code, SPHGal calculates the gasdynamics using the Lagrangian SPH technique (see e.g.,Monaghan 1992) and it ensures the conservation of en-ergy and entropy (Springel & Hernquist 2002).To complement the original
Gadget , in which star for-mation and cooling was computed for a primordial com-position of hydrogen and helium (Katz et al. 1996), weincluded metal-line cooling using the method and coolingrates described in Aumer et al. (2013), where net cool-ing rates are tracked for 11 species: H, He, C, Ca, O, N,Ne, Mg, S, Si, and Fe. This method exposes the gas tothe redshift-dependent UV/X-ray and cosmic microwavebackground modeled in Haardt & Madau (2001) and al-lows for mixing and spreading of metal abundances in theenriched gas. The code assumes the gas to begin withzero metallicity and to be optically thin and in ionizationequilibrium.
Star Formation Model tellar Feedback in Galaxy Simulations 3We included a chemical evolution and stochastic starformation model (Aumer et al. 2013) that considersenrichment by type II SNe, type Ia SNe, and AGB-type winds following the chemical yield prescriptions ofWoosley & Weaver (1995), Karakas (2010), and Iwamotoet al. (1999), respectively (see Secs. 2.5.2–2.6). As de-scribed in Aumer et al. (2013), SFR is set by dρ ∗ /dt = ερ gas /t dyn , where ρ ∗ and ρ gas are the stellar and gas den-sities, respectively, ε is the star formation efficiency, and t dyn = 1 / p πGρ gas is the local dynamical time for thegas particle. This is applied to gas particles in a conver-gent flow that are denser than a density threshold ρ th ,which we define as ρ th = ρ (cid:18) T gas T (cid:19) (cid:18) M M gas (cid:19) , (1)where ρ is a critical threshold density value, T is a criti-cal temperature threshold, M is a low-resolution particlemass fiducial value, and T gas and M gas are the temper-ature and mass of the gas particles, respectively. Thevalues of ρ , T , and M used in our simulations arelisted in Table 2. This scaling is modeled on the require-ment that the density must exceed the value for the Jeansgravitational instability of a mass M gas at temperature T gas . Star particles are then created stochastically withone gas particle being turned into one star particle of thesame mass. Radiative Recombination
To account for radiative recombination processes in thecool ( < ∼ K) ISM, where H and He line cooling dom-inate, we implemented a simple cooling limiting mecha-nism based on the recombination time of dense H ii re-gions. The recombination time for a gas particle is de-fined by t − = n H α B , where n H is the hydrogen numberdensity of the ISM in cm − and α B ≈ . × − cm s − is the effective radiative recombination rate for hy-drogen, assuming a sphere temperature of 10 K (Draine2011).We calculate the gas temperature resulting from thecooling limiting effect at every time step using the fol-lowing definition T n +1 , rec = T n + ∆ t ∆ t + t rec ( T n +1 − T n ) , (2)where T n and T n +1 are the gas particle temperaturesas determined by the cooling rates described in Aumeret al. (2013) (see Sec. 2.1) for the previous and currenttime steps n and n + 1, respectively, and ∆ t is the timestep size. If the resulting recombination-limited temper-ature T n +1 , rec is below 10 K, then the gas particle getsassigned T n +1 , rec instead of T n +1 . The purpose of thislimitation on the cooling rate is to allow for the fact thatthe relevant atomic processes (i.e., H and He line cooling)cannot proceed faster than the recombination time evenif the code assumes (for simplicity) equilibrium states ofionization. Early Stellar Feedback
As described by several stellar evolution models (e.g.,Leitherer et al. 1999), almost as much momentum isemitted into the ISM by winds from young massive stars as by type II SN explosions. Furthermore, radiation fromthese massive stars heats and ionizes the surrounding gas.The effect of this ionizing radiation can be characterizedby the Strömgren sphere (Strömgren 1939) surroundingmassive stars, within which the ISM is ionized and heatedto ≈ K (Hopkins et al. 2012; Renaud et al. 2013).Stinson et al. (2013) presented an “early stellar feed-back” model that adds thermal energy to gas surroundingstars, but only considering the energy in the UV stellarflux —approximately 10% of the total stellar luminosity(Leitherer et al. 1999). To have a more complete modelfor pre-SN stellar feedback, we included two distinct feed-back mechanisms from young massive stars before theyexplode as SN: stellar winds and Strömgren sphere heat-ing (see e.g., Agertz et al. 2013; Renaud et al. 2013).We implemented a simple momentum conservation en-ergy transfer of winds from massive stars, to first approxi-mation using the same amount of ejected mass and ejectavelocity as those of type II SN explosions (see Sec. 2.5.2)evenly spread in time from star birth up to the momentof SN explosion t SN (for a detailed description of stellarwinds see Kudritzki & Puls 2000).We also included Strömgren sphere heating by ad-justing the temperature of cold gas particles within theStrömgren radius of a star particle, defined as R Str = (cid:18) π ˙ Qn α B (cid:19) / , (3)where ˙ Q is the ionizing photon flux and equals 10 ∆ N ,with ∆ N being the number of stars to explode as typeII SN. To first approximation, the temperature of coldgas particles (i.e., with < ∼ K) within the Strömgrensphere of a star particle is adjusted such that T new = T old + (10 − T old ) ∆ tt ion , (4)where T new and T old are the updated and original gasparticle temperatures in Kelvin, respectively. Assumingionization equilibrium, we set t ion = t rec . SN Energy Feedback
SNR Phases in Detail
Our motivation for developing a snowplow phasedriven SN feedback is based on the current spatial andtemporal resolutions of simulations, which are far frombeing able to resolve individual SN events and the en-ergy transfer that occurs in their immediate vicinity intime and space (for a general review of the SNR phases,see Ostriker & McKee 1988). The typical duration ofthe first phase of a SNR —the FE phase— is limited toa time t ST ≈
200 yr, beyond which the Sedov-Taylorphase begins (see equation 39.3 of Draine 2011). Duringtime t < t ST the ejected mass is much greater than theswept up mass and to lowest order it expands at nearlyconstant velocity, conserving momentum as it hits thesurrounding medium.The radius at which the FE phase ends lies at the pointwhere the total mass swept by the SNR blast wave equalsthe mass of the SNR itself. This occurs at a radius r ST = (cid:18) M ej πρ (cid:19) / , (5) A. Núñez et al.where M ej is the mass ejected in the SN event and ρ isthe mass density of the surrounding medium. SPH simu-lations with particles masses ≈ M (cid:12) numerically resultin r ST in the order of tens of parsecs. Clearly, implement-ing SN feedback based on pure momentum conservationof the SN ejecta, as it applies to the FE phase, would behitting the limits of spatial resolution, which is typicallylarger than 100 pc.During the ST phase the blast wave expands adiabat-ically, experiences negligible radiative loses, and is de-scribed by a self-similar solution (Taylor 1950; Sedov1959; Bandiera 1984). The kinetic energy of the expand-ing material is fixed at ≈
30% of the total energy. Ra-diative cooling becomes significant when the blast waveexpands beyond the radius r cool ’ (19 . E / n − / (1 − f hot ) − / , (6)where E is the kinetic energy of the ejecta normalizedto 10 erg, n is the number density of the surround-ing (primarily) hydrogen medium in cm − , as describedin Blondin et al. (1998), and f hot is the volume fractionof the surrounding medium that is hot and has a verylow density. The dependence of r cool on (1 − f hot ) − / is obtained from analysis of the high resolution simula-tions of Li et al. (2015). That work also found that inself-consistently generated multiphase media there was astrong correlation between the volume averaged temper-ature h T i and the filling factor of the hot gas, which iswell approximated by f hot = 1 − exp [ − ( h T i /T f ) / ] , (7)where T f = 2 × K. With this formulation, f hot tendsto zero in cold regions and to near unity in hot regions.Correspondingly, r cool ∝ exp ( h T i / ), and so this tran-sition radius increases with increasing temperature. SPHsimulations with particles masses ≈ M (cid:12) numericallyresult in r cool values at least ∼ r ST but can increase byorders of magnitude for low n or large f hot values.Beyond r cool the SNR enters the SP phase, in which theblast wave experiences a significant energy loss via radia-tive cooling and expands with relatively constant radialmomentum (Cox 1972; Cioffi et al. 1988; Petruk 2006).Residual thermal energy within the blast wave allowsthe momentum to increase slightly in this phase and, ac-cording to Blondin et al. (1998), is proportional to t / .Based primarily on the detailed hydrodynamic simula-tions of Li et al. (2015), we describe the thermal energyto transfer from the SNR to the surrounding medium as E Th = 0 . E SN . r/r cool ) , (8)where E SN is the total original kinetic energy of the SNRand r is the distance in pc to the medium. The 0.7 factoraccounts for the fact that, at the start of the SP phase,the thermal fraction of the total energy is ≈ Blondin et al. (1998) tested SN remnant propagation in the n range 0.084–8400 cm − . They found unrealistic results only forthe highest-density case, which was driven by their assumption ofno magnetic fields. This approximation is valid mostly for gas with T > ∼ × K. Similarly, we describe the kinetic energy to transferfrom the SNR to the surrounding medium as E K = 0 . E SN . r/r cool ) . (9)Note that since the exponent of the radius fraction issmaller in Eq. 9 than in Eq. 8, the late stages of the SPphase are dominated by its kinetic energy.The resulting deceleration parameter is vt/r ≈ ≈ v is the velocity of the blast wave. Thevalues that we have adopted are consistent with thoserecently found by Kim & Ostriker (2015) (see also Walch& Naab 2015; Martizzi et al. 2015), who also found amore rapid deceleration and energy loss than the valueof 0.33 found by Blondin et al. (1998) (see also Walch &Naab 2015; Girichidis et al. 2016, for ISM simulations ac-counting for the effect of feedback during the SP phase). Construction of a New SN Feedback Model
We construct a feedback formulation that considers thethree different SNR phases described in Sec. 2.5.1, andwe call it
SP-model feedback, since it includes the SPphase of SNR evolution. This model can be applied toboth type II and Ia SN explosions. We allow the modelto transfer SN ejecta energy to the closest 10 neighboringgas particles according to the physics of the SNR phase inwhich the neighboring particles happens to lie accordingto the following recipe: • If the distance between the gas particle and theSN event is smaller than r ST (Eq. 5), SN energyis transferred by conserving momentum, which re-sults in a small fraction of the original SN energybeing transferred as kinetic energy and the rest asthermal energy; we call this FE feedback . • If the distance lies between r ST and r cool (Eq. 6),30% of the SN energy is transferred as kinetic en-ergy and the rest as thermal energy; we call this ST feedback . • If the distance is larger than r cool , SN energy istransferred as both kinetic and thermal energy fol-lowing Eqs. 8 and 9; we call this SP feedback .The total energy liberated in a typical Type II SN ison the order of 10 erg, but much of the loss is emittedas neutrinos and only 1% goes out as kinetic energy withthe ejecta, or ∼ erg. This canonical value for SNenergy in the ejected material is described as E SN = 12 M ej v , (10)where v ej is the velocity of the ejecta (Janka et al. 2007).In our simulations we calculate M ej following the massyield prescriptions for type II SNe by Woosley & Weaver(1995) and for type Ia SNe by Iwamoto et al. (1999), theformer resulting in a typical ejected mass of ≈
13% of theexploding star. Since only stars with masses between 8and 100 M (cid:12) experience type II SNe, the expected valuefor v ej falls in the range 3000–10000 km s − . We choose v ej = 4500 km s − as our fiducial value for both type Iaand II SNe.tellar Feedback in Galaxy Simulations 5Note that each SN event in our simulation correspondsto the energy of many simultaneous physical SN events,as mandated by simulation resolution (see, however, Huet al. 2016, for dwarf galaxy simulations resolving indi-vidual SNe). For a star particle of mass ∼ M (cid:12) , thetypical value of E SN is in the order of 10 erg. Fortu-nately, the same way that Eq. 10 scales with particlemass, so do Eqs. 5–9, making the SP-model formulationformally independent of simulation resolution. Winds From AGB Stars
As part of the SP-model feedback, we also include feed-back from the low-mass end of the stellar population inthe form of slow winds from AGB stars. Similar to ourtreatment of winds and metals from young massive stars(see Sec. 2.4), we transfer energy from the star particlesto the surrounding gas particles by conserving momen-tum of the ejected material, which is calculated using themass yield prescriptions for AGB by Karakas (2010). Weassume an AGB wind velocity of v AGB = 10 km s − . Isolated Galaxy Simulation
With our simulated galaxy we are recreating a MilkyWay-like spiral galaxy isolated from any outside inter-ference to isolate the effects of stellar feedback on theevolution of the galaxy dynamics and structure. Thegalaxy consists of a DM halo, a rotationally supporteddisk of gas and stars, a central stellar bulge, and a hotgas halo. We followed the method described in Springelet al. (2005) and extended by Moster et al. (2011) for ahot gas halo to set up close to equilibrium initial condi-tions.There are 3 × DM particles with mass 4 . × M (cid:12) , for a total DM mass of 1 . × M (cid:12) and theyare arranged so as to follow the mass distribution pro-file described by Hernquist (1990) with a concentrationparameter c = 9 (Navarro et al. 1997).The total baryonic mass is 5 . × M (cid:12) and is rep-resented by three different galactic structures: a disk ofgas and stars, a stellar bulge, and a hot ( > ∼ K) gashalo. Almost 80% of the disk is in a stellar disk, whilethe other 20% is represented by 9 . × collisional SPHgas particles each of mass 1 . × M (cid:12) . The total massof the stellar bulge is 1 . × M (cid:12) . The total mass ofthe hot gas halo is 2 × M (cid:12) , approximately 4% of thetotal baryonic mass of the galaxy.We include this diffuse, rotating, hot halo —as de-scribed in Moster et al. (2011) (see also Choi et al.2014)— to recreate observed, X-ray emitting, extendedhot haloes around moderate mass, spiral galaxies (An-derson et al. 2013). The gas halo assumes that thereis no angular momentum transport between it and thespherical DM halo. The median temperature of the hothalo gas is 1 . × K and the mean hydrogen densityis 6 . × − cm − , similar to the values found for theMilky Way Galaxy by Faerman et al. (2016).The total galactic mass is 1.4 x 10 M (cid:12) . Up to 95%of the total galactic mass is contained within 1 Mpc fromthe galactic center, and 100%, within 8 Mpc. The sidelength of the cube box is 50 Mpc. Table 1 summarizesSPH particle characteristics for our fiducial initial condi-tions.The galaxy is evolved for 6 Gyr using the fixed grav- Table 1
SPH Particles Description of Initial Condition for our FiducialMilky-Way Like Isolated Galaxy.Mass MassParticle Number per particle Total (cid:15) a type (10 ) (10 M (cid:12) ) (10 M (cid:12) ) (kpc)Halo gas 1 . . . . . . . . . . . . . . . . . . . . . . a Gravitational softening length. itational softening lengths specified in Table 1 for eachparticle type. The gas is converted to stars accordingto the prescription in Sec. 2.2 and using the fiducial pa-rameter values described in Table 2. The same softeninglength that applies to gas particles also applies to newstar particles (see Table 1). Star particles then stochasti-cally experience one type II SN event after an age t SN = 3Myr. They also experience continuous type Ia SN eventsand AGB feedback starting at age 50 Myr and every 50Myr after that (see Secs. 2.5.2–2.6). The disk and bulgestars from the initial conditions do not contribute to anyfeedback at any point in the simulation. All input pa-rameters for the simulation are summarized in Table 2.Lastly, mass and energy emanating from a star par-ticle from one of the events described in Sec 2.4–2.6 isdistributed among neighboring gas particles according totheir distance from the star particle using weights calcu-lated by the Wendland C kernel (Dehnen & Aly 2012).All velocities added to neighbor gas particles point ra-dially away from the location of the star particle. Thisis true for stellar winds from O and B stars, Strömgrensphere heating, type II and Ia SNe, and AGB winds. Wefind that at the end of our simulations the mass rangeof the newly formed star particles is 0.7–13.4 × M (cid:12) ,with a median value of 0.9 × M (cid:12) . The details of theimplementation of the SP-model to SPH simulations areincluded in Appendix A. RESULTS
In addition to a simulation using the SP-model feed-back formulation, we ran two other simulations with thesame initial conditions but different feedback schemes:
NO feedback , in which we apply mass transfer butno energy transfer at all, and
Hot bubble feedback, inwhich we transfer SN energy in purely thermal form tothe five closest neighboring particles. Note that, by con-struct, the amount of total kinetic energy transferredis null in these two alternative scenarios. We comparethe results from these two alternative scenarios to theSP-model results to better evaluate the effect of theSP-model. The most relevant results from our isolatedgalaxy simulations are summarized in Table 3.
SFR and Galactic Characteristics
The most important effect of the SP-model feedbackcan be observed in the star formation history. Fig. 1shows the SFR for our isolated galaxy using all threeSN feedback schemes considered. The SFR histories for A. Núñez et al.
Table 2
Parameters Used in Simulation.Symbol Fiducial Value Note
Star Formation ε ρ − Critical density threshold T M × M (cid:12) Low resolution particle mass fiducial value
Radiative Recombination and Early Stellar Feedback α B × − cm s − Effective radiative recombination rate for hydrogen t SN × yr Lifetime of massive stars before they explode as type II SN SN and AGB Feedback v ej − Velocity of ejected SN mass and stellar wind mass v AGB
10 km s − Velocity of ejected AGB mass T f × K Temperature factor for the calculation of f hot . S F R ( M fl y r − ) NO feedbackHot bubbleSP-model
Figure 1.
SFR as a function of time for our Milky Way-like iso-lated galaxy using three different SN feedback schemes: no energyfeedback (NO feedback), purely thermal feedback to the five closestneighboring particles (Hot bubble), and the model SN energy feed-back (SP-model), which transfers energy according to the physicsof the SN phase (free expansion, Sedov-Taylor, or snowplow) inwhich each of the closest 10 neighboring particles lies. The verticalarrows indicate the age at which half the new galactic stellar masshas been assembled. the three feedback schemes differ more significantly dur-ing the early phases of the simulations, when the pres-ence of feedback in kinetic form suppress SFR more effi-ciently. The SP-model feedback, which most effectivelysuppresses early star formation, has the highest rate oflate-time star formation, as material partially ejected bya central “fountain” returns at later times to be incor-porated in late forming stars. Under NO feedback, thecurrent SFR is 0.23 M (cid:12) yr − , compared to a value of0.94 under Hot bubble feedback and 1.15 M (cid:12) yr − un-der SP-model feedback. The estimated value for theMilky Way is 0.68–1.45 M (cid:12) yr − (Robitaille & Whitney2010). Had we allowed for growing torques from non-axisymmetric tidal forces, these effects would have beenmore pronounced (see e.g., Brook et al. 2014; Übler et al.2014).The differences in star formation affect several othergalactic characteristics. We therefore measure other im-portant galactic quantities to better understand the im-pact of each feedback scheme. These measurements areshown in Table 3.The strongly suppressed early star formation resultsin an overall younger stellar population, as indicated bythe age at which half the new galactic stellar mass is assembled: 1.48 Gyr for SP-model feedback as opposedto 0.37 Gyr for NO feedback and 0.63 for Hot bubblefeedback. These half ages are indicated by vertical arrowsin Fig. 1.The gas outflow rate (OFR) —measured as the cumu-lative outflow rate of gas beyond ± (cid:12) yr − underNO feedback, 0.37 M (cid:12) yr − under Hot bubble feedback,and 0.76 M (cid:12) yr − under SP-model feedback.The mass loading factor η —measured as the OFR di-vided by the cumulative stellar mass produced over 6Gyr— is a very important metric for assessing the ap-plicability of the models. For the NO feedback case itis 0.01, for the Hot bubble feedback case it is 0.36, andfor the SP-model feedback case it is 0.72. Muratov et al.(2015) found in their study of galactic winds using high-resolution cosmological simulations —part of the Fire project— that galaxies similar in size to our fiducial ex-hibit gusty outflows at high redshift ( z > .
0) but thenreach a moderate, steady SFR and subdued, weak out-flows at low redshift ( z < . η values close to 10 at high redshift and lower than 1 atlow redshift .The different SFR histories of the three feedback sce-narios result in different stellar populations after 6 Gyr.Table 3 includes final halo and galactic characteristics.We define the virial radius to be r , the radius wherethe mean density drops below 200 times the critical den-sity of the universe. We then define the galactic radiusto be r = 0 . r . The variations across the feedbackscenarios of the virial mass and virial radius are small, asexpected. The distinguishable impact of the SP-modellies in the age and mass distribution of the stellar popu-lation.The SP-model feedback generates a flatter disk struc-ture in the inner 4 kpc of the galaxy than NO feedbackand Hot bubble feedback. Fig. 2 shows surface density—using dm/dr = 2 πr Σ( r )— at different galactic radiiat the end of the simulations for all three SN feedbackscenarios. Another indicator of the spread of stellar con-tent is the effective radius r eff —the radius of a spherecircumscribing half of the stellar mass (vertical arrows Muratov et al. (2015) presents several different methods tomeasure η for their simulated galaxies. Our method is most similarto their Cross T method, albeit with different geometry in thethreshold boundary to categorize the outflow. tellar Feedback in Galaxy Simulations 7 Σ ( M fl p c − ) NO feedbackHot bubbleSP-model
Figure 2.
Current surface density curves for our simulated galaxyusing three different SN feedback schemes. The vertical arrowsindicate the effective radius r eff of a sphere circumscribing half ofthe new stellar mass. in Fig. 2). r eff is larger under SP-model feedback NOfeedback by a factor of 3.2, but it is smaller than underHot bubble feedback by a factor of 0.8. We found no barformation in our simulations.Fig. 3 shows the theoretical circular velocities v c = p Gm ( r ) /r of each galactic component for the initialconditions (IC) (top panel), the NO feedback scenario(upper middle panel), the Hot bubble feedback scenario(lower middle panel), and the SP-model feedback sce-nario (bottom panel). The SP-model feedback creates asignificantly flatter rotational curve for the new stellarcomponent (black solid line) compared to NO feedbackand Hot bubble feedback. Under the SP-model feedback,the mean v c for gas and star particles within r at theend of the simulation are 24 km s − and 99 km s − ,respectively, compared to 55 km s − and 69 km s − inthe initial conditions. Interestingly, the old bulge stellarcomponent (orange dashed line) ends up with a signif-icantly flatter rotational curve under the NO feedbackand Hot bubble scenarios. We suspect this is caused bya transfer of angular momentum from the bulge star par-ticles to the new star particles, as a significant portion ofthe star formation occurs near the galactic bulge. Thelatter is a direct consequence of the absence of energytransfer in the NO feedback scenario and lack of kineticenergy feedback in the Hot bubble scenario, which wouldotherwise help push out hot gas from star forming re-gions.To compare how stellar populations assemble in thegalaxy, Fig. 4 shows a history of stellar (top panel) andgas (bottom panel) masses within the galactic radius r for all three feedback scenarios. The vertical arrows in-dicate the time when the galaxy has assembled half of itsnew stellar mass. In the NO feedback scenario (dottedline), the new stellar content of the galaxy accumulatesquickly, whereas any form of energy feedback (Hot bub-ble or SP-model) hinders this process significantly. Con-sequently, gas is very quickly depleted under NO feed-back, as seen in the bottom panel. Both the Hot bubbleand SP-model feedback schemes allow more gas to re-main by either heating it or kicking some of it out of thegalaxy early on and allowing it to later rain back intothe disk. After 6 Gyr the gas-to-stellar-mass ratio is sig-nificantly higher when the SP-model is in place (0.16) v c ( k m s − ) r IC new starsgasDMbulge starsdisk starsTotal v c ( k m s − ) r NO feedback v c ( k m s − ) r Hot bubble v c ( k m s − ) r SP-model
Figure 3.
Rotational curves for each galactic component for theinitial conditions (IC) (upper panel), at end of simulation usinga NO feedback scheme (upper middle panel), using a Hot bub-ble feedback scheme (lower middle panel), and using the SP-modelscheme (bottom panel). The galactic radius r = 0 . r is in-dicated with a vertical grey line. The vertical arrows indicate theeffective radius r eff of a sphere circumscribing half of the new stel-lar mass. compared to NO feedback (0.04). Energy Transfer
Figs. 5–7 present a breakdown of the energy transferredin type II SN events. They help to illustrate in moredetail how the SP-model feedback distributes energy tothe surrounding medium.According to the description of the SP-model feedbackin Sec. 2.5.1, some SN energy is radiatively lost dur-ing the SP phase (outside r cool ), whereas no energy islost during the FE and ST phases. Therefore, the total A. Núñez et al. Table 3
Simulations Results.Halo Galaxy Galactic PropertiesSN Feedback r
200 a m
200 b r eff c m
10 d m ∗
10 e m g10 f f ∗ g L X h t / SFR j OFR k η l Type (kpc) (kpc) (%) (erg s − ) (Gyr) (M (cid:12) yr − )NO feedback 162.14 94.89 0.53 18.88 0.75 0.03 4.7 5.1 × × × M (cid:12) . a virial radius; b virial mass; c radius of a sphere circumscribinghalf of the stellar mass; d galactic mass; e galactic stellar mass; f galactic gas mass; g percentage of galactic baryons converted into stars: m ∗ /( m × universal baryon fraction); h current X-ray luminosity; i age at which half the new galactic stellar mass is assembled; j current star formation rate; k cumulative gas outflow rate from galactic disk over 6 Gyr; l mass loadingfactor: cumulative gas outflow from galactic disk divided by the cumulative stellar mass produced over 6 Gyr. N e w M ∗ ( M fl h − ) M g a s ( M fl h − ) NO feedbackHot bubbleSP-model
Figure 4.
New stellar mass (top panel) and gas mass (bottompanel) within the galactic radius r = 0 . r as a function oftime for three different SN feedback schemes. The vertical arrowsindicate the time when the galaxy assembled half of its new stellarmass. amount of energy transferred to gas particles is not al-ways equal to the original SN energy of the ejecta. Fig. 5shows a normalized histogram of the fraction of type IISN energy released by a star particle that is actuallytransferred to neighboring gas particles. In ≈
40% of thecases, the fraction of the original SN energy transferred is80% or greater, either because neighboring gas particlesare in FE or ST phase, where no SN energy is lost, or be-cause they are within the SP phase but not too far fromthe ST–SP transition radius r cool , thus allowing only asmall fraction of the original SN ejecta energy to be lostin cooling. By construct, most of the energy in the FEphase, ST phase, and SP phase near r cool is in thermalform. On the other hand, ≈
10% of the cases resultedin an energy transfer less than 20% of the original SN E SN transferred00.030.060.090.12 N o r m a li z e d c o un t Figure 5.
Normalized histogram of the fraction of type II SNfeedback energy transferred to neighboring gas particles. energy transferred to neighboring gas particles. This is aconsequence of gas particles being far into the SP phase,i.e., well outside the ST–SP transition radius r cool , wherea significant portion of the original SN energy is lost viaradiative cooling. Most of the energy in the SP phasefar from r cool is in kinetic form, as the thermal portionof the remaining energy dissipates more quickly than thekinetic portion (see Eqs. 8 and 9). Thus, for the caseswhere most of the original SN energy is dissipated, al-most all of the energy transferred is in kinetic form.Fig. 6 illustrates the fraction of energy transferred atthe different SNR phases as a function of time. Most SNenergy is transferred to gas particles lying within the STphase (red region) for three main reasons: first, no energyloss occurs during the ST phase; second, neighboring gasparticles in the ST phase have larger SPH kernel weightsthan those in the SP phase; and third, r cool increasesas ISM densities drop, the latter driven by SN feedback.Also, not many gas particles are found in the FE phase,which is subject to simulation resolution. Finally, Fig. 6shows that during the first ≈ n values and low f hot values — and thus,small r cool — must have shifted to lower n and higher f hot values rather rapidly, which in turn increased r cool ,such that after ≈ ≈
50% of the total.Fig. 7 shows the fraction of total energy transferredin the form of kinetic and thermal energy as a functiontellar Feedback in Galaxy Simulations 9 F r a c t i o n o f t o t a l E S N t r a n s f e rr e d FE phaseST phaseSP phase
Figure 6.
Fraction of total transferred type II SN energy at eachSNR phase as a function of time when the SP-model energy feed-back scheme is implemented on our Milky Way-like isolated galaxy. F r a c t i o n o f t o t a l E S N t r a n s f e rr e d KineticThermal
Figure 7.
Fraction of total transferred type II SN energy as eitherkinetic (blue) or thermal (orange) form when the SP-model energyfeedback scheme is implemented on our Milky Way-like isolatedgalaxy. The gray dashed line indicates the limit where 30% ofenergy is transferred in kinetic form, which is the case during theST phase of the SP-model feedback scheme. of time. During the first ≈ f hot , the amount of transferred thermaland kinetic energy approximates the type of feedbackpresent in the ST phase, not only because the ST–SPtransition radius r cool increases in size (Eq. 6), but alsobecause the radial factor in Eqs. 8 and 9 approachesunity. We tested a simulated galaxy using exclusivelyST feedback on all neighboring gas particles and foundthat the results were very similar to the results using theSP-model feedback. X-rays
As noted by Boroson et al. (2011), a correlation is ex-pected between the hot gas content and gravitationalpotential of the galaxy. These authors measured X-rayemission by hot and diffuse gas in galactic halos in therange 0.3–8 keV as proxy for the hot gas content, and cen-tral velocity dispersion σ as proxy for the gravitationalpotential, on a sample of 30 early-type galaxies observedwith Chandra .To verify this correlation in our simulated galaxies, weestimate the total galactic X-ray luminosity L X using anemissivity prescription, described in Choi et al. (2014),that includes bremsstrahlung radiation and line emis-sions of all tracked species within the 0.3–8 keV band pro-duced by hot ( T ≥ K) and diffuse ( ρ ≤ . × − g cm − ) gas halo particles. We then measure circularvelocity h v c i as the mass-weighted mean rotational ve-locity for star particles within the galactic disk r . Wecompare our h v c i values against the measured σ in Boro-son et al. (2011). We convert their σ values into h v c i using h v c i ≈ σ . We find that, given the h v c i valuesof our simulated galaxies, L X falls within the observedlocus defined in Boroson et al. (2011). Fig. 8 shows L X versus h v c i of our isolated galaxy for all three feedbackscenarios. Only small L X differences are seen betweenthe three models, with the SP-model feedback resultingin a L X value higher than the Hot bubble feedback bya factor of 2. This is probably caused by the large frac-tion of kinetic energy transferred to the disk gas nearSN events at early epochs (see Fig. 7), which pushes itfarther into the galaxy halo.Also, if we consider the mass– L X relation found byAnderson et al. (2015) (see bottom panel of their figure5) on a large sample of bright galaxies from the SloanDigital Sky Survey (SDSS), we find that given the cur-rent total stellar mass content of our simulated galaxy( ∼ × M (cid:12) ), its halo L X falls within the trend ob-served in the range 10 − M (cid:12) . Metallicity
Since the different SN feedback mechanisms affect dif-ferently the rate and time of star formation, they alsoaffect the final metal content of both gas and star par-ticles in the galaxy. As we described in Sec. 2.2, ejectedmaterial from type Ia and II SNe, as well as from AGBand young massive stellar winds, is composed of severalmetal species in addition to hydrogen and helium.Under NO feedback, the current median metallicityof the new stellar mass log(Z/Z (cid:12) ) of our simulatedgalaxy is 0.11 +0 . − . (upper and lower bounds indicatethe 16th and 84th percentiles); under Hot bubble feed-back, it is 0.26 +0 . − . ; and under the SP-model feedback,0 A. Núñez et al.
100 200 300 › v c fi (km s − )10 L X ( . - k e V ) ( e r g s − ) NO feedback Hot bubbleSP-model
Boroson et al. 2011
Figure 8.
X-ray luminosity L X of our Milky Way-like isolatedgalaxy versus circular velocity h v c i for three different feedbackschemes. We calculate h v c i as the mass-weighted mean rota-tional velocity for star particles within the galactic radius r .To calculate L X we use an emissivity prescription that includesbremsstrahlung radiation and line emissions of all tracked specieswithin the 0.3–8 keV band produced by hot ( T ≥ K) and dif-fuse ( ρ ≤ . × − g cm − ) gas halo particles. The observedrelation (in light gray) is from Boroson et al. (2011). We translatetheir measurements of central velocity dispersion σ into h v c i using h v c i ≈ σ . it is 0.34 +0 . − . . The SP-model feedback increases stellarmetallicity by a factor of 1.7 compared to NO feedback,and by a factor of 1.2 compared to Hot bubble feedback.These values, however, do not account for the fact thatwe do not track metallicity for primordial disk and bulgestars in our galaxy. Therefore, to determine the medianmetallicity for the entire stellar content, we estimate themedian metallicity of this primordial stellar mass usingthe empirical stellar mass–metallicity relation in Gallazziet al. (2005). These authors calculated the median and16 th and 84 th percentiles of the distributions in stellarmetallicity as a function of stellar mass for a sample ofhigh-quality SDSS spectra of galaxies that includes bothearly- and late-type galaxies. Given the total mass of theprimordial bulge and disk stars in our galaxy (5.3 × M (cid:12) ), we find that their metallicity ought to be log(Z/Z (cid:12) )= 0 . +0 . − . according to this relation. Using this value,the mass-weighted stellar metallicity for the entire stellarmass content of our galaxy under SP-model feedback is0 . +0 . − . . This value is statistically similar to the oneobtained using the Gallazzi et al. (2005) relation for agalaxy with the same total stellar content as ours. Weconsider this result only as guidance, since in our isolatedgalaxy scenario we do not account for external galacticevents —such as mergers or infall of primordial cosmolog-ical gas— that could alter the metal content of galaxies.Fig. 9 shows the stellar radial metallicity distributionin the galaxy under NO feedback (blue dotted line), Hotbubble feedback (red dashed line), and SP-model feed-back (solid black line). A linear regression analysis inthe inner part of the disk ( r/r eff <
5) indicates thatthe stellar metallicity gradient becomes shallower in theSP-model feedback case ( − . ± .
01) compared to theNO feedback case ( − . ± .
01) and Hot bubble case( − . ± . − . ± . r/r eff l o g ( Z / Z fl ) S ´ a n c h ez - B l ´ a z q u ez + NO feedbackHot bubbleSP-model
Figure 9.
Radial distribution of the stellar metallicity contentin the galaxy under three different SN feedback schemes. Resultsfrom Sánchez-Blázquez et al. (2014) on samples of galaxies fromthe CALIFA survey are indicated with a gray line. We used thezero point from our own linear regression analysis to indicate theirresult. the zero point from our fit to over-plot their result (grayline) in Fig. 9.The median galactic gas metallicity increases from12 + log(O/H) = 8.91 +0 . − . under NO feedback, to9.53 ± .
01 under Hot bubble feedback, to 9.61 ± .
02 un-der the SP-model feedback. Tremonti et al. (2004) founda median gas metallicity for galaxies similar in stellarcontent to ours of 12 + log(O/H) = 8 . +0 . − . in theirsample of SDSS galaxies at 0 . < z < .
25. Similarly,Andrews & Martini (2013) measure 12 + log(O/H) =8.75 ± .
02 for galaxies similar in stellar content to ours.The higher-than expected gas metallicity in our sim-ulated galaxy with SP-model feedback may also reflectthe inclusion of galaxies at higher redshifts and of early-type galaxies in these empirical samples. Nevertheless,the trend of higher gas metallicity with the SP-modelsuggests that more metal content is reaching more gasin the disk of the galaxy. The gradient of gas metallic-ity decreases from − . ± .
01 (NO feedback) to ≈ -0.01(SP-model feedback). Sánchez et al. (2014) found fortheir sample of 306 galaxies from the CALIFA survey acharacteristic slope of − . ± .
09 (with a zero-point of12+log(O/H) = 8 . ± .
16 between 0.3 and 2 r eff ).There is an important effect that may also explain ourhigher than expected gas metal content. In our ownMilky Way galaxy the deuterium abundance is close to80% of the primordial abundance, indicating that inflow-ing pristine and metal-free gas makes up a non negligiblefraction of the observed local ISM. Thus, if we made al-lowance for cosmological infall of gas in our simulations,the metal abundance in our galaxies would be signifi-cantly lower. Sensitivity to Ejecta Velocity Parameter
To explore the sensitivity of the SP-model scheme tothe ejecta velocity parameter v ej , we re-ran our simula-tions with the SP-model feedback changing v ej from ourfiducial value of 4500 km s − to 3000, 6000, and 7500km s − . We find that the baryon conversion efficiency f ∗ — defined as m ∗ / ( m × f bar ), where f bar ≈ . v ej increases: f ∗ goes from 0.042 to 0.032and 0.024 for v ej = 3000, 6000, and 7500 km s − , respec-tively. Under the NO feedback scenario —which servestellar Feedback in Galaxy Simulations 11 Table 4
Comparison of SPH Particles Parameters at Different ResolutionLevels.Low Res High ResPar. N a Mass b (cid:15) c N a Mass b (cid:15) c Type (10 ) (10 ) (kpc) (10 ) (10 ) (kpc)Gas Halo . .
59 0.13 0.02
Disk . .
46 0.13 0.02Stars
Disk . . Bulge . . . . a Number of particles; b per particle, in units of M (cid:12) ; c gravitational softening length. as a proxy for v ej = 0 km s − , f ∗ = 0 . η increases with increasing v ej , from0.39 to 0.95 and 1.72 for v ej = 3000, 6000, and 7500km s − , respectively. Under the NO feedback scenario, η = 0 .
01. This highlights the sensitivity of the fountaineffect to the parameter v ej , particularly when v ej has avalue greater than our fiducial. Lastly, the galactic halois significantly affected as we increase v ej beyond our fidu-cial value: the galactic halo L X increases by a factor of11 when we increase v ej from 4500 to 6000 km s − , andby a factor of 19 when we increase it to 7500 km s − .There appears to be an energy feedback threshold —dependent on v ej according to Eq. 10— beyond which anascent galaxy at high redshift is strongly compromisedby SNR winds. Clearly, this threshold is the point atwhich the SN ejecta velocity starts reaching the escapevelocity of the system. The ejecta velocity would thus setthe limit of the halo mass below which the resulting stel-lar mass would be very small. The peak ratio of stellarto DM mass (Guo et al. 2010) occurs at stellar mass of2.2 × M (cid:12) , with a strongly observed decline for lowermass systems, which have lower escape velocities. Comparison of Results at Different ResolutionLevels
To test convergence of results at different numericalresolution levels, we ran low and high resolution simu-lations in addition to our fiducial resolution simulationusing the SP-model feedback with 10 times lower and 10times higher mass resolution, respectively. Table 4 de-scribes the SPH particle parameters for the low and highresolution levels.Fig. 10 compares the SFR amongst the three resolu-tion levels tested. We find that to first order, SFR con-verges with resolution, with the high-resolution case hav-ing slightly more spikes of star formation at very earlyepochs than the other two cases. Overall, however, thehigh resolution case results in higher stellar and gas con-tent in the galaxy: gas-to-star mass ratio goes from 0.3in the low resolution case to 0.6 in the high resolutioncase. The high resolution case also results in youngerstellar content: t / increases from 1.4 Gyr in the lowresolution case to 1.8 Gyr in the high resolution case.Similarly, to first order, surface density convergeswith resolution. However, higher resolution results in aslightly more spread-out, flatter stellar content. Fig. 11shows the surface density curves for the three resolutioncases becoming slightly flatter with higher resolution. S F R ( M fl y r − ) low resolutionfiducialhigh resolution Figure 10.
SFR as a function of time for our Milky Way-like iso-lated galaxy using the model SN energy feedback SP-model schemeat low resolution (dotted orange line), our fiducial mid resolution(solid black line), and high resolution (dashed purple line) levels.The vertical arrows indicate the age at which half the new galacticstellar mass has been assembled. Σ ( M fl p c − ) low resolutionfiducialhigh resolution Figure 11.
Current surface density curves for our simulatedgalaxy using the model SN energy feedback SP-model scheme atlow, mid (fiducial), and high resolution levels. The vertical arrowsindicate the effective radius r eff of a sphere circumscribing half ofthe new stellar mass. Also, r eff (vertical lines in Fig. 11) increases from 1.3kpc in the low resolution case to 2.1 kpc in the high res-olution case. Lastly, the mass loading factor η increasesto 5.4 in the low resolution case, and to 1.7 in the highresolution case. Several parameter values are comparedin Table 5 for the three resolution levels. CONCLUSIONS
We implement a physically motivated stellar feedbackmodel, including winds from young and old stars, heat-ing from young massive stars, and a SN energy feedbackformulation (“SP-model”) that allows for cooling of theoriginal SN ejecta, with the thermal portion dissipatingaway more quickly than the kinetic portion. The SP-model mimics the known physics that takes place duringall phase of SNRs —free expansion (FE) phase, Sedov-Taylor (ST) phase, and snowplow (SP) phase— as re-vealed by both observations and high spatial resolutionsimulations. It also allows for the important effects ofsignificantly enhanced SN remnant propagation in a mul-tiphase medium. We implement this on an isolated MilkyWay-type galaxy with a hot gas halo assembled followingthe methods in Springel et al. (2005) using the improved2 A. Núñez et al.
Table 5
Simulations Results at Different Resolution Levels Using the SP-model Feedback.Halo Galaxy Galactic PropertiesResolution r
200 a m
200 b r eff c m
10 d m ∗
10 e m g10 f f ∗ g t / SFR j OFR k η l Level (kpc) (kpc) (%) (Gyr) (M (cid:12) yr − )Low Res 161.53 94.13 1.30 18.93 0.48 0.15 3.0 1.36 1.00 4.29 5.39Fiducial 162.09 94.80 1.71 18.86 0.63 0.10 3.9 1.48 1.15 0.76 0.72High Res 162.80 96.42 2.10 19.45 0.64 0.40 3.9 1.79 1.62 1.78 1.67NOTE. – All masses in units of 10 M (cid:12) . Notes (a) to (l) are defined in Table 3. TreeSPH code SPHGal (Hu et al. 2014), and we evolveit for 6 Gyr.We also implement wind feedback from both youngmassive stars and older AGB stars by conserving mo-mentum of the ejected material. For the former, we ap-proximate the ejected material to be commensurate bothin mass and velocity with that of type II SN explosions,except that we spread it evenly during the first 3 Myr ofthe life of the star particle. For the latter, we use massyield prescriptions found in the literature and an ejectavelocity of 10 km s − , and we implement it throughoutseveral Gyr following type II SN events.We also implement two mechanisms that limit coolingof the ISM gas, one produced by the Strömgren sphereheating surrounding young massive stars, and the otherone produced by the radiative recombination processesin dense H ii regions.In our simulations the SP-model feedback addressesthe sub-resolution physics problem in low resolution SPHsimulations in which particle distances are barely re-solved at the level necessary to allow for FE feedback, thetype that many galaxy simulations include. By formu-lating a type of SN feedback that occurs mainly duringthe later stages of an SNR, namely ST and SP phases,the SP-model feedback allows for a more accurate repre-sentation of the physics involving SN energy feedback.Given that both ST and SP phases allow for ≥ (cid:12) yr − — as we go from Hot bubble feedback to theSP-model feedback.The assembly of the galactic stellar content is morespread out in time in a galaxy with SP-model feedback,as much of the gas kicked out by the SNRs at early timeseventually falls back to the galaxy, the so-called galac-tic “fountain” effect, where it cools off and is subject tostar formation. The SP-model feedback results in a gasoutflow rate from the galaxy of 0.76 M (cid:12) yr − , comparedto 0.01 M (cid:12) yr − under no stellar feedback and 0.37 un- der Hot bubble feedback. This fountain created by theSP-model feedback helps alleviate the issue of over clus-tering of stellar content at the galactic center, resultingin flatter rotational curves for the stellar component ofthe galaxy, flatter surface density curves for the galacticdisk, and an overall larger galactic size, with the effectiveradius increasing by a factor of 3.2 as we go from havingno stellar feedback to the SP-model feedback. Further-more, the SP-model feedback results in a mass loadingparameter η of 0.72, versus 0.36 under Hot bubble feed-back and 0.01 under no feedback.The integrated stellar metallicity of our galaxy in-creases by a factor of 1.7 [measured as log(Z/Z (cid:12) )],and gas metallicity, by a factor of 5.0 [measured as12+log(O/H)] as we go from having no stellar feedback tothe SP-model feedback. The stellar metallicity disk gra-dients turn shallower under SP-model feedback ( − . − . ACKNOWLEDGMENTS
We thank Miao Li for valuable comments and feedbackand Michaela Hirschmann for her insights into observedgalactic metallicities. The simulations used here wererun on Columbia University’s Yeti cluster and the MaxPlanck Computing and Data Facility.
REFERENCESAgertz, O., & Kravtsov, A. V. 2016, ApJ, 824, 79 tellar Feedback in Galaxy Simulations 13
Agertz, O., Kravtsov, A. V., Leitner, S. N., & Gnedin, N. Y.2013, ApJ, 770, 25Anderson, M. E., Bregman, J. N., & Dai, X. 2013, ApJ, 762, 106Anderson, M. E., Gaspari, M., White, S. D. M., Wang, W., &Dai, X. 2015, MNRAS, 449, 3806Andrews, B. H., & Martini, P. 2013, ApJ, 765, 140Aumer, M., White, S. D. M., Naab, T., & Scannapieco, C. 2013,MNRAS, 434, 3142Bandiera, R. 1984, A&A, 139, 368Bertschinger, E. 1998, ARA&A, 36, 599Blondin, J. M., Wright, E. B., Borkowski, K. J., & Reynolds,S. P. 1998, ApJ, 500, 342Boroson, B., Kim, D.-W., & Fabbiano, G. 2011, ApJ, 729, 12Brook, C. B., Stinson, G., Gibson, B. K., Shen, S., Macciò, A. V.,Obreja, A., Wadsley, J., & Quinn, T. 2014, MNRAS, 443, 3809Chevalier, R. A. 1974, ApJ, 188, 501Choi, E., Naab, T., Ostriker, J. P., Johansson, P. H., & Moster,B. P. 2014, MNRAS, 442, 440Choi, E., Ostriker, J. P., Naab, T., Somerville, R. S., Hirschmann,M., Núñez, A., Hu, C.-Y., & Oser, L. 2016, ArXiv e-printsCioffi, D. F., McKee, C. F., & Bertschinger, E. 1988, ApJ, 334,252Cox, D. P. 1972, ApJ, 178, 159Cullen, L., & Dehnen, W. 2010, MNRAS, 408, 669Dalla Vecchia, C., & Schaye, J. 2012, MNRAS, 426, 140Davé, R., Katz, N., Oppenheimer, B. D., Kollmeier, J. A., &Weinberg, D. H. 2013, MNRAS, 434, 2645Dehnen, W., & Aly, H. 2012, MNRAS, 425, 1068Dekel, A., & Silk, J. 1986, ApJ, 303, 39Draine, B. 2011, Physics of the Interstellar and IntergalacticMedium, Princeton Series in Astrophysics (PrincetonUniversity Press)Dubois, Y., & Teyssier, R. 2008a, A&A, 477, 79Dubois, Y., & Teyssier, R. 2008b, in Astronomical Society of thePacific Conference Series, Vol. 390, Pathways Through anEclectic Universe, ed. J. H. Knapen, T. J. Mahoney, &A. Vazdekis, 388Faerman, Y., Sternberg, A., & McKee, C. F. 2016, ArXiv e-printsFurlong, M. et al. 2015, MNRAS, 450, 4486Gallazzi, A., Charlot, S., Brinchmann, J., White, S. D. M., &Tremonti, C. A. 2005, MNRAS, 362, 41Girichidis, P. et al. 2016, MNRAS, 456, 3432Guedes, J., Callegari, S., Madau, P., & Mayer, L. 2011, ApJ, 742,76Guo, Q., White, S., Li, C., & Boylan-Kolchin, M. 2010, MNRAS,404, 1111Haardt, F., & Madau, P. 2001, in Clusters of Galaxies and theHigh Redshift Universe Observed in X-rays, ed. D. M.Neumann & J. T. V. TranHernquist, L. 1990, ApJ, 356, 359Hopkins, P. F. 2013, MNRAS, 428, 2840Hopkins, P. F., Kereš, D., Murray, N., Hernquist, L., Narayanan,D., & Hayward, C. C. 2013, MNRAS, 433, 78Hopkins, P. F., Kereš, D., Oñorbe, J., Faucher-Giguère, C.-A.,Quataert, E., Murray, N., & Bullock, J. S. 2014, MNRAS, 445,581Hopkins, P. F., Quataert, E., & Murray, N. 2012, MNRAS, 421,3522Hu, C.-Y., Naab, T., Walch, S., Glover, S. C. O., & Clark, P. C.2016, MNRAS, 458, 3528Hu, C.-Y., Naab, T., Walch, S., Moster, B. P., & Oser, L. 2014,MNRAS, 443, 1173Iwamoto, K., Brachwitz, F., Nomoto, K., Kishimoto, N., Umeda,H., Hix, W. R., & Thielemann, F.-K. 1999, ApJS, 125, 439Janka, H.-T., Langanke, K., Marek, A., Martínez-Pinedo, G., &Müller, B. 2007, PhysRep, 442, 38Karakas, A. I. 2010, MNRAS, 403, 1413Katz, N., Weinberg, D. H., & Hernquist, L. 1996, ApJS, 105, 19Keller, B. W., Wadsley, J., Benincasa, S. M., & Couchman,H. M. P. 2014, MNRAS, 442, 3013Keller, B. W., Wadsley, J., & Couchman, H. M. P. 2015,MNRAS, 453, 3499Kereš, D., Katz, N., Davé, R., Fardal, M., & Weinberg, D. H.2009, MNRAS, 396, 2332Kim, C.-G., & Ostriker, E. C. 2015, ApJ, 802, 99 Kimm, T., Cen, R., Devriendt, J., Dubois, Y., & Slyz, A. 2015,MNRAS, 451, 2900Kudritzki, R.-P., & Puls, J. 2000, ARAA, 38, 613Larson, R. B. 1974, in IAU Symposium, Vol. 58, The Formationand Dynamics of Galaxies, ed. J. R. Shakeshaft, 191Leitherer, C. et al. 1999, ApJS, 123, 3Li, M., Ostriker, J. P., Cen, R., Bryan, G. L., & Naab, T. 2015,ApJ, 814, 4Martin, C. L. 1999, ApJ, 513, 156Martizzi, D., Faucher-Giguère, C.-A., & Quataert, E. 2015,MNRAS, 450, 504Monaghan, J. J. 1992, ARA&A, 30, 543Moster, B. P., Macciò, A. V., Somerville, R. S., Naab, T., & Cox,T. J. 2011, MNRAS, 415, 3750Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428,3121Muratov, A. L., Kereš, D., Faucher-Giguère, C.-A., Hopkins,P. F., Quataert, E., & Murray, N. 2015, MNRAS, 454, 2691Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490,493Navarro, J. F., & White, S. D. M. 1993, MNRAS, 265, 271Newman, S. F. et al. 2012, ApJ, 761, 43Oppenheimer, B. D., & Davé, R. 2008, MNRAS, 387, 577Ostriker, J. P., & McKee, C. F. 1988, Reviews of Modern Physics,60, 1Petruk, O. 2006, ArXiv Astrophysics e-printsRead, J. I., & Hayfield, T. 2012, MNRAS, 422, 3037Renaud, F. et al. 2013, MNRAS, 436, 1836Ritchie, B. W., & Thomas, P. A. 2001, MNRAS, 323, 743Robitaille, T. P., & Whitney, B. A. 2010, ApJ, 710, L11Rupke, D. S., Veilleux, S., & Sanders, D. B. 2005, ApJS, 160, 115Saitoh, T. R., & Makino, J. 2013, ApJ, 768, 44Sánchez, S. F. et al. 2014, AAP, 563, A49Sánchez-Blázquez, P. et al. 2014, AAP, 570, A6Scannapieco, C., Tissera, P. B., White, S. D. M., & Springel, V.2006, MNRAS, 371, 1125Schaye, J. et al. 2015, MNRAS, 446, 521Sedov, L. 1959, Similarity and Dimensional Methods inMechanics (New York: Academic Press)Simpson, C. M., Bryan, G. L., Hummels, C., & Ostriker, J. P.2015, ApJ, 809, 69Somerville, R. S., & Davé, R. 2015, ARAA, 53, 51Springel, V. 2005, MNRAS, 364, 1105—. 2010, ARA&A, 48, 391Springel, V., Di Matteo, T., & Hernquist, L. 2005, MNRAS, 361,776Springel, V., & Hernquist, L. 2002, MNRAS, 333, 649—. 2003, MNRAS, 339, 289Steidel, C. C., Erb, D. K., Shapley, A. E., Pettini, M., Reddy, N.,Bogosavljević, M., Rudie, G. C., & Rakic, O. 2010, ApJ, 717,289Stinson, G., Seth, A., Katz, N., Wadsley, J., Governato, F., &Quinn, T. 2006, MNRAS, 373, 1074Stinson, G. S., Brook, C., Macciò, A. V., Wadsley, J., Quinn,T. R., & Couchman, H. M. P. 2013, MNRAS, 428, 129Strickland, D. K., Heckman, T. M., Weaver, K. A., & Dahlem, M.2000, AJ, 120, 2965Strömgren, B. 1939, ApJ, 89, 526Su, K.-Y., Hopkins, P. F., Hayward, C. C., Faucher-Giguere,C.-A., Keres, D., Ma, X., & Robles, V. H. 2016, ArXiv e-printsTaylor, G. 1950, Proceedings of the Royal Society of LondonSeries A, 201, 159Teyssier, R., Pontzen, A., Dubois, Y., & Read, J. I. 2013,MNRAS, 429, 3068Tremonti, C. A. et al. 2004, ApJ, 613, 898Übler, H., Naab, T., Oser, L., Aumer, M., Sales, L. V., & White,S. D. M. 2014, MNRAS, 443, 2092Vogelsberger, M., Genel, S., Sijacki, D., Torrey, P., Springel, V.,& Hernquist, L. 2013, MNRAS, 436, 3031Vogelsberger, M. et al. 2014, MNRAS, 444, 1518Walch, S., & Naab, T. 2015, MNRAS, 451, 2757White, S. D. M., & Frenk, C. S. 1991, ApJ, 379, 52Woosley, S. E., & Weaver, T. A. 1995, ApJS, 101, 181
APPENDIX
DESCRIPTION OF THE SP-MODEL FOR SPH SIMULATIONS.
As described in Sec. 2.5.2, the SP-model feedback formulation states that a neighboring gas particle j can receiveenergy from the SN event in one of three different ways, depending on its distance to the star particle relative to thephase transition radii r ST and r cool (Eqs. 5 and 6).The FE-ST transition radius r ST depends on the density of the medium, and the ST-SP transition radius r cool depends on both the density and the temperature of the medium, the latter via the dependency of the hot volumefraction f hot on the gas temperature (Eq. 7). To allow for sharp density and temperature gradients around the starparticle, we calculate both r ST and r cool for each of the closest 10 neighboring gas particles using the gas particledensity ρ j and temperature T j . The total amount of energy assigned to particle j from the SN event is∆ E j = w j E SN , (A1)where w j is the SPH kernel weight assigned to particle j (in our simulations we use the Wendland C kernel). Thetransfer of this energy affects three basic quantities in particle j : its mass, its velocity, and its entropy (via additionof thermal energy). For the mass, one simply adds to the mass of particle j , M j , the value ∆ m j = w j M ej . How theother two quantities change varies from one SNR phase to the other. FE Phase
If particle j is in the FE phase, then momentum conservation applies. The velocity vector to add to particle j ,assuming a perfectly non-elastic collision with the ejecta, is ∆ v FE = (cid:15) (cid:15) ( v ej − v j ) , (A2)where (cid:15) = ∆ m j /M j , v ej is the ejecta velocity and v j is the velocity of particle j before the collision. All three velocitiesare with respect to the star particle’s rest frame. ∆ v FE is added to v j radially away from the star particle.To conserve energy, one must increase the thermal energy of particle j by∆ E Th , FE = 12 µ ( v ej − v j ) , (A3)where µ = ∆ m j M j / (∆ m j + M j ). ST Phase
If particle j is in the ST phase, then the analytical self-similar solution applies, which gives a constant breakdownof SN energy between ∼
30% kinetic and ∼
70% thermal. This means that one must increase the thermal energy ofparticle j by ∆ E Th , ST = 0 . E j . In general, the change in kinetic energy for particle j is ∆ E K = E K , new − E K , old , or∆ E K = 12 ( M j + ∆ m j )( v j + ∆ v ) − M j v j , (A4)where ∆ v is some velocity added to particle j . In the case of ST phase, the velocity vector ∆ v ST to add to particle j is such that ∆ E K , ST = 0 . E j = 12 ∆ m j ∆ v . (A5)Combining Eqs. A4 and A5 —with ∆ v ST in lieu of ∆ v — and working out the vector algebra, one finds that ∆ v ST = α ST v ej ; (A6)here ∆ v ST has the same direction as v ej , and α ST is α ST = q | . (cid:15) (cid:15) v + ( (cid:15) − sin θ ) v j | − v j cos θv ej , (A7)where θ is the angle between v ej and v j . As in the FE phase, All velocities are with respect to the star particle’s restframe, and ∆ v ST is added to v j radially away from the star particle.tellar Feedback in Galaxy Simulations 15 SP Phase
If particle j is in the SP phase, then Eqs. 8 and 9 apply. One must increase the thermal energy of particle j by∆ E Th , SP = 0 . E j . r j /r cool ) , (A8)where r j is the distance between the star particle and particle j . Similarly, one must increase the kinetic energy ofparticle j by ∆ E K , SP = 0 . E j . r j /r cool ) . (A9)Mimicking the steps for the ST phase and defining f SP = 11 + 0 . r j /r cool ) , (A10)one combines Eqs. A1, A4, A9, A10 and 10 to find that the velocity vector ∆ v SP to add to particle j is ∆ v SP = α SP v ej , (A11)where ∆ v SP has the same direction as v ej , and α SP is α SP = q | . f SP (cid:15) (cid:15) v + ( (cid:15) − sin θ ) v j | − v j cos θv ej . (A12)Again, all velocities are with respect to the star particle’s rest frame, and ∆ v SP is added to v jj