TLA_large - MNRAS-2014-Hopkins-581-603.pdf

MNRAS

445,

581–603 (2014)

doi:10.1093/mnras/stu1738

Galaxies on FIRE (Feedback In Realistic Environments): stellar feedback

explains cosmologically inefficient star formation

Philip F. Hopkins,

1,2

‹

san Kere

Jos

norbe,

Claude-Andr

e Faucher-Gigu

ere,

2,5

Eliot Quataert,

Norman Murray

†

and 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, USA

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

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

Accepted 2014 August 22. Received 2014 August 21; in original form 2013 November 10

ABSTRACT

We present a series of high-resolution cosmological simulations

of galaxy formation to

spanning halo masses

∼

–10

, and stellar masses

∼

–10

. Our simulations

include fully explicit treatment of the multiphase interstellar medium and stellar feedback. The

stellar feedback inputs (energy, momentum, mass, and metal fluxes) are taken directly from

stellar population models. These sources of feedback, with

zero

adjusted parameters, reproduce

the observed relation between stellar and halo mass up to

halo

∼

. We predict weak

redshift evolution in the

∗

–

halo

relation, consistent with current constraints to

6. We

find that the

∗

–

halo

relation is insensitive to numerical details, but is sensitive to feedback

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

particularly in dwarf and high-redshift galaxies: radiative feedback (photoheating and radiation

pressure) is necessary to destroy giant molecular clouds and enable efficient coupling of later

supernovae to the gas. Star formation rates (SFRs) agree well with the observed Kennicutt

relation at all redshifts. The galaxy-averaged Kennicutt relation is very different from the

numerically imposed law for converting gas into stars, and is determined by self-regulation via

stellar feedback. Feedback reduces SFRs and produces reservoirs of gas that lead to rising late-

time star formation histories, significantly different from halo accretion histories. Feedback

also produces large short-time-scale variability in galactic SFRs, especially in dwarfs. These

properties are not captured by common ‘sub-grid’ wind models.

Key words:

stars: formation – galaxies: active – galaxies: evolution – galaxies: 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 time-scales of

∼

50 dynamical times (Kennicutt

1998), while the total fraction of giant molecular cloud (GMC)

mass converted into stars is only a few per cent (Zuckerman &

E-mail:

phopkins@caltech.edu

†

Canada Research Chair in Astrophysics.

Evans

1974; Williams & McKee

1997;Evans1999; Evans et al.

2009). Without strong stellar feedback, however, gas inside galax-

ies cools efficiently and collapses on a dynamical time, predicting

order-unity star formation efficiencies on all scales (Bournaud et al.

2010; Dobbs, Burkert & Pringle

2011; Hopkins, Quataert & Murray

2011; Harper-Clark & Murray

2011; Krumholz, Klein & McKee

2011;Tasker

2011).

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

mological models cools rapidly and inevitably turns into stars,

predicting galaxies with far larger masses than are observed

(e.g. Katz, Weinberg & Hernquist

1996; Somerville & Primack

1999;Coleetal.

2000; Springel & Hernquist

2003b;Kere

setal.

2009, and references therein). Decreasing the instantaneous star

formation efficiency does not eliminate this integral problem: the

2014 The Authors

Published by Oxford University Press on behalf of the Royal Astronomical Society

at California Institute of Technology on December 1, 2014

http://mnras.oxfordjournals.org/

Downloaded from

582

P. F. Hopkins et al.

amount of baryons in real galactic discs is much lower than that

predicted in models absent strong feedback (essentially, the Uni-

versal baryon budget; see White & Frenk

1991;Kere

setal.2009).

Constraints from intergalactic medium (IGM) enrichment require

that many of those baryons must have entered galaxy haloes and

discs at some point to be enriched, before being expelled (Aguirre

et al.

2001; Pettini et al.

2003; Songaila

2005;Martinetal.

2010).

Galactic superwinds 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

e 2006). Such

winds have been observed ubiquitously in local and high-redshift

star-forming galaxies (Martin

1999,

2006; Heckman et al.

2000;

Sato et al.

2009;Chenetal.2010; Steidel et al.

2010;Coiletal.

2011;Newmanetal.2012).

However, until recently, numerical simulations have been unable

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

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

galaxy mass or other properties), nor to simultaneously predict

the instantaneous inefficiency of star formation within galaxies.

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; Nagamine

2010; Brook et al.

2011; Bournaud et al.

2011; Powell, Slyz & Devriendt

2011, and references therein). More

recent simulations, using higher resolution and invoking stronger

feedback prescriptions, have seen strong winds, but have generally

found it necessary to include simplified prescriptions for ‘turning

off cooling’ in the SNe-heated gas and/or include some adjustable

parameters representing ‘pre-SNe’ feedback (see Governato et al.

2010; Macci

oetal.2012;Agertzetal.2013; Stinson et al.

2013;

Teyssier et al.

2013). This is physically motivated since feedback

processes other than SNe – protostellar jets, H

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 heating can have much larger effects (Evans et al.

2009;

Hopkins et al.

2011; Lopez et al.

2011;Tasker

2011; Stinson et al.

2013; Kannan et al.

2014).

And in fact, there have been many studies with enormously higher

resolution (enough to evolve each star explicitly) and a full treat-

ment of the radiation-magnetohydrodynamics and time dependence

of these multiple feedback mechanisms. Because of computational

limitations, however, these have necessarily been restricted to very

small systems, either single molecular clouds/star clusters (e.g.

Krumholz, Klein & McKee

2007;Offneretal.

2009,

2011; Harper-

Clark & Murray

2011; Krumholz et al.

2011;Bate

2012), or the ‘first

stars’ (e.g. Wise et al.

2012; Muratov et al.

2013; Pawlik, Milosavl-

jevi

c & Bromm

2013). But these studies,

without exception

,have

found that the non-linear interaction of the feedback mechanisms

above – especially the dual roles of H

photoionization and radia-

tion pressure in concert with SNe – is absolutely critical to explain

the generation of large local outflows, the self-regulation of star

formation, and the shape of the stellar initial mass function (IMF).

Despite these breakthroughs, given limited resolution and the

complexity of the baryonic physics, many cosmological models

have treated galactic wind generation and the inefficiency of star

formation in a tuneable, ‘sub-grid’ fashion. This is not to say that

the models have not tremendously 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 allowed feedback models, made predic-

tions for the critical role of outflows and recycling in enriching

the IGM, provided possible baryonic solutions to apparent dark

matter ‘problems’ (e.g. Pontzen & Governato

2012), demonstrated

the need for ‘early’ feedback from radiative mechanisms beyond

SNe alone, and generally created the framework for our interpre-

tation 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 particularly true for studies of

gas in the circumgalactic 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’ hydrodynamics or cooling, or

mimic strong feedback via pure thermal energy injection or ‘particle

kicks’ (see e.g. Hummels et al.

2013, for an explicit demonstration

of this).

Accurate treatment of star formation and galactic winds ulti-

mately requires realistic treatment of the stellar feedback processes

that maintain the multiphase interstellar medium (ISM). Motivated

by this philosophy (and building on the studies with single-star res-

olution), in Hopkins et al. (

2011, hereafter Paper

I) and Hopkins,

Quataert & Murray (

2012a, hereafter Paper

II), we developed a new

set of numerical models to follow stellar feedback on scales from

sub-GMC star-forming regions through galaxies. These simulations

include the energy, momentum, mass, and metal fluxes from stellar

radiation pressure, H

photoionization and photoelectric heating,

Types I and II SNe, and stellar winds (O-star and AGB). Critically,

the feedback is directly tied to the young stars, with the energet-

ics and time-dependence taken from stellar evolution models. In

our previous work, we showed, in isolated galaxy simulations, that

these mechanisms produce a quasi-steady ISM in which GMCs

form and disperse rapidly, with phase structure, turbulence, and

disc and GMC properties in good agreement with observations (for

various comparisons, see Hopkins, Kere

s & Murray

2013b; Hopkins

et al.

2012c,

2013d; Narayanan & Hopkins

2013). In Paper

I, Hop-

kins, Narayanan & Murray (2013c), and Hopkins et al. (2013a), we

showed that this leads naturally to ‘instantaneously’ inefficient star

formation (predicting the KS-law), regulated self-consistently by

feedback and

independent

of the numerical prescription for star for-

mation in very dense gas. In Hopkins, Quataert & Murray (

2012b,

hereafter Paper

III) and Hopkins et al. (2013e) 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 mechanisms is critical to produce

massive, multiphase 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 se-

ries, we present the FIRE (Feedback In Realistic Environments)

simulations:

a suite of fully cosmological ‘zoom-in’ simulations

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

the models and understand feedback in a wide range of environ-

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

opposed to focusing just on MW-like systems), and follow evolu-

tion fully to

0. Our suite of calculations includes several of the

highest resolution galaxy formation simulations that have been run

Movies and summaries of key simulation properties are available at

http://www.tapir.caltech.edu/

∼

phopkins/Site/Movies_cosmo.html

and the

FIRE project website:

http://fire.northwestern.edu

MNRAS

445,

581–603 (2014)

at California Institute of Technology on December 1, 2014

http://mnras.oxfordjournals.org/

Downloaded from

Galaxies on FIRE: feedback and star formation

583

Figure 1.

Gas in a representative simulation of a Milky Way-mass halo (m12i in Table ). 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 centred 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 disc has formed, with molecular gas

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

0. Our simulations utilize a significantly improved numerical

implementation 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

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

of feedback in determining galaxy structure, and many other open

questions.

In Sections 2–4, we describe our methodology. Section 2 de-

scribes the initial conditions for the simulations; Section 3 out-

lines the implementation of the key baryonic physics of cooling,

star formation, and feedback (a much more detailed description

is given in Appendix A); Section 4 briefly describes the improve-

ments in the numerical method compared to past work (again, more

details are in Appendix B). And in Appendix C, we test and com-

pare these algorithms with higher resolution simulations of isolated

(non-cosmological) galaxies.

We describe our results in Section 5. We examine the pre-

dicted galaxy stellar masses (Section 5.1), and how this depends

on both numerical algorithms (Section 5.3) and feedback physics

(Section 5.4), as well as how it compares to previous theoretical

work (Section 5.5). We show that the treatment of feedback physics

overwhelmingly dominates these results, and discuss the distinct

roles of multiple independent feedback mechanisms. We also ex-

plore the predictions for the KS relation (Section 5.6), the shape

of galaxy SFHs (Section 5.7), the star formation ‘main sequence’

(Section 5.8), and the ‘burstiness’ of star formation (Section 5.9).

We summarize our important conclusions and discuss future work in

Section 6.

2 INITIAL CONDITIONS AND 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 haloes 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 (with gas added) 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 haloes 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 haloes of the same

0 mass. The simulations m09 and m10 are constructed using

the methods from Onorbe et al. (

2014); they are isolated dwarfs.

Simulations m11, m12q, m12i, and m13 are chosen to match a

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

Both gas and stellar images are true three-colour volume renderings gen-

erated by ray-tracing lines of sight through the simulation (with every gas or

star particle a source, respectively). For the stars, the physical luminosities

and dust opacities in each band are used to generate the observed intensity

map. For the gas, we construct synthetic ‘bands’ where the particle emis-

sivity is uniform if it falls within the temperature range specified, and zero

otherwise, and the particle opacity is uniform across bands.

Initial conditions were generated with the

MUSIC

code (Hahn & Abel

2011

using second-order Lagrangian perturbation theory.

MNRAS

445,

581–603 (2014)

at California Institute of Technology on December 1, 2014

http://mnras.oxfordjournals.org/

Downloaded from

584

P. F. Hopkins et al.

Figure 2.

Stars in the m12i simulation at

∼

0, in a box 50 kpc on a side near

present-time. Image is a mock

composite. The disc is approximately

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

STARBURST

99 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

), atten-

uating with a MW-like reddening curve with constant dust-to-metals ratio

for the abundances at each point.)

2014), 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 conditions studied in Kere

s & Hernquist (2009) and Faucher-

Gigu

ere & Kere

s(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 molecular cloud

complexes. This is necessary to resolve a genuine multiphase ISM

and for our ISM feedback physics to be meaningful. Fortunately,

because most of the mass and star formation in GMCs in both

observations (Evans

1999; Blitz & Rosolowsky

2005) and simulated

systems (Paper

II) is concentrated in the most massive GMCs, the

resolution studies in Paper

I and Paper

II confirm that resolving

small molecular clouds makes little difference. We refer interested

readers to Paper

II for a detailed discussion of the scales that must

be resolved for feedback to operate appropriately, but note here that

all our simulations are designed to be approximately comparable to

the ‘high-resolution’ simulations of isolated galaxies and the ISM

in Paper

I and Paper

II, within the range of resolution where the

results in those studies (SFRs, wind outflow rates, GMC lifetimes,

etc.) were numerically converged (unfortunately, it is not possible

to evolve cosmological simulations to

0 with the ‘ultrahigh’

sub-pc resolution therein).

In terms of the Jeans mass/length of the galaxies, our resolution

is broadly comparable between different simulations. Our worst

Figure 3.

Gas, as Fig.

, for a dwarf galaxy (m10 in Table

). Top: 40 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.

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.

Every galaxy identified in this paper contains at

least

bound particles.

The approximate Jeans (GMC) mass/length for the

0 discs, assuming

Toomre

∼

1, increases from

∼

(

∼

10–30 pc) in the

haloes to

∼

(

∼

100–200 pc) in the

haloes. If

1, or

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

are larger as well.

MNRAS

445,

581–603 (2014)

at California Institute of Technology on December 1, 2014

http://mnras.oxfordjournals.org/

Downloaded from

Galaxies on FIRE: feedback and star formation

585

Table 1.

Simulation initial conditions.

Name

halo

Merger

Notes

(

−

)(

−

)(

−

pc) (

−

)(

−

pc)

history

m09

1.9e9

1.8e2

1.0

8.93e2

Normal

Isolated dwarf

m10

0.8e10

1.8e2

2.0

8.93e2

Normal

Isolated dwarf

m11

1e11

5.0e3

5.0

2.46e4

Quiescent –

m12v

5e11

2.7e4

7.0

1.38e5

100

Violent

Several

2mergers

m12q

1e12

5.0e3

7.0

1.97e5

100

Late merger –

m12i

1e12

3.5e4

1.97e5

100

Normal

Large (

∼

vir

) box

m13

1e13

2.6e5

1.58e6

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)

: minimum baryonic gravity/force softening (minimum SPH smoothing lengths are comparable or smaller).

Recall, force softenings are adaptive (mass resolution is fixed); for more details see Appendix B.

(3)

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

(4)

: minimum dark matter force softening (fixed in physical units at all redshifts).

We adopt a ‘standard’ flat

CDM cosmology with

≈

0.7,

−

≈

0.27, and

≈

0.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 pa-

pers studying isolated galaxies (Paper

III; Hopkins et al.

2012c,

2013a,d), adapted for fully cosmological simulations. We summa-

rize their properties below, but refer to Appendix A for a more

detailed explanation and list of improvements. Readers interested

in further details (including resolution studies and a range of tests of

the specific numerical methodology) should see Paper

I and Paper

II.

3.1 Cooling

Gas follows an ionized

atomic

molecular cooling curve from 10

to 10

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

ular 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 ionization

states and cooling rates from a compilation of

CLOUDY

runs, includ-

ing the effect of the photoionizing background, accounting for gas

self-shielding. Photoionization and photoelectric 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

∼

10 to 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 calcu-

late the molecular fraction

in dense gas as a function of lo-

cal column density and metallicity, and allow star formation only

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

2013

) and produce

negligible effects compared to differences between randomly chosen haloes.

from molecular gas. We also follow Hopkins et al. (

2013c)and

restrict star formation to gas which is locally self-gravitating, i.e.

has

≡

δv

gas

(

<δ

)

1 on the smallest available scale

(

being our force softening or smoothing length). This forms stars

at a rate ̇

∗

mol

(i.e. 100 per cent 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). Because of this, in Paper

I, Paper

II, and Hopkins et al.

(2013c) we show that the galaxy structure and SFR are basically

independent of the small-scale star formation law, density threshold

(provided it is high), and treatment of molecular chemistry.

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

STARBURST

99, assuming a Kroupa (2002)IMF.

(1)

Radiation pressure:

gas illuminated by stars feels a momen-

tum flux

rad

≈

−

exp (

−

optical

)) (1

)

incident

along

the optical depth gradient, where 1

gas

accounts

for the absorption of the initial UV/optical flux and multiple scat-

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

particle is optically thick in the IR (see Appendix A). We assume

that the opacities scale linearly with gas metallicity.

(2)

Supernovae:

we tabulate the SNe Type I and Type II rates

from Mannucci, Della Valle & Panagia (2006)and

STARBURST

99,

respectively, as a function of age and metallicity for all star parti-

cles and stochastically determine at each timestep if an individual

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

and Davis et al.

2014

, who find much stronger effects in the infrared). We have considered

alternatives, discussed in Paper

. However, in the simulations here we never

resolve the extremely high densities where

1 (where this distinction

is important), and so if anything we

under

estimate the IR radiation pressure,

even compared to the most conservative studies.

MNRAS

445,

581–603 (2014)

at California Institute of Technology on December 1, 2014

http://mnras.oxfordjournals.org/

Downloaded from

586

P. F. Hopkins et al.

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 appropriate fraction of the SNe energy into momentum). We

discuss this in detail in Appendix A, but emphasize that it is partic-

ularly important that SNe momentum not be neglected in massive

haloes 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 mechanical 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

∼

0.3.

(4)

Photoionization and photoelectric heating:

knowing the ion-

izing photon flux from each star particle, we ionize each neighbour-

ing neutral gas particle (provided there are sufficient photons, given

the gas density, metallicity, and prior ionization state), moving out-

wards until the photon budget is exhausted; this alters the heating

and cooling rates appropriately. The UV fluxes are also used to de-

termine photoelectric heating rates following Wolfire et al. (

1995).

Extensive numerical tests of the feedback models are presented

in Paper

II.

4 SIMULATION NUMERICAL DETAILS

All simulations are run using a newly developed version of TreeSPH

which we refer to as ‘

SPH

’ (Hopkins

2013), in the code

GIZMO

This adopts the Lagrangian ‘pressure-entropy’ formulation of the

smoothed particle hydrodynamics (SPH) equations developed in

Hopkins (2013); this eliminates the major differences 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;Sijackietal.

2012). The

gravity solver is a heavily modified version of the

GADGET

-3 code

(Springel

2005); but

GIZMO

also includes substantial improvements

in the artificial viscosity, entropy diffusion, adaptive timestepping,

smoothing kernel, and gravitational softening algorithm, as com-

pared to the ‘previous generation’. These are all described in detail

in Appendix B.

We emphasize that our version of SPH has been tested extensively

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; oscillating poly-

tropes; hydrostatic equilibrium ‘deformation’/surface tension tests

(Saitoh & Makino

2013); Kelvin–Helmholtz and Rayleigh–Taylor

instabilities; the ‘blob test’ (Agertz et al.

2007); supersonic and

sub-sonic turbulence tests (from Mach

∼

0.1–10

); Keplerian gas

ring, disc shear, and shearing shock tests (Cullen & Dehnen

2010);

the Evrard test; the Gresho–Chan vortex; spherical collapse tests;

and non-linear galaxy formation tests. For additional tests showing

Details of the

GIZMO

code, together with a limited public version,

user’s guide, movies, and test problem examples, are available at

www.tapir.caltech.edu/

∼

phopkins/Site/GIZMO.html

the improvements relative to previous-generation SPH, see Hu et al.

(2014). Since it is critical for the problems addressed here that a

code be able to handle high dynamic range situations, the numerical

method and parameters such as SPH ‘neighbour number’ were not

modified for these tests individually, but are similar to what we use

in our production runs in this paper.

In Appendix B, we note that we have explicitly tested many of

the purely numerical elements of the gravity and hydrodynamic

solvers in the simulations shown here: for example, whether to

use adaptive or fixed gravitational softenings, the choice of SPH

smoothing kernel, and the timestepping algorithm. However, we do

not discuss these in the main text because they produce extremely

small (

10 per cent-level) differences in the quantities plotted 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 sev-

eral smaller-mass, independent haloes also in that region. We there-

fore identify and plot all such haloes.

We exclude those haloes

that are outside the high-resolution region (more than 1 per cent

mass-contaminated by low-resolution particles; although varying

this between 0.5 and 10 per cent makes little difference to our com-

parisons here) or insufficiently resolved (

0.01 times the primary

halo mass, or with

dark matter particles). We also exclude

sub-haloes/satellite galaxies.

The known sources of stellar feedback we include, with

ad-

justment, automatically reproduce a relation between galaxy stellar

and halo mass consistent with the observations

from

halo

∼

–

. Specifically, the distribution of points for all

halo

≤

is statistically consistent (in a

sense) with having been

drawn from the Moster et al. (2013) curve; if we allow for the

observed scatter (

∼

0.15 dex at

∼

, the width between the

We use the

HOP

halo finder (Eisenstein & Hut

1998

) to automatically

identify haloes (which combines an iterative overdensity identification with

a saddle density threshold criterion to merge sub-haloes and overlapping

haloes). 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 centre of the central galaxy in the halo (we do not include satellite

galaxy masses). However, we have compared with the results of a basic

friends-of-friends routine or simple by-eye identification, and find that for

the results here (focused on simple, integral halo quantities, and ignoring

sub-haloes), this makes no significant difference. Likewise defining the mass

within

∼

0.1

vir

instead of 20 kpc makes no significant difference.

Note that Behroozi et al. (

2013

) and Moster et al. (

2013

) use definitions

of halo mass which differ slightly (by

≈

10 per cent). For our purposes, this

produces negligible differences in our comparison.

We can also fit the points here to the same power-law functional form

used empirically: if we do so, the best-fitting slope and normalization are

both within 0.5

of the fit to observations in Moster et al. (

2013

)–the

error bar is dominated by the small-number statistics in our halo sampling.

The simulations with

halo

are statistically inconsistent with

the extrapolation of the flatter slope from Behroozi et al. (

2013

), but this is

entirely below the region actually observed, where Behroozi et al. (

2013

)and

MNRAS

445,

581–603 (2014)

at California Institute of Technology on December 1, 2014

http://mnras.oxfordjournals.org/

Downloaded from

Galaxies on FIRE: feedback and star formation

587

Figure 4.

Galaxy stellar mass-halo mass relation at

0. Top:

∗

(

halo

Bottom:

∗

relative to the Universal baryon budget of the halo (

halo

Each simulation (points) from Table

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 relationship

as determined in Moster, Naab & White (

2013

, magenta) and Behroozi,

Wechsler & Conroy (

2013

, cyan) – dashed lines denote extrapolation beyond

the observed range; see footnote 9. 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

∼

∗

plotted lines, increasing to

∼

0.3 dex at the lowest observed masses),

then all our primary galaxies lie within the 2

scatter.

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 pa-

rameters at different masses. This is particularly impressive at low

masses, where the integrated stellar mass must be suppressed by

factors of

∼

1000 relative to the Universal baryon fraction. Unfor-

tunately, at high masses (

), the large Lagrangian regions

(hence large number of required particles) limit the resolution we

can achieve; we have experimented with some low-resolution test

Moster et al. (

2013

) agree well, and there the simulations do not significantly

‘prefer’ either fit.

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

may decrease weakly

with mass, from

∼

0.5 dex in dwarf galaxies (

halo

)

∼

0.1–0.2 dex in massive (

∼

∗

) galaxies. But given the limited

number of haloes 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 (given, of course, the limited

statistics by which we are ‘sampling’ the MF). At high redshifts, the

haloes we simulate are of course lower mass, so eventually we have

no high-mass galaxies; this limits the extent to which our results

can be compared to observations above

∼

5.2 Other (basic) galaxy properties

We wish to focus here on galaxy masses and star formation histories

(SFHs). Companion papers (in preparation) will examine the galaxy

morphologies and other observables in more detail. It is important,

when studying those properties, to construct a meaningful compar-

ison (e.g. using the same methods and wavelengths observed), and

this is non-trivial. Moreover, it is by no means obvious that these

properties are as robust to numerical details as the galaxy stellar

masses (discussed further below), and it is completely outside the

scope of this paper to fairly explore those dependences.

That said, we can briefly note the basic properties of the specific

simulations in Table

1 at

0, with the caveat that these

may

not be robust

to changes in either the initial conditions (the partic-

ular halo simulated) or our numerical methods. Morphologically,

0, run m09 resembles an ultrafaint dwarf; m10 a thick, but

rotating dwarf irregular; and m11 a more ‘fluffy’ dwarf spheroidal.

Runs m12v, m12q, m12i produce bulge

disc systems, with m12v

showing a prominent bulge at all times

2; m12q is more disc-

dominated until a late major merger at

0.5 destroys the disc;

and run m12i produces a stellar disc with little bulge. Run m13

is totally bulge dominated. Each galaxy has an approximately flat

rotation curve outside of the central couple kpc; those with

∗

slowly rise with radius to

max

, and the more massive sys-

tems are flat to within the central

∼

kpc, except for m12v (where

the compact bulge leads to a central-kpc spike at

∼

250 km s

−

The galaxy sizes, measured as the half-stellar mass effective radii,

are (0.3, 0.52, 3.5, 2.8, 3.2, 4.2, 4.8) kpc for (m09, m10, m11, m12v,

m12q, m12i, m13), consistent with the observed stellar size–mass

relation (Shen et al.

2003;Wolfetal.

2010).

5.3 (Lack of) dependence on numerical methods

In Fig.

6, we investigate how the

∗

–

halo

relation depends on

numerical parameters and feedback. First, we repeat Fig.

4 for

simulations with different

purely numerical

parameters. These can

and do, indeed, have significant quantitative effects – they can easily

shift the predicted stellar masses by factors

∼

2–3. However, we

stress that they do not

qualitatively

change our conclusions.

Modest changes in resolution (our ‘low-resolution’ runs corre-

spond to one power of two step in spatial resolution, and a cor-

responding factor of 2

8 change in mass resolution) lead to

significant, but not order-of-magnitude, changes in

∗

: generally

we obtain larger

∗

by factors of

∼

1.5 at high masses (

halo

)and

∼

2–3 at the lowest masses (

halo

)at

MNRAS

445,

581–603 (2014)

at California Institute of Technology on December 1, 2014

http://mnras.oxfordjournals.org/

Downloaded from

588

P. F. Hopkins et al.

Figure 5.

∗

–

halo

relation as Fig.

(points follow the legend therein), at different redshifts. Observational constraints are also shown at each redshift (each

pair of lines shows the

fit to the observations at that 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

ere

2012

, and references therein).

lower resolution, owing to a combination of (a) artificially enhanced

mixing and thus cooling of diffuse gas, since ISM phases are less

well-resolved, and (b) the fact that the coupling of feedback energy

and momentum is necessarily spread over larger mass elements.

If we downgrade our resolution more substantially – by a factor

∼

100 in mass, or

10 in spatial scale (i.e. using the

100 pc

spatial resolution which is typical of most previous cosmological

simulations), the results diverge more substantially: galaxy masses

∼

0 are a factor of

∼

3–5 higher at high masses and

∼

10 higher

at low masses. This makes sense, because at that resolution, we sim-

ply cannot meaningfully resolve even the most massive structures

in the ISM.

Some of our numerical tests are not plotted here because their

effects are not significant. We have, for example, re-run several

simulations with twice and five times larger dark matter soften-

ing lengths (same baryonic softening); using or de-activating adap-

tive gravitational softenings (which ensure there are always

∼

100

neighbour particles in the softening kernel); varying the number

of SPH ‘neighbours’ in the hydrodynamic kernel and number of

SPH particles to which energy and momentum are coupled; using

a single timestep or Strang-split integration scheme in the code;

varying the Courant factor of the hydrodynamic solver; changing

the order of operator-splitting for the cooling and feedback steps;

or forcing equal versus allowing separate gravitational softenings

We have run a couple tests with 30 times higher particle numbers than our

production-quality runs (for m12i and m10), to

2, and found that the

stellar masses at this time and earlier vary by

∼

10–50 per cent from those

quoted here. However, this appears to be primarily stochastic, rather than

systematic, so we suspect the masses will not change much further at still

higher resolution.

for baryons and dark matter. These produce very small (

10 per

cent) differences. We also varied the sizes of the high-resolution

‘zoom-in’ Lagrangian regions of the haloes; the results here are

insensitive to the region size if we choose sizes

vir

(at the

redshift of interest), but the cooling of halo gas and star forma-

tion are artificially suppressed if the high-resolution region is much

smaller.

Changing the small-scale star formation prescription in the

simulations has very little effect on our predictions. This is ex-

pected based on all of our previous studies using isolated (non-

cosmological) simulations (for explicit examples where we vary the

density threshold and instantaneous ‘efficiency’ of star formation in

dense gas by factors of

1000, as well as the density, temperature,

chemical, and virial-state dependence of star formation, see Paper

Paper

II, and Hopkins et al.

2013c). Globally, star formation is

feedback-regulated

: a certain number of young stars are required to

balance gravitational collapse/dissipation, independent of how those

stars form. So long as cooling can proceed, they will form (what will

change, via this self-regulation, if we change e.g. the density thresh-

old above which stars form, is the amount of gas which ‘builds up’

above that threshold; see Hopkins et al.

2013d). In Fig.

6,weshow

examples where we vary the density threshold for star formation

by factors of

∼

100 (producing

0.2 dex random/non-systematic

changes in stellar mass), or impose a much stricter local virial

criterion for star formation (local virial parameter

0.5 instead of

1; producing an

∼

20 per cent difference in stellar mass); we have

also experimented with removing the virial parameter entirely or

dividing the instantaneous efficiency of star formation in dense gas

by a factor of 100 (both produce

10 per cent changes).

We have also investigated different purely algorithmic methods

for coupling the same feedback physics. Subtle differences in the

algorithmic implementation of feedback have little

systematic

effect

MNRAS

445,

581–603 (2014)

at California Institute of Technology on December 1, 2014

http://mnras.oxfordjournals.org/

Downloaded from

Galaxies on FIRE: feedback and star formation

589

Figure 6.

∗

–

halo

relation at

0, as Fig.

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

idealized hydrodynamic test problems. These have little effect on our pre-

dictions. Bottom: effect of physical variation in stellar feedback properties.

We compare runs with no stellar feedback, with no SNe (but stellar winds,

radiation pressure, and photoionization heating included), or with no ra-

diative feedback (radiation pressure and local H

-heating). ‘No feedback’

runs generally predict

∗

∼

halo

, in severe conflict with the observa-

tions; see footnote 12. Removing radiative or 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.

on the stellar mass, provided the same mechanisms are included;

however, they can only be compared statistically, since the stochastic

nature of feedback means that even very subtle changes can produce

significant differences in the exact time history of bursts, for exam-

ple. 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 un-

resolved (see Appendix A). For example, one of the ‘mod. SNe

coupling algorithm’ examples changes the particle weights (using

a standard SPH kernel weight – effectively mass-weighted in the

smoothing kernel – instead of a volumetric weighting) used to de-

termine the coupling of SNe energy and momentum in the kernel.

This can have dramatic effects on test problems: for a SNe in an

infinitely thin, adiabatic disc with a low-density exterior, a mass-

weighting couples all the momentum in the disc plane, instead of

the vertical direction (the correct solution). Nevertheless, we see

that this has relatively weak (

∼

20 per cent) effects on the stellar

mass and SFH (in part because, in the average over many SNe over

large volumes, all that matters is the total feedback input); however,

it can significantly affect the morphological structure of e.g. the

dense gas in a thin disc. We have also experimented with differ-

ent functional forms for the ratio of the SNe energy and momentum

coupling (producing small effects). The ‘mod. RP algorithm’ choice

discretizes our radiation pressure term (which is usually a continu-

ous force) into intentionally very large (

500 km s

−

) ‘kicks’ (this

keeps the same total momentum flux, but makes each such particle

‘kicked’ unbound) – unsurprisingly this suppresses star formation

further, but only by a factor of

∼

3. In the ‘mod. RP

SNe’ choice,

we discretize the radiation pressure into smaller kicks (

5kms

−

)

and see this has little effect (as expected). In all cases, the results

lie within the (rather large) range allowed by observations.

In a companion paper, Kere

s et al. (in preparation) consider the

detailed effects of substantial changes to each aspect of our numer-

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

basic comparisons. Considerable attention has recently been paid to

differences between the results of grid codes and older SPH methods

(such as that in Springel & Hernquist

2002) for certain problems

(especially sub-sonic fluid mixing instabilities; see Agertz et al.

2007; Kitsionas et al.

2009; Bauer & Springel

2012;Kere

setal.

2012;Sijackietal.2012; Vogelsberger et al.

2012). The numeri-

cal method used for our standard simulations has been specifically

shown to resolve most of these discrepancies (giving results quite

similar to grid codes in test problems); this is verified in Hopkins

(2013) for standard test problems and Kere

s et al. (in preparation)

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.

(

top

panel) that it makes little difference for the predicted galaxy

masses. The older SPH method gives slightly lower

∗

(

halo

)(by

about

≈

0.15 dex), primarily because cooling of diffuse ‘hot halo’

gas is suppressed by less-efficient mixing. But for this specific ques-

tion, the effect is quite small compared to the effects of including the

appropriate stellar feedback physics. We also show an experiment

where we adopt an entirely distinct artificial viscosity prescription

(see Appendix B for details), which produces negligible differences.

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 lowest order, integrated quantities. In future

work, we will study other properties of the simulations, such as the

galaxy morphologies, which can (and do) depend on some param-

eters much more sensitively. For example, the modifications to the

SNe coupling algorithm described above, which produce very little

systematic change in our predicted stellar masses and SFHs, pro-

duce surprisingly large changes to the angular momentum content

and thickness of discs in the more massive galaxies.

Finally, we show these results in part to stress, emphatically, that

while there are

always

numerical choices in any code, there has

been no ‘tuning’ of these parameters for our study here. Certainly

none of these has been ‘fit’ or ‘adjusted’ to match any observations,

and all the choices above are held constant across our standard set

of simulations, using values calibrated from simple test problems

(e.g. Hopkins

2013).

MNRAS

445,

581–603 (2014)

at California Institute of Technology on December 1, 2014

http://mnras.oxfordjournals.org/

Downloaded from

590

P. F. Hopkins et al.

5.4 (Strong) dependence on feedback

The lower panel in Fig.

6 shows the effect of varying the

physics

feedback: now, we see dramatic differences in

∗

(

halo

). Removing

all feedback (every mechanism listed in Section 3.3), gas cools and

collapses on a dynamical time

dyn

within the disc, forming stars at

arate

∗

∼

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 feedback,

the results are nearly as bad: again,

∗

is severely overpredicted

in both dwarfs and Milky Way (MW)-mass systems. In the

haloes, with no SNe, other forms of feedback may still suppress

star formation significantly (so

∗

halo

), but the masses are

still much too large relative to those observed by factors of

∼

100.

We also note that, as many previous studies have pointed out (Mur-

ray, Quataert & Thompson

2005; McKee & Ostriker

2007; Shetty

&Ostriker2008; Faucher-Gigu

ere, Quataert & Hopkins

2013), it

is ultimately the

momentum

injected by SNe, not just the thermal

energy, which regulates star formation. So as expected, if we artifi-

cially turn off the SNe momentum (coupling only thermal energy,

as is common in many cosmological simulations), then in our sim-

ulations of massive (

) haloes, this is nearly as bad as

removing SNe entirely. In the lowest mass dwarfs, the discrepancy

is not so severe (factor

2 changes in the SFH), because the mass

resolution (

∼

100 M

) is such that the early expansion phases of

SN remnants (in which the thermal energy begins to be converted

into momentum) are well resolved.

If we remove radiative feedback entirely (both radiation pres-

sure and local photoionization and photoelectric heating, as de-

scribed in Section 3.3), but retain SNe (and stellar winds), we see

a nearly identical failure (to the no-SNe case) in both dwarfs and

massive galaxies: while

∗

halo

, far too many stars form. As

we showed in Paper

II, these mechanisms are critical to disrupt the

dense regions of GMCs in which young stars are born,

before

SNe

explode, and thus allowing the SNe to heat larger, lower density

volumes of gas (which can both avoid overcooling and feel the

collective effects of many SNe rather than just one), and therefore

actually generate significant galactic outflows. The same result is

found (on smaller scales) in much higher resolution simulations

of either single star clusters or the first stars, which directly treat

the radiation-hydrodynamics with each single star as a source (e.g.

Offner et al.

2009; Krumholz et al.

2011;Tasker

2011;Wiseetal.

2012).

Interestingly, in the dwarfs, if we turn off

only

radiation pressure,

only

photoionization heating, the effect is much less severe: the

predicted stellar mass is still significantly larger, but it is

100 times

larger when both are removed. Radiation pressure can, to some

extent, ‘make up for’ the loss of photoheating, and vice versa.

This should actually not be surprising: the most massive GMCs

in dwarf galaxies have local characteristic velocities

10 km s

−

thus either H

heating or UV radiation pressure alone can disrupt

them (though we expect, under these conditions, H

heating should

dominate), and this is completely consistent with both observa-

tions of star-forming regions (e.g. Lopez et al.

2011) and numerical

radiation-hydrodynamic simulations of low-density, low-velocity

clouds (Harper-Clark & Murray

2011; Sales et al.

2014). And in-

deed this tradeoff between photoheating and radiation pressure in

Even with no feedback, at very low masses

halo

,some

suppression of star formation occurs after reionization because we still

include a photoionizing background. However, the predicted stellar mass is

still larger than observed by at least an order of magnitude.

small clouds is exactly what we saw in our ultrahigh resolution sim-

ulations of isolated dwarfs of the same mass (discussed extensively

in Paper

II; see figs 7, 9, and 14–19 therein).

In the massive systems, on the other hand, the radiation pressure

term becomes more important than the H

heating. We see this

in tests with both m12q and m12v. Even when the difference in

stellar mass is not large (e.g. the m12v case), the lack of radiation

pressure feedback is particularly evident in the dense, early-forming

centre of the galaxy, where in the runs without radiation pressure

feedback an enormous central density ‘spike’ appears, leading to a

very large circular velocity of

∼

400 km s

−

in the central regions of

these systems. At these densities, H

photoheating is dynamically

insignificant.

If we disable stellar wind feedback (specifically, retaining stellar

winds as a source of mass and metals, but associating no energy or

momentum with that mass), and retain all other feedback, we see

relatively weak effects. This is not surprising: their momentum flux

is comparable to but not larger than other sources, and their ener-

getics are much less than SNe. But they are obviously an extremely

important source of mass and metals in the ISM.

5.5 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, van der Werf

&Icke

2004; Stinson et al.

2007,

2013; Mashchenko, Wadsley &

Couchman

2008; Valcke, de Rijcke & Dejonghe

2008; Feldmann

et al.

2010;Governatoetal.2010,

2012;Oseretal.2010; Brooks

et al.

2011; Guedes et al.

2011;Sawalaetal.2011; Scannapieco et al.

2011; De Rossi et al.

2013; Okamoto

2013; Kannan et al.

2014).

All

of these simulations include some form of sub-grid model designed

to mimic the ultimate effects of stellar feedback, although the pre-

scriptions adopted differ substantially between each. Most of these

models are specifically tuned to reproduce reasonable scaling for

MW-mass systems at

∼

0. However, two discrepancies are imme-

diately evident. First, nearly all the previous models predict much

larger stellar masses in dwarf galaxies with

halo

11.5

compared to either our simulations or the observational constraints.

Secondly, even simulations which produce excellent agreement with

the observations at

0 tend to predict far too much star forma-

tion at high redshift (take e.g. the simulation in Guedes et al.

2011,

which produces a MW-like system with many properties consistent

with observations at

0, but has turned nearly all its baryons into

stars at

2).

These are similar to the discrepancies that appear when we re-run

our simulations excluding radiative feedback. And indeed, nearly

all of the models from the literature in Fig.

7, even given various

freely adjustable parameters, are designed and motivated only to

reproduce the effects of SN feedback, which we have shown is

insufficient to explain the observations.

In fact, the only sub-grid models, to our knowledge, which cur-

rently do not produce such discrepancies (and agree broadly with

our simulations both at low masses and high redshifts) are the recent

generation of models in Stinson et al. (

2013), Aumer et al. (

2013),

Ceverino et al. (2014), and Trujillo-Gomez et al. (

2013); for some

additional results from these see Kannan et al.

2014. These new

models (all of which have been developed recently) are

specifically

designed/tuned to mimic the effects of radiative feedback (albeit in-

directly), and to reproduce via simple sub-resolution prescriptions

(including turning off cooling) some of the most important effects

of radiation pressure and photoheating which were studied in our

MNRAS

445,

581–603 (2014)

at California Institute of Technology on December 1, 2014

http://mnras.oxfordjournals.org/

Downloaded from

Galaxies on FIRE: feedback and star formation

591

Figure 7.

Comparison of the

∗

(

halo

) relation (as Fig.

) predicted by

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

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

(

0–0.5; top) and high-

(

2–3; bottom). We compare against the

simulations presented here (grey points) with explicit feedback, and obser-

vational constraints (lines). Even sub-grid models which are ‘successful’

near

∼

∗

∼

0 overpredict

∗

(

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 exceptions appear

to be the newest generation of sub-grid models which have been explicitly

adjusted to mimic the effects of radiative feedback as well as SNe, seen

in our explicit-feedback models: this includes Stinson et al. (

2013

), Aumer

et al. (

2013

), Ceverino et al. (

2014

) and Trujillo-Gomez et al. (

2013

previous work (Paper

II).

Whether this is unique or not remains to

be tested; the phase structure and other properties of outflows and

the CGM in such models can be very different from those predicted

here, even for the same mass-loading efficiencies (discussed fur-

ther below). It will be particularly interesting to see whether other

recently developed sub-grid models such as that in Agertz et al.

(2013), also incorporating the effects of radiative feedback but via

very different prescriptions, will also agree well with observations

There have also been interesting results from the re-tuned wind model

of Oppenheimer & Dav

2006

) used more recently in slightly different

forms in Torrey et al. (

2014

), Marinacci, Pakmor & Springel (

2014

), and

Hirschmann et al. (

2013

). However, in this model, the wind outflow rates

are set

explicitly

by-hand (and in fact the most recent scalings used were

adjusted based on comparison to the sub-grid models including radiative

feedback), and then tuned to reproduce the observed mass function. So this

is essentially what we attempt to

predict

here.

Figure 8.

KS-law, observed (Kennicutt

1998

; Bigiel et al.

2008

; Genzel

et al.

2010

; Daddi et al.

2010

, yellow shaded range) and simulated (points

as Fig.

). We emphasize that this is a prediction: the instantaneous star

formation efficiency per dynamical time in dense gas is 100 per cent in the

simulations, but the emergent KS-law, as a consequence of feedback, has

an efficiency a factor

∼

50 lower. As shown in Paper

and Hopkins et al.

(

2013c

), this is insensitive, with resolved feedback models, to the small-

scale star formation law, and entirely determined by stellar feedback. ‘No

feedback’ models lie a factor

∼

50 above the observations; ‘no radiation’

and ‘no SNe’ models (Fig.

) lie a factor

∼

10 above observations.

at both low and high redshifts. In any case, these comparisons – and

the results from this new generation of sub-grid models – highlight

that some accounting for non-SNe feedback is critical.

5.6 Instantaneous suppression of star formation (at fixed gas

densities)

We now examine galaxy SFRs. In the previous section, we showed

that the

integrated

star formation is suppressed with feedback. But

equally important is that feedback suppresses

instantaneous

SFRs

in galaxies. This is manifest in the KS relation, shown in Fig.

We plot the simulations at all redshifts (the redshift evolution is

insignificant), and compare to observations at a range of redshifts

(which also find little or no evolution).

The predicted KS-law agrees well with observations at all red-

shifts. As shown in Paper

I–Paper

III, this emerges naturally as a

consequence of feedback, and is

not

put in by hand. Recall that the

instantaneous star formation efficiency (star formation per dynam-

ical time) in dense gas in the simulations is 100 per cent; however,

the global star formation efficiency is

∼

2 per cent. This difference

arises because at

∼

2 per cent efficiency, feedback injects sufficient

momentum to offset dissipation (indeed, given the same feedback,

we obtain the identical KS-law independent of the details of our

small-scale star formation law; see Paper

I; Hopkins et al.

2013c).

We define

SFR

∗

/π R

SFR

(where

∗

is the total SFR and

SFR

the half-SFR radius) and

gas

/π R

SFR

(where

gas

is the gas mass

within the 90 per cent SFR radius). Defining both

∗

and

gas

within

SFR

or the stellar effective radius shifts the points along the relation.

We compile the observed local galaxies in Kennicutt (

1998

) and Bigiel

et al. (

2008

), and high-redshift galaxies in Genzel et al. (

2010

) and Daddi

et al. (

2010

); shaded region shows the 90 per cent inclusion range at each

gas

from the compilation. As discussed in those papers, there is no signif-

icant offset between the high and low-redshift systems at fixed

gas

MNRAS

445,

581–603 (2014)

at California Institute of Technology on December 1, 2014

http://mnras.oxfordjournals.org/

Downloaded from