1311.2073v1.pdf

Mon. Not. R. Astron. Soc.

000

, 000–000 (0000)

Printed 12 November 2013

(MN L

X style file v2.2)

Galaxies on FIRE (Feedback In Realistic Environments): Stellar

Feedback Explains Cosmologically Inefficient Star Formation

Philip F. Hopkins

∗

, Dušan Kereš

, José Oñorbe

, Claude-André Faucher-Giguère

Eliot Quataert

, Norman Murray

, & James S. Bullock

TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA

Department of Astronomy and Theoretical Astrophysics Center, University of California Berkeley, Berkeley, CA 94720

Department of Physics, Center for Astrophysics and Space Science, University of California at San Diego, 9500 Gilman Drive, La Jolla, CA 92093

Department of Physics and Astronomy, University of California at Irvine, Irvine, CA 92697, USA

Department of Physics and Astronomy and CIERA, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA

Canadian Institute for Theoretical Astrophysics, 60 St. George Street, University of Toronto, ON M5S 3H8, Canada

Canada Research Chair in Astrophysics

Submitted to MNRAS, November, 2013

ABSTRACT

We present a series of high-resolution cosmological zoom-in simulations

of galaxy formation to

=

0, spanning

halo masses

halo

∼

−

and stellar masses

∗

∼

−

. Our simulations include a fully explicit

treatment of both the multi-phase ISM (molecular through hot) and stellar feedback. The stellar feedback inputs

(energy, momentum, mass, and metal fluxes) are taken directly from stellar population models. These sources of

stellar feedback, with

zero

adjusted parameters, reproduce the observed relation between stellar and halo mass up

halo

∼

(including dwarfs, satellites, MW-mass disks, and small groups). By extension, this leads to

reasonable agreement with the stellar mass function for

∗

. We predict weak redshift evolution in the

∗

−

halo

relation, consistent with current constraints up to

We find that the

∗

−

halo

relation in our calculations is relatively insensitive to numerical details, but is sensi-

tive to the feedback physics. Simulations with only supernova feedback fail to reproduce the observed stellar masses,

particularly for dwarf and/or high-redshift galaxies: radiative feedback (photo-heating and radiation pressure) is nec-

essary to disrupt molecular clouds and enable efficient coupling of later supernovae explosions to the gas.

Instantaneous star formation rates agree well with the observed Kennicutt relation, with weak redshift evolution.

The galaxy-averaged Kennicutt relation is very different from the numerically imposed law for converting gas into

stars on small scales in the simulation and is instead determined by self-regulation via stellar feedback. We find

that feedback reduces star formation rates considerably and produces a reservoir of gas that leads to flatter or rising

late-time star formation histories significantly different from the halo accretion history. Feedback also produces large

short-timescale variability in galactic star formation rates, especially in dwarf galaxies. Many of these properties of

galaxy formation with explicit feedback are not captured by common “sub-grid” galactic wind models.

Key words:

galaxies: formation — galaxies: evolution — galaxies: active — stars: formation — cosmology: theory

1 INTRODUCTION

It is well-known that feedback from stars is a critical, yet poorly-

understood, component of galaxy formation. Within galaxies, star

formation is observed to be inefficient in both an instantaneous and

an integral sense.

Instantaneously, the Kennicutt-Schmidt (KS) relation implies

gas consumption timescales of

∼

50 dynamical times (Kennicutt

1998), while the total fraction of GMC mass converted into stars is

only a few percent (Zuckerman & Evans 1974; Williams & McKee

1997; Evans 1999; Evans et al. 2009). Without strong stellar feed-

back, however, gas inside galaxies cools efficiently and collapses

on a dynamical time, predicting order-unity star formation efficien-

cies on all scales (Hopkins et al. 2011; Tasker 2011; Bournaud et al.

2010; Dobbs et al. 2011; Krumholz et al. 2011; Harper-Clark &

Murray 2011).

In an integral sense, without strong stellar feedback, gas in

cosmological models cools rapidly and inevitably turns into stars,

predicting galaxies with far larger masses than are observed (e.g.

∗

E-mail:phopkins@caltech.edu

Katz et al. 1996; Somerville & Primack 1999; Cole et al. 2000;

Springel & Hernquist 2003b; Kereš et al. 2009, and references

therein). Decreasing the instantaneous star formation efficiency

does not eliminate this integral problem: the amount of baryons

in real galactic disks is much lower than that predicted in models

absent strong feedback (essentially, the Universal baryon budget;

see White & Frenk 1991; Kereš et al. 2009). Constraints from in-

tergalactic medium (IGM) enrichment require that many of those

baryons must have entered galaxy halos and disks at some point

to be enriched, before being expelled (Aguirre et al. 2001; Pettini

et al. 2003; Songaila 2005; Martin et al. 2010). Galactic super-

winds with mass-loading

̇

wind

of many times the star formation

rate (SFR) are therefore generally required to reproduce observed

galaxy properties (e.g. Oppenheimer & Davé 2006). Such winds

have been observed ubiquitously in local and high-redshift star-

forming galaxies (Martin 1999, 2006; Heckman et al. 2000; New-

man et al. 2012; Sato et al. 2009; Chen et al. 2010; Steidel et al.

2010; Coil et al. 2011).

However, until recently, numerical simulations have been un-

able to produce winds with large-mass loading factors from an a

priori model (let alone the correct scalings of wind mass-loading

0000 RAS

arXiv:1311.2073v1 [astro-ph.CO] 8 Nov 2013

Hopkins et al.

Figure 1.

Gas in a representative simulation of a Milky Way-mass halo (

m12q

in Table 1). Image shows the projected gas density, log-weighted (

∼

4 dex

stretch). Magenta shows cold molecular/atomic gas (

1000 K). Green shows warm ionized gas (10

K). Red shows hot gas (

K). Each

image shows a box centered on the main galaxy.

Left:

Box 200 kpc (physical) on a side at high redshift. The galaxy has undergone a violent starburst, leading

to strong outflows of hot and warm gas that have blown away much of the surrounding IGM (even outside the galaxy). Note that the “filamentary” structure

of cool gas in the IGM is clearly affected by the outflows.

Right:

Near present-day, with a

∼

50 kpc box. A more relaxed, well-ordered disk has formed, with

molecular gas tracing spiral structure, and a halo enriched by diffuse hot outflows.

with galaxy mass or other properties), nor to simultaneously pre-

dict the instantaneous inefficiency of star formation within galax-

ies. This is particularly true of models which invoke only energetic

feedback via supernovae (SNe), which is efficiently radiated in the

dense gas where star formation actually occurs (see e.g. Guo et al.

2010; Powell et al. 2011; Brook et al. 2011; Nagamine 2010; Bour-

naud et al. 2011, and references therein). More recent simulations,

using higher resolution and invoking stronger feedback prescrip-

tions, have seen strong winds, but have generally found it neces-

sary to include simplified prescriptions for “turning off cooling” in

the SNe-heated gas and/or include some adjustable parameters rep-

resenting “pre-SNe” feedback (see Governato et al. 2010; Macciò

et al. 2012; Teyssier et al. 2013; Stinson et al. 2013; Agertz et al.

2013). This is physically motivated since feedback processes other

than SNe – protostellar jets, HII photoionization, stellar winds, and

radiation pressure – both occur and are critical in suppressing star

formation in dense gas, as well as “pre-processing” gas prior to

SNe explosions so that SNe occur at densities where thermal heat-

ing can have much larger effects (Evans et al. 2009; Hopkins et al.

2011; Tasker 2011; Stinson et al. 2013; Kannan et al. 2013).

Given limited resolution and the complexity of the baryonic

physics, many cosmological models have treated galactic wind gen-

eration and the inefficiency of star formation in a tuneable, “sub-

grid” fashion. This is not to say that the models have not tremen-

dously improved our understanding of galaxy formation! They have

demonstrated that stellar feedback can

plausibly

lead to (globally)

inefficient star formation, constrained the parameter space of al-

lowed feedback models, made predictions for the critical role of

outflows in enriching the IGM, provided possible baryonic solu-

tions to apparent dark matter “problems” (e.g. Pontzen & Gover-

nato 2012), demonstrated the need for “early” feedback from ra-

diative mechanisms beyond SNe alone, and generally created the

framework for our interpretation of observations. However, with

wind models often relying on adjustable parameters, the integrated

efficiency of star formation in galaxies is to some extent tuned “by

hand” and the predictive power is inherently limited. This is partic-

ularly true for studies of gas in the circum-galactic medium (CGM),

a current area of much observational progress – measurements of

the CGM are sensitive to the phase structure of the gas, which is

not faithfully represented in models which simply “turn off” hy-

drodynamics or cooling, or mimic strong feedback via pure thermal

energy injection or “particle kicks.”

Accurate treatment of star formation and galactic winds ul-

timately requires realistic treatment of the stellar feedback pro-

cesses that maintain the multi-phase ISM. Motivated by this philos-

ophy, in Hopkins et al. (2011) (Paper

) and Hopkins et al. (2012d)

(Paper

), we developed a new set of numerical models to fol-

low stellar feedback on scales from sub-GMC star-forming regions

through galaxies. These simulations include the energy, momen-

tum, mass, and metal fluxes from stellar radiation pressure, HII

photo-ionization and photo-electric heating, SNe Types I & II, and

stellar winds (O-star and AGB). Critically, the feedback is directly

tied to the young stars, with the energetics and time-dependence

taken from stellar evolution models. In our previous work, we

showed, in isolated galaxy simulations, that these mechanisms pro-

duce a quasi-steady ISM in which GMCs form and disperse rapidly,

with phase structure, turbulence, and disk and GMC properties in

good agreement with observations (for various comparisons, see

Narayanan & Hopkins 2012; Hopkins et al. 2012b,a, 2013c). In Pa-

per

, Hopkins et al. (2013d), and Hopkins et al. (2013a) we showed

that this leads naturally to “instantaneously” inefficient SF (predict-

ing the KS-law), regulated self-consistently by feedback and

inde-

0000 RAS, MNRAS

000

, 000–000

Galaxies on FIRE: Feedback & Inefficient Star Formation

Figure 2.

Stars in the

m12q

simulation at

∼

0, in a box 50 kpc on a

side near present-time. Image is a mock

composite. The disk is ap-

proximately face-on, and the spiral structure is visible. (The image uses

STARBURST99 to determine the SED of each star particle given its known

age and metallicity, then ray-traces the line-of-sight flux following Hopkins

et al. (2005), attenuating with a MW-like reddening curve with constant

dust-to-metals ratio for the abundances at each point.)

pendent

of the numerical prescription for star formation in very

dense gas. In Hopkins et al. (2012c) (Paper

III

) and Hopkins et al.

(2013b) we showed that the same feedback models reproduce the

galactic winds invoked in previous semi-analytic and cosmological

simulations, and that the

combination

of multiple feedback mecha-

nisms is critical to produce massive, multi-phase galactic winds.

However, our simulations have thus far been limited to ideal-

ized studies of isolated galaxies and galaxy mergers. These previ-

ous calculations thus cannot follow accretion from or interaction

of outflows with the IGM, realistic galaxy merger histories, and

many other important processes. In this paper, the first of a series,

we present the FIRE (Feedback In Realistic Environments) simula-

tions:

a suite of fully cosmological “zoom-in” simulations devel-

oped to study the role of feedback in galaxy formation. To test the

models and understand feedback in a wide range of environments,

we study a wide range in galaxy halo and stellar mass (as opposed

to focusing just on MW-like systems), and follow evolution fully

0. Our suite of calculations includes several of the highest-

resolution galaxy formation simulations that have been run to

Our simulations utilize a significantly improved numerical imple-

mentation of SPH (which has resolved historical discrepancies with

grid codes), as well as the full physical models for feedback and

ISM physics introduced and tested in Paper

-Paper

III

. Here, we

explore the consequences of stellar feedback for the inefficiency

of star formation, perhaps the most basic consequence of stellar

feedback for galaxy formation. In companion papers, we will in-

vestigate the properties of outflows and their interactions with the

IGM, the effect of those outflows on dark matter structure, the dif-

ferences between numerical methods in treating feedback, the role

Movies and summaries of key simulation properties are avail-

able at

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

Movies_cosmo.html

Figure 3.

Gas, as Fig. 1, for a dwarf galaxy (

m10

in Table 1).

Top:

100 kpc

(physical) box, at high redshift.

Bottom:

20 kpc box at intermediate redshift.

Strong outflows are still present, though they are more spherical, because

the galaxy halo is itself small and embedded

within

a much larger filament.

of feedback in determining galaxy structure, and many other open

questions.

In § 2-4, we describe our methodology. § 2 describes the initial

conditions for the simulations; § 3 outlines the implementation of

the key baryonic physics of cooling, star formation, and feedback

(a much more detailed description is given in Appendix A); § 4

briefly describes the improvements in the numerical method com-

pared to past work (again, more details are in Appendix B). And

in Appendix C we test and compare these algorithms with higher-

resolution simulations of isolated (non-cosmological) galaxies.

We describe our results in § 5. We examine the predicted

galaxy stellar masses (§ 5.1), and how this depends on both nu-

merical algorithms (§ 5.2) and feedback physics (§ 5.3), as well

as how it compares to previous theoretical work (§ 5.4). We show

that the treatment of feedback physics overwhelmingly dominates

these results, and discuss the distinct roles of multiple indepen-

0000 RAS, MNRAS

000

, 000–000

Hopkins et al.

dent feedback mechanisms. We also explore the predictions for the

Kennicutt-Schmidt relation (§ 5.5), the shape of galaxy star forma-

tion histories (§ 5.6), the star formation “main sequence” (§ 5.7),

and the “burstiness” of star formation (§ 5.8). We summarize our

important conclusions and discuss future work in § 6.

2 INITIAL CONDITIONS & GALAXY PROPERTIES

The simulations presented here are a series of fully cosmological

“zoom-in” simulations of galaxy formation; some images of the

gas and stars in representative stages are shown in Figs. 1-3. The

technique is well-studied; briefly, a large cosmological box is sim-

ulated at low resolution to

0, and then the mass within and

around halos of interest at that time is identified, traced back to the

starting redshift, and the Lagrangian region containing this mass is

re-initialized at much higher resolution for the ultimate simulation

(Porter 1985; Katz & White 1993).

We consider a series of systems with different masses. Table 1

describes the initial conditions. All simulations begin at redshifts

∼

100

−

125, with fluctuations evolved using perturbation theory

up to that point.

The specific halos we re-simulate are chosen to represent a

broad mass range and be “typical” in most properties (e.g. sizes,

formation times, and merger histories) relative to other halos of the

same

0 mass. The simulations

m09

and

m10

are constructed

using the methods from Onorbe et al. (2013); they are isolated

dwarfs. Simulations

m11

m12q

, and

m13

are chosen to match a

subset of initial conditions from the AGORA project (Kim et al.

2013), which will enable future comparisons with a wide range of

different codes. These are chosen to be somewhat quiescent merger

histories, but lie well within the typical scatter in such histories at

each mass (and each has several major mergers). Simulation

m12v

for contrast, is chosen to have a relatively violent merger history

(several major mergers since

∼

2), and is based on the initial con-

ditions studied in Kereš & Hernquist (2009) and Faucher-Giguère

& Kereš (2011).

In each case, the resolution is scaled with the simulated mass,

so as to achieve the optimal possible force and mass resolution. It is

correspondingly possible to resolve much smaller structures in the

low-mass galaxies. The critical point is that in all our simulations

with mass

, we resolve the Jeans mass/length of gas in

the galaxies, corresponding to the size/mass of massive molecu-

lar cloud complexes. This is necessary to resolve a genuine multi-

phase ISM and for our ISM feedback physics to be meaningful

(see Paper

for a detailed discussion of the scales that must be re-

solved for feedback to operate appropriately). Fortunately, most of

the mass and star formation in GMCs in both observations (Evans

1999; Blitz & Rosolowsky 2004) and simulated systems (Paper

)

is concentrated in the most massive GMCs, so resolution studies

in Paper

-Paper

confirm that resolving small molecular clouds

makes little difference. In terms of the Jeans mass/length of the

galaxies, our resolution is broadly comparable between different

simulations. Our worst resolution in units of the Jeans length/mass

occurs in the more massive galaxies at late times, when they are

relatively gas poor, and so (despite the large total galaxy mass) the

Jeans length can become relatively small.

Initial conditions were generated with the MUSIC code (Hahn & Abel

2011), using second-order Lagrangian perturbation theory.

The approximate Jeans (GMC) mass/length for the

0 disks, assuming

Toomre

∼

1, increases from

∼

(

∼

−

30 pc) in the

halos to

∼

(

∼

100

−

200 pc) in the

halos. If

1, or

Figure 4.

Galaxy stellar mass-halo mass relation at

Top

∗

(

halo

)

Bottom: M

∗

relative to the Universal baryon budget of the halo (

halo

Each simulation (points) from Table 1 is shown; large point denotes the

most massive halo in each box. We compare the relation if all baryons

became stars (

∗

halo

; dotted) and the observationally inferred re-

lationship as determined in Behroozi et al. (2012, magenta) and Moster

et al. (2013, cyan) (dashed lines denote extrapolation beyond the observed

range).

The agreement with observations is excellent at

halo

including dwarf though MW-mass galaxies. We stress that there are

zero

adjusted parameters here: stellar feedback, with known mechanisms taken

from stellar population models, is sufficient to explain galaxy stellar masses

at/below

∼

∗

We adopt a “standard” flat

CDM cosmology with

≈

Ω

−

Ω

Λ

≈

27, and

Ω

≈

046 for all runs.

3 BARYONIC PHYSICS

The simulations here use the physical models for star formation

and stellar feedback developed and presented in a series of papers

studying isolated galaxies (Hopkins et al. 2012c, 2013a, 2012b,a),

if the gas fractions are higher (at higher redshifts), the Jeans masses/lengths

are larger as well.

Because of our choice to match some of our ICs to widely-used examples

for numerical comparisons, they feature very small cosmological parameter

differences. These are percent-level, smaller than the observational uncer-

tainties in the relevant quantities (Planck Collaboration et al. 2013) and pro-

duce negligible effects compared to differences between randomly chosen

halos.

0000 RAS, MNRAS

000

, 000–000

Galaxies on FIRE: Feedback & Inefficient Star Formation

Table 1.

Simulation Initial Conditions

Name

halo

Merger

Notes

[

−

]

[

−

]

[

−

pc]

[

−

]

[

−

pc]

History

m09

3e9

1.8e2

0.7

8.93e2

normal

isolated dwarf

m10

1e10

1.8e2

1.4

8.93e2

normal

isolated dwarf

m11

1e11

5.0e3

5.0

2.46e4

quiescent

–

m12v

5e11

2.7e4

7.0

1.38e5

violent

several

2 mergers

m12q

1e12

2.0e4

7.0

1.97e5

normal

large (

∼

vir

) box

m13

1e13

2.0e5

1.5e6

150

normal

“small group” mass

Parameters describing the initial conditions for our simulations (units are physical):

(1)

Name: Simulation designation.

(2)

halo

: Approximate mass of the

0 “main” halo (most massive halo in the high-resolution region).

(3)

: Initial baryonic (gas and star) particle mass in the high-resolution region, in our highest-resolution

simulations.

(4)

: Baryonic force softening (minimum SPH smoothing lengths are comparable or smaller).

(3)

: Dark matter particle mass in the high-resolution region, in our highest-resolution simulations.

(4)

: Dark matter force softening (minimum softenings are fixed over the entire simulation duration).

Figure 5.

∗

−

halo

relation as Fig. 4, at different redshifts. Observational constraints are also shown at each redshift. With no tuned parameters, the

simulations predict

∗

−

halo

and, by extension, the stellar MF and galaxy clustering, at all

. Redshift evolution in

∗

−

halo

is weak, with the sense that

low-mass dwarfs become higher-

∗

, leading to a steeper faint-end galaxy MF, in agreement with constraints from reionization (see Kuhlen & Faucher-Giguère

2012, and references therein).

adapted for fully cosmological simulations. We summarize their

properties below, but refer to Appendix A for a more detailed ex-

planation and list of improvements. Readers interested in further

details (including resolution studies and a range of tests of the spe-

cific numerical methodology) should see Paper

& Paper

3.1 Cooling

Gas follows an ionized+atomic+molecular cooling curve from 10

−

K, including metallicity-dependent fine-structure and molecu-

lar cooling at low temperatures, and high-temperature (

metal-line cooling followed species-by-species for 11 separately

tracked species. At all times, we tabulate the appropriate ioniza-

tion states and cooling rates from a compilation of

CLOUDY

runs,

including the effect of the photo-ionizing background, accounting

for gas self-shielding. Photo-ionization and photo-electric heating

from local sources are accounted for as described below.

3.2 Star Formation

Star formation is allowed only in dense, molecular, self-gravitating

regions above

crit

(

crit

100 cm

−

for our primary runs, but

we also tested from

∼

−

1000 cm

−

). This threshold is much

higher than that adopted in most “zoom-in” simulations of galaxy

formation (the high value allows us to capture highly clustered star

formation). We follow Krumholz & Gnedin (2011) to calculate the

molecular fraction

in dense gas as a function of local column

density and metallicity, and allow SF only from molecular gas. We

0000 RAS, MNRAS

000

, 000–000

Hopkins et al.

also follow Hopkins et al. (2013d) and restrict star formation to

gas which is locally self-gravitating, i.e. has

≡

G m

gas

(

)

1 on the smallest available scale (

being our force soften-

ing or smoothing length). This forms stars at a rate

̇

∗

mol

(i.e. 100% efficiency per free-fall time); so that the galaxy and even

kpc-scale star formation efficiency is

not

set by hand, but regulated

by feedback (typically at much lower values). In Paper

, Paper

and Hopkins et al. (2013d) we show that the galaxy structure and

SFR are basically independent of the small-scale SF law, density

threshold (provided it is high), and treatment of molecular chem-

istry.

3.3 Stellar Feedback

Once stars form, their feedback effects are included from several

sources. Every star particle is treated as a single stellar popula-

tion, with a known age, metallicity, and mass. Then all feedback

quantities (the stellar luminosity, spectral shape, SNe rates, stellar

wind mechanical luminosities, metal yields, etc.) are tabulated as

a function of time directly from the stellar population models in

STARBURST99, assuming a Kroupa (2002) IMF.

(1)

Radiation Pressure:

Gas illuminated by stars feels a

momentum flux

̇

rad

≈

(

−

exp

(

−

optical

)) (

)

incident

along the optical depth gradient, where 1

+ Σ

gas

ac-

counts for the absorption of the initial UV/optical flux and multiple

scatterings of the re-emitted IR flux if the region between star and

gas particle is optically thick in the IR (see Appendix A). We as-

sume that the opacities scale linearly with gas metallicity.

(2)

Supernovae:

We tabulate the SNe Type-I and Type-II rates

from Mannucci et al. (2006) and STARBURST99, respectively, as

a function of age and metallicity for all star particles and stochas-

tically determine at each timestep if an individual SNe occurs. If

so, the appropriate mechanical luminosity and ejecta momentum is

injected as thermal energy and radial momentum in the gas within

a smoothing length of the star particle, along with the relevant mass

and metal yield (for all followed species). When the Sedov-Taylor

phase is not fully resolved, we account for the work done by hot

gas inside the unresolved cooling radius (converting the appropri-

ate fraction of the SNe energy into momentum). We discuss this in

detail in Appendix A, but emphasize that it is particularly important

in massive halos whose mass resolution

∼

is much larger

than the ejecta mass of a single SNe.

(3)

Stellar Winds:

Similarly, stellar winds are assumed to

shock locally and so we inject the appropriate tabulated mechan-

ical power

(

)

, wind momentum, mass, and metal yields, as a

continuous function of age and metallicity into the gas within a

smoothing length of the star particles. The integrated mass fraction

recycled in winds (including both fast winds from young stars and

slow AGB winds) and SNe is

∼

(4)

Photo-Ionization and Photo-Electric Heating:

Knowing

the ionizing photon flux from each star particle, we ionize each

neighboring neutral gas particle (provided there are sufficient pho-

tons, given the gas density, metallicity, and prior ionization state),

moving outwards until the photon budget is exhausted; this alters

the heating and cooling rates appropriately. The UV fluxes are also

There has been some debate in the literature regarding whether or not the

full

“boost” applies to the infrared radiation pressure when

1 (see

e.g. Krumholz & Thompson 2012, but also Kuiper et al. 2012). We have

considered alternatives, discussed in Paper I. However, in the simulations

here we never resolve the extremely high densities where

1 (where

this distinction is important), and so if anything are

under

-estimating the IR

radiation pressure.

used to determine photo-electric heating rates following Wolfire

et al. (1995).

Extensive numerical tests of the feedback models are pre-

sented in Paper

4 SIMULATION NUMERICAL DETAILS

All simulations are run using a newly developed version of

TreeSPH which we refer to as “

P-SPH

” (Hopkins 2013). This adopts

the Lagrangian “pressure-entropy” formulation of the SPH equa-

tions developed in Hopkins (2013); this eliminates the major differ-

ences between SPH, moving mesh, and grid (adaptive mesh) codes,

and resolves the well-known issues with fluid mixing instabilities

in previously-used forms of SPH (e.g. Agertz et al. 2007; Sijacki

et al. 2012). The gravity solver is a heavily modified version of

the

GADGET-3

code (Springel 2005); but

P-SPH

also includes sub-

stantial improvements in the artificial viscosity, entropy diffusion,

adaptive timestepping, smoothing kernel, and gravitational soften-

ing algorithm, as compared to the “previous generation.” These are

all described in detail in Appendix B.

We emphasize that our version of SPH has been tested exten-

sively and found to give good agreement with analytic solutions as

well as well-tested grid codes on a broad suite of test problems.

Many of these are presented in Hopkins (2013). This includes Sod

shock tubes; Sedov blastwaves; wind tunnel tests (radiative and

adiabatic, up to Mach

∼

); linear sound wave propagation; os-

cillating polytropes; hydrostatic equilibrium “deformation”/surface

tension tests (Saitoh & Makino 2013); Kelvin-Helmholtz and

Rayleigh-Taylor instabilities; the “blob test” (Agertz et al. 2007);

super-sonic and sub-sonic turbulence tests (from Mach

∼

−

); Keplerian gas ring, disk shear, and shearing shock tests

(Cullen & Dehnen 2010); the Evrard test; the Gresho-Chan vor-

tex; spherical collapse tests; and non-linear galaxy formation tests

such such as the Santa Barbara cluster comparison. Since it is crit-

ical for the problems addressed here that a code be able to handle

high dynamic range situations, the numerical method and parame-

ters such as SPH “neighbor number” were not modified for these

tests individually, but are similar to what we use in our production

runs in this paper.

5 RESULTS

5.1 Galaxy Masses as a Function of Redshift

Fig. 4 plots the

0 stellar mass-halo mass relation for our main

set of simulations from Table 1 (highest-resolution, with all physics

enabled). Note that although each high-resolution region at

contains one “primary” halo (the focus of that region), there are

several smaller-mass, independent halos also in that region. We

therefore identify and plot all such halos.

We exclude those ha-

los that are outside the high-resolution region (more than 1% mass-

We use the HOP halo finder (Eisenstein & Hut 1998) to automatically

identify halos (which combines an iterative overdensity identification with

a saddle density threshold criterion to merge subhalos and overlapping ha-

los). Halo masses are defined as the mass within a spherical aperture about

the density maximum with mean density

200 times the critical density at

each redshift (this is chosen to be similar to the choice used in abundance-

matching models, which define the observations to which we compare).

Stellar mass plotted is the total stellar mass within 20 kpc of the center of the

central galaxy in the halo (we do not include satellite galaxy masses). How-

ever, we have compared with the results of a basic friends-of-finds routine

or simple by-eye identification, and find that for the results here (focused

on simple, integral halo quantities, and ignoring subhalos), this makes no

0000 RAS, MNRAS

000

, 000–000

Galaxies on FIRE: Feedback & Inefficient Star Formation

contaminated by low-resolution particles; although varying this be-

tween 0

−

10% makes little difference to our comparisons here) or

insufficiently resolved (

01 times the primary halo mass, or with

dark matter particles). We also exclude subhalos/satellite

galaxies.

The known sources of stellar feedback we include, with

ad-

justment, automatically reproduce a relation between galaxy stel-

lar and halo mass consistent with the observations

from

halo

∼

−

. Despite the fact that this relation implies a non-

uniform (and even non-monotonic) efficiency of star formation as a

function of galaxy mass, we do

not

need to invoke different physics

or distinct parameters at different masses. This is particularly im-

pressive at low masses, where the integrated stellar mass must be

suppressed by factors of

∼

1000 relative to the Universal baryon

fraction. Unfortunately, at high masses (

), the large La-

grangian regions (hence large number of required particles) limit

the resolution we can achieve; we have experimented with some

low-resolution test runs which appear to produce overly massive

galaxies, but higher-resolution studies are required to determine if

that owes to a need for additional physics or simply poor numerical

resolution.

Interestingly, the scatter in

∗

at fixed

halo

appears to de-

crease with mass, from

∼

5 dex in dwarf galaxies (

halo

) to

∼

−

2 dex in massive (

∼

∗

)

galaxies. But given

the limited number of halos we study here, further investigation

allowing more diverse merger/growth histories is needed.

Fig. 5 shows the

∗

−

halo

relation at various redshifts. At

each

, we compare with observationally constrained estimates of

the

∗

−

halo

relation. Implicitly, if they agree in

∗

(

halo

)

, our

models are consistent with the observed stellar MF. At high red-

shifts, the halos we simulate are of course lower-mass, so eventu-

ally we have no high-mass galaxies; this limits the extent to which

our results can be compared to observations above

∼

5.2 (Lack of) Dependence on Numerical Methods

In Fig. 6 we investigate how the

∗

−

halo

relation depends on nu-

merical parameters and feedback. First we repeat Fig. 4 for simu-

lations with different numerical parameters. These can have signif-

icant quantitative effects, but do not qualitatively change our con-

clusions. Modest changes in resolution (our “low-resolution” runs

correspond to one power of two step in spatial resolution, and a

corresponding factor of 2

8 change in mass resolution) leads to

small changes in

∗

(generally we obtain larger

∗

at lower res-

olution, owing to artificially enhanced mixing and thus cooling of

diffuse gas since the ISM phases are less well-resolved, as well as

failure to resolve the clustering of star formation). At sufficiently

poor resolution (

100 times the particle masses used here) the re-

sults do diverge substantially (galaxy masses at

∼

0 are a factor

of several higher). We also varied the sizes of the high-resolution

Lagrangian regions of the halos; the results here are insensitive to

the region size if we choose sizes

vir

(at the redshift of inter-

est), but the cooling of halo gas and star formation are artificially

suppressed if the high-resolution region is much smaller.

We have also investigated different algorithmic methods

(varying the number of SPH particles to which energy and momen-

tum are coupled; the functional calculation of the local density and

significant difference. Likewise defining the mass within

∼

vir

instead

of 20 kpc makes no significant difference.

Note that Behroozi et al. (2012) and Moster et al. (2013) use definitions of

halo mass which differ slightly (by

≈

10%). For our purposes, this produces

negligible differences in our comparison.

Figure 6.

∗

−

halo

relation at

0, as Fig. 4.

Top:

Simulations with

different numerical parameters: we show the effects of varied resolution,

artificial viscosity, and the algorithmic implementation of feedback. We

also compare a completely different version of SPH (with a different set

of hydrodynamic equations), which is known to differ significantly in cer-

tain idealized hydrodynamic test problems. These have little effect on our

predictions.

Bottom:

Effect of physical variation in stellar feedback prop-

erties. We compare runs with no stellar feedback, with no supernovae (but

stellar winds, radiation pressure, and photo-ionization heating included), or

with no radiative feedback (radiation pressure and local HII-heating). “No

feedback” runs generally predict

∗

∼

halo

, in severe conflict with the

observations.

Removing radiative

SNe feedback also produce order-of-

magnitude too-large stellar masses. The non-linear

combination

of feedback

mechanisms (not any one in isolation) is critical to drive winds and regulate

galaxy masses.

pressure used for those couplings; where in the timestep the cou-

pling occurs; the use of adaptive versus fixed gravitational soften-

ings; equal versus distinct softenings for baryons and dark matter).

Subtle differences in the algorithmic implementation of feedback

have little

systematic

effect on the stellar mass, provided the same

mechanisms are included; however they can only be compared sta-

tistically, since the stochastic nature of feedback means that even

very subtle changes can produce significant differences in the exact

time history of bursts, for example. Of what we have considered,

the most important parameter is how we implement the momentum

gained during the Sedov-Taylor phase of SNe remnant expansion

when the cooling radius is unresolved (see Appendix A). Purely nu-

merical parameters such as the choice of softening (at fixed mass

0000 RAS, MNRAS

000

, 000–000

Hopkins et al.

resolution), and ensuring a sufficiently large number of SPH par-

ticles are sampled in the kernel so that local quantities (densities,

density gradients) needed to couple feedback appropriately are sta-

ble produce

∼

50% changes in the

0 stellar mass for Milky

Way-mass systems; a more detailed numerical study will be pre-

sented in future work. In Fig. 6 we see that in dwarf systems, these

can produce up to factor of several changes in the

0 stellar

mass, although the more typical variation is factor

∼

2. However

in all cases the results lie within the (rather large) range allowed by

observations.

In a companion paper, Kereš et al. (2013) consider the de-

tailed effects of substantial changes to each aspect of our numerical

method described in Appendix B. Here, we simply show a few ba-

sic comparisons. Considerable attention has recently been paid to

differences between the results of grid codes and older SPH meth-

ods (such as that in Springel & Hernquist 2002) for certain prob-

lems (especially sub-sonic fluid mixing instabilities; see Agertz

et al. 2007; Kitsionas et al. 2009; Bauer & Springel 2012; Vogels-

berger et al. 2012; Sijacki et al. 2012; Kereš et al. 2012). The nu-

merical method used for our standard simulations has been specif-

ically shown to resolve most of these discrepancies (giving results

quite similar to grid codes in test problems); this is verified in Hop-

kins (2013) for standard test problems and Kereš et al. (2013) for

cosmological simulations. However we have re-run some of our

simulations using the Springel & Hernquist (2002) formulation of

SPH (described in Appendix B), which shows the most pronounced

forms of these discrepancies. Despite the known differences be-

tween such methods for certain test problems, we find in Fig. 6

(

top

panel) that it makes little difference for the predicted galaxy

masses. The older SPH method gives slightly lower

∗

(

halo

)

(by about

≈

1 dex), primarily because cooling of hot gas is sup-

pressed by less-efficient mixing, but this is small compared to the

effects of including the appropriate stellar feedback physics.

It is important to stress that our conclusion here – that our

results depend only weakly on the numerical details – applies to

the galaxy stellar masses and other integrated quantities. It is less

clear that all properties of the simulations are so robust. This will

be studied in future work.

5.3 (Strong) Dependence on Feedback

The lower panel in Fig. 6 shows the effect of varying the physics of

feedback: we see dramatic differences in

∗

(

halo

)

. Removing all

feedback, gas cools and collapses on a dynamical time

dyn

within

the disk, forming stars at a rate

̇

∗

∼

gas

dyn

∼

̇

gas

where

̇

gas

is the inflow rate from the halo. Most of the baryons are turned into

stars.

If we turn off SNe feedback, but retain all other forms of feed-

back, the results are nearly as bad: again,

∗

is severely overpre-

dicted. In dwarfs, other forms of feedback may still be able to sup-

press SF significantly (giving

∗

halo

), but the masses are

still much too large relative to those observed, by factors

∼

−

50.

More interestingly, if we instead remove radiative feedback (photo-

heating and radiation pressure), but retain SNe (and stellar winds),

we see a nearly identical failure: while

∗

halo

, far too many

stars form.

Even with no feedback, at very low masses

halo

, some sup-

pression of SF occurs after reionization because we still include a photo-

ionizing background. However the predicted stellar mass is still larger than

observed by at least an order of magnitude.

Figure 7.

Comparison of the

∗

(

halo

)

relation (as Fig. 4) predicted by

other published simulations in the literature using sub-grid stellar feedback

models. We compile these results (colored points), where available, at low-

(

−

top

) and high-

(

−

bottom

). We compare against

the simulations presented here (gray points) with explicit feedback, and ob-

servational constraints (lines). Even sub-grid models which are “success-

ful” near

∼

∗

∼

0 over-predict

∗

(

halo

)

by an order-of-magnitude

relative to our explicit feedback simulations and observations at both low

masses (

halo

) and/or high redshifts (

2). The only excep-

tion appears to be the recent sub-grid model of Stinson et al. (2013), which

is explicitly adjusted to mimic the effects of radiative (“early”) feedback as

well as SNe, seen in our explicit feedback models.

5.4 Comparison to Previous Work

In Fig. 7, we compare our results (grey points) at low and high red-

shifts, to those from previous simulations spanning a wide range

of galaxy properties and numerical methods (Pelupessy et al. 2004;

Stinson et al. 2007, 2013; Mashchenko et al. 2008; Valcke et al.

2008; Governato et al. 2010, 2012; Oser et al. 2010; Brooks et al.

2011; Guedes et al. 2011; Sawala et al. 2011; Scannapieco et al.

2011; Okamoto 2013; Kannan et al. 2013).

All

of these simulations

include some form of sub-grid model designed to mimic the ulti-

mate effects of stellar feedback, although the prescriptions adopted

differ substantially between each. Most of these models are specif-

ically tuned to reproduce reasonable scaling for MW-mass sys-

tems at

∼

0. However, two discrepancies are immediately evi-

dent. First, nearly all the previous models predict much larger stel-

lar masses in dwarf galaxies with

halo

, compared

to either our simulations or the observational constraints. Second,

0000 RAS, MNRAS

000

, 000–000