An instability of feedback regulated star formation in galactic nuclei

5.2 MB·

by  Paul Torrey

Mon.Not.R.Astron.Soc.000,000–000(0000) Printed28January2016 (MNLATEXstylefilev2.2)
An instability of feedback regulated star formation in
galactic nuclei
Paul Torrey(cid:63)1,2, Philip F. Hopkins2, Claude-Andr´e Faucher-Gigu`ere3,
Mark Vogelsberger1, Eliot Quataert4, Duˇsan Kereˇs5, & Norman Murray6,7
1MIT Kavli Institute for Astrophysics & Space Research, Cambridge, MA, 02139, USA
6 2TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA
1 3Department of Physics and Astronomy and CIERA, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA
0 4Department of Astronomy and Theoretical Astrophysics Center, University of California Berkeley, Berkeley, CA 94720, USA
2 5Department of Physics, University of California at San Diego, 9500 Gilman Drive, La Jolla, CA 92093, USA
6Canadian Institute for Theoretical Astrophysics, 60 St. George Street, University of Toronto, ON M5S 3H8, Canada
n
7Canada Research Chair in Astrophysics
a
J
28January2016
6
2

ABSTRACT
]
A We examine the stability of feedback-regulated star formation (SF) in galactic nuclei
andcontrastittoSFinextendeddiscs.Ingalacticnucleithedynamicaltimebecomes
G
shorterthanthetimeoverwhichfeedbackfromyoungstarsevolves.Weargueanalyt-
.
h icallythatthebalancebetweenstellarfeedbackandgravityisunstableinthisregime.
p We study this using numerical simulations with pc-scale resolution and explicit stel-
- lar feedback taken from stellar evolution models. The nuclear gas mass, young stellar
o
mass, and SFR within the central ∼100 pc (the short-timescale regime) never reach
r
t steady-state, but instead go through dramatic, oscillatory cycles. Stars form until a
s
criticalsurfacedensityofyoungstarsispresent(suchthatfeedbackoverwhelmsgrav-
a
[ ity), at which point they begin to expel gas from the nucleus. Since the dynamical

times are shorter than the stellar evolution times, the stars do not die as the gas is
1
expelled, but continue to push, triggering a runaway quenching of star formation in
v
thenucleus.Howevertheexpelledgasislargelynotunboundfromthegalaxy,butgoes
6
into a galactic fountain which re-fills the nuclear region after the massive stars from
8
1 the previous burst cycle have died off (∼50 Myr timescale). On large scales (>1 kpc),
7 the galaxy-scale gas content and SFR is more stable. We examine the consequences
0 of this episodic nuclear star formation for the Kennicutt-Schmidt (KS) relation: while
. a tight KS relation exists on ∼1 kpc scales in good agreement with observations, the
1
scatter increases dramatically in smaller apertures centered on galactic nuclei.
0
6
Key words: methods: numerical – stars: formation – galaxies: evolution – galaxies:
1
formation – galaxies: high- redshift – galaxies: ISM – galaxies: starburst
:
v
i
X
r
a 1 INTRODUCTION lation.TheKSrelationconsidersthecurrentstarformation
rate surface density as a function of the available fuel for
Observations of star formation in our Galaxy and in extra-
star formation measured through the total (or molecular)
galactic systems show that stars are formed inefficiently. A
gassurfacedensity.Althoughthereisstilldebateaboutthe
commonapproachforparameterizingthestarformationrate
detailed slope and normalization of the KS relation, it is
(SFR) efficiency is to express ρ˙ =(cid:15)ρ /t where ρ˙ is the
∗ gas ff ∗ generallyagreedthat–whenmeasured/averagedoverlarge
starformationrateperunitvolume,ρ isthevolumeden-
gas areas within galaxies – only a few percent of gas is able to
sityofgas,t isthelocalgasfreefalltime,and(cid:15)isthestar
ff be converted into stars per free fall time (e.g., Zuckerman
formation rate efficiency (i.e. the fraction of gas that turns
& Evans 1974; Leroy et al. 2008; Bigiel et al. 2008). Un-
into stars per galactic free fall time). Since the volumetric
derstanding and explaining this inefficiency forms a central
gas and star formation rate densities are not easily observ-
componentofstarformationandgalaxyformationresearch.
able, this inefficiency is observationally probed through the
Kennicutt-Schmidt(KS;Schmidt1959;Kennicutt1998)re- Twomainclassesofmodelshavebeenpresentedtoex-
plainthelowefficiencyofstarformationandtheformofthe
KSlaw.Inoneclassitisarguedthatthesmall-scaleproper-
(cid:63) E-mail:[email protected] ties of supersonic ISM turbulence are able to disrupt cloud
(cid:13)c 0000RAS
2 Torrey et al.
collapse, giving rise to low SFR efficiencies (e.g., Krumholz
& McKee 2005; Federrath & Klessen 2012). These models 104 Stellar Winds
relate the low efficiency of star formation to the statistical Supernovae
properties of the density fluctuations in supersonic turbu- Photons
lence.However,sincetherate-limitingstepforgravitational
collapse is the formation of dense gas clouds in the first
place, this is an incomplete explanation of the galaxy-scale )fl103
KS law. M
Another approach has been to develop analytic global /c
“equilibrium” models that explain the low efficiency of star fl
L
formation and form of the KS law (Silk 1997; Thompson (

et al. 2005; Ostriker & Shetty 2011; Faucher-Gigu`ere et al. m∗102
2013). The common component of the equilibrium models /
˙P
is that the interstellar medium (ISM) remains in a quasi
stable state by balancing feedback from young stellar pop-
ulations with the global pressure of the self gravity of the
disk (Silk 1997; Hopkins et al. 2011a). While thermal pres- 101
sure or radiation pressure contribute to the disk support,
play a role in disrupting the giant molecular cloud (GMC)
complexeswherestarsareformed(e.g.,Murrayetal.2010), 105 106 107 108
and could be a source of turbulence, turbulent pressure is
Time (years)
thought to provide the dominant source of support against
the vertical collapse in gas-rich galaxies (Ostriker & Shetty
Figure 1. The direct momentum input taken from the STAR-
2011; Faucher-Gigu`ere et al. 2013), except possibly in the
BURST99stellarevolutionmodelsusingaKroupaIMF.Thefirst
most optically thick inner parts (Thompson et al. 2005). ∼3×106yearsarecharacterizedbynearlyconstantfeedbacklev-
Equilibrium star formation models are able to derive star elsdominatedbyradiationpressurefromyoungstars.Thedirect
formationefficienciesbyassertingthattheturbulentenergy momentum injections levels begin to evolve only after massive
injection rate from young stellar populations balances the starsbeginmovingoffthemainsequencewhichleadstoacontin-
turbulent energy dissipation rate. The derived star forma- ualdecreaseinthedirectradiationpressuremomentuminjection
tionefficienciesscalewithmeasurablepropertiesofgalactic and an increase in the supernova momentum injection. We note
that the momentum injection values shown are direct momen-
disks and are consistent with observed star formation effi-
tum injection values and do not account for any gains owing to
ciencies (Faucher-Gigu`ere et al. 2013).
trapping of photons (i.e. multiple scattering events), gains from
These models implicitly assume that (1) gas in the lo-
pressure confinement of hot gas (the “b” factor, as defined in
cal disk responds to feedback, (2) the strength of feedback
Equation 10), or losses owing to inefficient coupling. We note
scales with the star formation rate, and (3) the star forma- that the boost to the direct momentum injection is likely mod-
tionrate–andthereforefeedback–respondsinturntothe est for radiation pressure owing to radiation hydro instabilities,
gas properties. If feedback is “too strong,” gas will be ex- butcanprovideanorderofmagnitudeincreaseinthesupernova
pelled or pushed out to such low densities that it cannot contributiontothetotalmomentuminjection.
formnewstars,sothefeedback“supply”willcease.Iffeed-
back is “too weak,” gas collapses rapidly and forms more
mechanisms depends on the physical conditions of the gas.
stars, injecting more feedback until collapse is halted.
Radiation pressure (Murray et al. 2010) and photoioniza-
This feedback cycle can only hold if the gas adjusts
tion (Whitworth 1979; Krumholz et al. 2006; Walch et al.
efficiently to the presence of feedback, and the feedback en-
2012;Salesetal.2014)areoftenidentifiedasbeingrespon-
ergy/momentum injection rate efficiently adjuststochanges
sible for disrupting the dense star forming gas, but super-
inthegasdensityviathestarformationrate.Aslongasthe
novaeprobablydominatetheturbulentmomentuminjection
feedback is sufficiently strong and efficiently coupled to the
that balances the disk vertical collapse under most condi-
gas, the gas will adjust to the presence of feedback. How-
tions (Ostriker & Shetty 2011).
ever, there are physical regimes where feedback evolves on
Figure 1 shows the direct momentum injection budget
timescales that are (much) longer than the local dynamical
associatedwithasingleagestellarpopulationasafunction
time. In these regimes, the feedback energy/momentum in-
oftimesincebirthtakendirectlyfromSTARBURST99(Lei-
jection rate will not respond to changes in the gas density,
therer et al. 2010). Radiation pressure dominates the early
potentiallybreakingthefeedbackloopinvokedinthesestar
stellar feedback budget until the death of the most mas-
formation models.
sive stars (∼3 Myr) and drops to ∼10% of its initial value
For stellar feedback, the primary sources of direct mo-
by ∼10 Myr. Supernovae dominate the direct momentum
mentumandenergyinjectionintothelocalenvironmentare
injection budget starting at ∼3 Myr and continue to domi-
radiation pressure, stellar winds, photoionization, and su-
nate until the time after which the least massive stars that
pernovae(Murrayetal.2010;Hopkinsetal.2011b),andcos-
die in supernovae move off the main sequence (∼40 Myr)1.
micrays(e.g.,Jubelgasetal.2008;Uhligetal.2012;Booth
From Figure 1, the total momentum injection rate remains
etal.2013;Salemetal.2014).Magneticfieldsmayalsoact
to suppress star formation regulation mechanisms (Dolag
et al. 1999; Wang & Abel 2009; Pakmor & Springel 2013; 1 For a Kroupa IMF with M = 8M(cid:12) assumed to be the lower
Marinacci et al. 2015). The relative importance of these masslimitforsupernovaprogenitors.
(cid:13)c 0000RAS,MNRAS000,000–000
Unstable feedback in galactic nuclei 3
nearly constant until ∼3 Myr, and remains variable but wherehisthediskscaleheightandρ isthegasdensity.We
g
order-of-magnitude constant until ∼40 Myr. The constant then assume that supersonic turbulent pressure is the pri-
early momentum injection rates allow us to define an in- mary source of pressure support in the disk with turbulent
teresting regime: regions of galaxies that possess dynamical gas velocity dispersion c , in which case the rate of (volu-
t
times less than the timescale over which stellar feedback metric) turbulent energy dissipation is
levels adjust. In these regimes, feedback can not efficiently
u˙ ≈ρ c2Ω. (2)
adjusttochangesinthegasdistributionand/orthecurrent g t
star formation rate. In this paper we focus on the nuclear Under the assumption that the system approaches an equi-
regionsofgalaxieswherethesmallspatialscalesandcentral librium,thisenergydissipationrateneedstobebalancedby
black hole ensure consistently short dynamical times. feedback from young stars. There are several possible ways
The understanding of star formation in galactic nuclei to write down the feedback injection rate. Here, we express
thatwepresentinthispaperwillsetthestageforfollowup the energy injection rate from young stars as
studiesfocusedontheaccretion,growth,andfeedbackasso-
p c f P /m Σ˙ c
ciatedwithsupermassiveblackholes.Previousstudieshave u˙ ≈ ∗ t = p ∗ ∗ ∗ t. (3)
considered the impact of stellar and black hole feedback in h 4h
concert on the central ∼100 pc nuclear disk (e.g., Hopkins The second equality follows from expressing the turbulent
et al. 2015). However, those studies lacked the large-scale pressure produced by stellar feedback as
galactic disk that is important for understanding the long
P
termbehaviorofsupermassiveblackholeaccretionandfeed- p∗ =fp4m∗ Σ˙∗ (4)
back. ∗
Thestructureofthispaperisasfollows:InSection2we whereP∗/m∗ istheradialmomentuminjectedperunitstel-
reviewthebasiccomponentsofequilibriumfeedbackmodels. lar mass formed (Ostriker & Shetty 2011; Faucher-Gigu`ere
This allows us to concretely explore the assumptions that et al. 2013) and fp parameterizes the efficiency with which
break down for the equilibrium model at short dynamical theradialmomentumdrivesturbulentenergyinjection,but
times. We propose alterations to the equilibrium feedback isoflikelyoforderunity(Ostriker&Shetty2011).Adopting
models that hold within short dynamical time regions and the above expression of the young star pressure/turbulent
conclude that star formation within short dynamical time energy injection rate does not require any assumptions
regions is likely to be inherently bursty. In Section 3 we about the origin of the momentum input, and can accom-
describethenumericalsimulationsthatweusetoexplorethe modate multiple momentum input sources (e.g., radiation
propertiesofnuclearstarformation.InSection4wepresent pressure, supernovae, stellar winds, etc.). For SNe, typical
theresultsofournumericalsimulations,withafocusonthe values for P∗/m∗ are of order ∼3000 km s−1 (see Ostriker
burstynatureofnuclearstarformationandimplicationsfor & Shetty 2011, for additional details).
observations of the KS relation. In Section 5 we present a Enforcing pressure balance between the self gravity of
discussion of our results with implications for observations the disk and turbulent pressure (i.e. between Equations 1
of nuclear stellar disks as well as the observed KS relation. and 4) we find the feedback self-regulated equilibrium con-
In Section 6 we summarize and conclude. dition
f P /m Σ˙
ρ c2 ≈ p ∗ ∗ ∗. (5)
g t 4
Using c = hΩ and Σ = 2hρ , this can be rearranged in
t g g
2 ANALYTIC ARGUMENTS
termsoftheequilibriumstarformationratesurfacedensity
2.1 The Equilibrium Model 2c
Σ˙ ≈ηΣ Ω=Σ Ω t (6)
∗ g g f P /m
Analytic models have been developed to explain the effi- p ∗ ∗
ciency of star formation and form of the KS law. In global whereη=2c /f (P /m )isanefficiencythatdescribesthe
t p ∗ ∗
equilibrium models the star formation efficiency is derived fraction of gas that can convert into stars in a dynamical
byassertingthatturbulentenergyinjectionfromyoungstel- time.Thisexpressionshowsthatthegalacticstarformation
lar populations balances the turbulent energy dissipation efficiency is directly related to the gas turbulent velocity
rate.OneofthefirstmodelstodothiswasThompsonetal. dispersion in the disk.
(2005),wheretheyexaminedradiationandmechanicalfeed-
backfromyoungstellarpopulationswithafocusonthecen-
tral100pcofastarburstdisk.Weadoptasimilarsetuphere, 2.2 Where the equilibrium models break down
which includes an isothermal spherical potential – which
The equilibrium model as described above breaks down in
could represent either a halo or bulge component. The pro-
two regimes.
jected surface density follows a power law Σ = σ2/πGr
tot
whereσisthevelocitydispersionofthepotential.Themass
profile is given√by Mtot = 2σ2r/G, with an angular fre- 2.2.1 Equilibrium Breakdown from Short Dynamical
quency of Ω = 2σ/r set to enforce radial centrifugal bal- Times
ance.Thespecifiedgasfractionisf =M /M =Σ /Σ .
g g tot g tot
The equilibrium model breaks down when feedback is un-
For such a model the disk mid plane pressure required to
abletomodulateitsstrengthonsufficientlyshorttimescales
balance gravity can be approximated by
(i.e.thegasdynamicaltimescale).Feedbacklevelsaresetby
p≈ρ h2Ω2. (1) youngstellarpopulations,whichevolvebasedontherateat
g
(cid:13)c 0000RAS,MNRAS000,000–000
4 Torrey et al.
1011 104
Total Mass Total Total
Disk (gas+stars) Disk (gas+stars) Disk (gas+stars)
1010 BH BH 103 BH
Bulge Bulge Bulge
Halo Halo Halo
109 102
()MMencfl108 ()Vkm/srot102 ()tMyrdyn 101
107 100
106 10-1
105 101 10-2
10-3 10-2 10-1 100 101 10-3 10-2 10-1 100 101 10-3 10-2 10-1 100 101
r (kpc) r (kpc) r (kpc)
Figure 2. The mass enclosed (left), rotational velocity (center), and dynamical time (right) are shown as a function of radius for the
D1(lowmass)initialdisk.Theblackhole’sdynamicalrelevanceatscalesr(cid:46)10-30parsecscanbeidentifiedfromanyoftheplots.The
horizontal lines in the dynamical time plot indicate 3 and 40 Myr. Dynamical times are calculated using the enclosed mass profiles to
determinethecircularvelocities.
which massive stars evolve. Here, we are specifically con- To explore this regime, we define the convenient quan-
cernedwiththetimescaleoverwhichradiationpressureand tity a (t) = P˙/m , which has units of acceleration. This
∗ ∗
energy/momentum injection from supernovae decline as an quantityisplottedinFigure1asafunctionoftime.Wecan
initially young stellar population ages. Radiation pressure express the turbulent pressure provided by young stars as
begins to drop when a stellar population is just ∼ 3 Myr
p =Σ a (t)·bf /4, (9)
oldwhilesupernovainjectsignificantmomentumuntil∼40 ∗,y ∗,y ∗ p
Myr (Figure 1). Both of these feedback mechanisms have where Σ is the surface density of “young stars” which
∗,y
fallen to a small fraction (∼ 0.01) of their initial value by dominatethestellarfeedbackbudget.Throughout,weadopt
the time a stellar population is ∼ 100 Myr old. In regimes adefinitionofyoungstarsasbeingthosestarswithagesless
where the dynamical time is much shorter than the stellar than 40 Myrs. The quantity b accounts for the increase in
evolution time (i.e. τSE ∼ 106−107 years; see Figure 1)2, the momentum provided by stars via trapping of IR pho-
feedbackwillnotreactsufficientlyquicklytochangesinthe tons (Thompson et al. 2005) or by SN acquired during the
currentgasdensityorstarformationrate,butratherremain Sedov-Taylor phase (e.g., Cioffi et al. 1988; Martizzi et al.
strong for a stellar evolution timescale. 2015; Gatto et al. 2015; Kim & Ostriker 2015). The total
Fortheonezonemodelusedinthissection,thedynam- momentum injection per unit stellar mass (P /m , used in
∗ ∗
ical time is the previous subsection) is related to a (t) via
∗
√
2πr (cid:18) r (cid:19)(cid:18) σ (cid:19)−1 (cid:90) ∞
τ = ∼106yr . (7) P /m = a (t)·bdt (10)
dyn σ 100pc 200km/s ∗ ∗ ∗
0
Comparing this with a stellar evolution timescale τ we where the boost factor b is folded into the total momentum
SE
find that the dynamical time will be much shorter than the injection. Setting the turbulent pressure from young stars
stellar evolution time (τ (cid:28)τ ) when equal to the mid-plane pressure given by equation (1), we
dyn SE
find the critical surface density of stars
(cid:18) (cid:19)(cid:18) (cid:19)
τ σ
r(cid:28)1kpc 10MSEyr 200km/s . (8) Σcrit ≈ 2ΣgctΩ = πGQΣ2g . (11)
∗,y f a (t)·b f a (t)·b
p ∗ p ∗
When this condition is met, massive stars contributing sig-
This result is similar to the argument in Fall et al. (2010)
nificantlytoSNeandradiationpressurelivelongerthanthe
and Murray et al. (2010) for when individual star clusters
gasdynamicaltimeandthusfeedbackfromyoungstarswill
disperse giant molecular clouds, but now we apply it to a
persist for several dynamical times. For reasonable τ val-
SE
galaxy nucleus which is allowed to have many independent
ues(e.g.∼10Myr)thedynamicaltimecanbemuchshorter
self-gravitating clouds.
than τ out to a significant fraction of the central 1 kpc.
SE
We briefly note that the same critical condition can be
Ifthefeedbackfromyoungstarsexceedsgravitywithinthis
derived by considering the energy – rather than pressure –
region,thenthesystembecomesunstable/unboundandre-
balance. Specifically, expressing the turbulent energy injec-
mains in that state for several dynamical times until the
tion rate as
young stars responsible for feedback move off the main se-
quence. p c f a (t)·bc Σ
u˙ ∼ ∗,y t = p ∗ t ∗,y. (12)
h 4h
andequatingthiswiththeturbulentenergydissipationrate
2 Inwhatfollows,ouranalysisdoesnotrequirethatweidentify given in eqn. 2 yields the same critical surface density of
aspecificvalueofτSE. young stars. This equivalence is built into the model with
(cid:13)c 0000RAS,MNRAS000,000–000
Unstable feedback in galactic nuclei 5
the assumption that the disk mid-plane pressure is related 2.2.2 Equilibrium Breakdown from Insufficient Feedback
directly to the turbulent energy injection and dissipation
Expressing the critical young star surface density as a frac-
rates.
tion of the gas mass yields
Defining the young stellar fraction as
Σcrit 2c Ω Qσf
Σ ∗,y = t =Ω√ g . (16)
f∗,y ≡ Σ∗,y (13) Σg fpa∗(t)·b 2fpa∗(t)·b
tot
An interesting limit is found by considering when the frac-
wecanwritedownthecriticalyoungstellarfractionorcrit- tionofgasthatmustbeturnedintostarsinordertoachieve
ical young star surface density in terms of the total surface feedback balance is of order unity (i.e. Σcrit/Σ ∼ 1). This
∗,y g
density limit occurs for
f a (t)·b M
fcrit ≡ Σc∗r,yit ∼ 2fgctΩ ∼ πGQfgΣg. (14) Σg = pπ∗QG ≈2.7×1011kpc(cid:12)2 (17)
∗,y Σ f a (t)·b f a (t)·b
tot p ∗ p ∗
for b = 20 and a (t) = 1000L /cM . At this gas surface
∗ (cid:12) (cid:12)
If we assume a boost factor of b = 20 and a (t) = density, feedback can efficiently oppose gas collapse only
∗
1000L /cM this reduces to by converting an order unity fraction of the gas mass into
(cid:12) (cid:12)
stars. This critical surface density value approximately cor-
(cid:18) Σ (cid:19) responds to the maximum observed stellar surface density
fcrit ≈1.5×10−2f Q g . (15)
∗,y g 103M pc−2 in a wide range of dense stellar systems (Hopkins et al.
(cid:12)
2010) and represents a feedback failure mode that will be
Ifthiscriticalsurfacedensityofyoungstarsisexactlymain- explored in more detail in a forthcoming paper (Grudic et
tained,feedbackdriventurbulenceanddiskself-gravitywill al., in prep).
beinequilibrium.However,itfollowsthatovershootingthis
critical young stellar surface density will lead to outflows,
which will lower the system’s gas fraction, and lead to a 3 NUMERICAL METHODS
runaway outflow condition that will persist for several local
dynamicaltimes.3 Forthisreason,wearguethatthecritical We investigate the galactic nuclear burst/quench cycles de-
scribed above using numerical simulations.
surface density of young stars is an unstable critical point
which will likely lead to a full shut down of star formation
with the system going into an outflow state for a stellar 3.1 Simulation Code
evolution feedback timescale.
Aftertheyoungstellarpopulationshaveagedandtheir Thesimulationspresentedinthispaperwereperformedus-
feedback levels have been reduced, gas will return to the ingtheN-bodyhydrodynamicscodeGIZMO(Hopkins2015).
nuclearregion.Gravitationaltorquingofgasviaasymmetric GIZMO is originally derived from GADGET (Springel 2005)
disk potential features can drive gas mass influx rates of and contains modifications centered around the ability to
order M˙ ∼ −|a|βΣ R2Ω where |a| is the asymmetric disk employseveralfundamentallydifferenthydrosolvers.Inthis
g
mode strength and β =1−2 depending on whether we are paper, we employ GIZMO using the meshless-finite-mass
in the linear or non-linear torquing regime (see Hopkins & method to solve the hydrodynamic equations of motion.
Quataert 2011, for more details). This can repopulate the In addition to gravity and hydrodynamics, all simula-
central gas reservoir within a few dynamical times, which tionsincludealistofphysicalprocessesimportanttogalaxy
will allow for a repetition of the central gas blowout cycle. formation.Thespecificgalaxyformationmodelusedinthis
This result stands in contrast to the “self-regulating” paper is the FIRE feedback model (Hopkins et al. 2014).
feedback balance that applies to most of the disk. The key This model has been explored extensively in previous pa-
difference is whether or not feedback driven turbulent pres- pers simulating the evolution of galaxies within cosmolog-
sureadjustsrapidlyorslowlycomparedtothelocaldynami- ical volumes. The FIRE model is capable of reproducing
caltime.Typicalgalacticdynamicaltimesat∼kpcdistances manyobservedgalaxypropertiesincludingthestellarmass-
are of order ∼100Myrs (i.e. longer than the timescale over halo mass relation, Kennicutt-Schmidt law (Hopkins et al.
whichstellarfeedbackevolves).Upwardperturbationstothe 2014), the covering fractions of neutral hydrogen in the ha-
SFRdonotleadtoblow-outundernormal‘equilibrium’con- losofz=2-3LymanBreakGalaxies(LBGs,Faucher-Gigu`ere
ditions because stars die on timescales much shorter than etal.2015),andcanself-consistentlygenerategalacticwinds
the local dynamical time. Where this is the case, there is thatregulatedgalaxymassgrowthconsistentwithobserva-
not enough time for minor SFR perturbations to affect the tional requirements (Muratov et al. 2015).
galaxy dynamics. Details of the FIRE feedback model are given in Hop-
kinsetal.(2014),whichwesummarizebrieflyhere.Allsimu-
lationsincluderadiativegascoolingandstarformationwith
associated feedback. Supermassive black hole particles are
included in the present simulations because their gravity is
dynamically important. However, no black hole feedback is
3 We caution that outflows can occur throughout our simulated
included. Gas cools radiatively under the assumption of lo-
disksevenwhenthisconditionisnotmet.Theconditionherewill
resultinastrong,nearlycomplete,depletionofthelocalcentral cal thermodynamic equilibrium down to 10 K. Dense/cold
gasreservoirwhereasasteadylevelofgasoutflowcanbepresent gas clouds are allowed to form stars if they are locally self-
elsewhereinthediskevenwhenthisconditionisnotmet. gravitating.Thestarformationrateisgivenbyρ˙∗ =ρmol/tff
(cid:13)c 0000RAS,MNRAS000,000–000
6 Torrey et al.
wherethemolecularfraction,f ,isinferredasafunctionof
H2
localgascolumndensityandmetallicityfollowingKrumholz
& Gnedin (2011). Newly formed stars provide feedback to
the ISM through thermal heating via supernovae, photo-
ionization,localandlongrangeradiationpressure,andstel-
larwinds.Eachofthesefeedbacksourceshasawelldefined
fiducial level that is adopted directly from STARBURST99
givenastellarparticle’sageandmetallicity(Leithereretal.
2010).
3.2 Initial Conditions
We employ initial conditions that are modeled after red-
shift z = 0 isolated star forming late type galaxies. The
initialdistributionofparticlesissampledbasedonthescal-
ing relations of Mo et al. (1998) via the procedure outlined Figure 3. The projected gas surface density is shown at three
in Springel & White (1999). We construct three isolated snapshots during the evolution of the D1 system. The top row
galaxiesthataresimilarintheirphysicalcharacteristicsbut shows synthetic SDSS-g, -r, -i band images made using a simple
varied in initial mass to use as initial conditions. line of sight attenuation model. The bottom row shows the gas
surface density, with magenta indicating atomic/molecular gas
ThethreegalaxieshavetotalmassesofM ={1.39×
tot (10K<T <8×103K),greenindicatingwarmionizedgas(8×103
1011,3.37×1011,1.39×1012}M .Eachofthesegalaxiescon-
(cid:12) K<T <8×105 K),andyellow(thoughlittleispresent)hotgas
sistsofadarkmatterhalo,stellardisk,stellarbulge,gaseous
(T >105K).Thetimeofeachsnapshotismarkedinthetoprow
disk(with20%gasfraction),andcentralsupermassiveblack
and the scale bar of 1 kpc is also indicated. From the initially
hole.Thethreesystemsaresimilarinthattheyallcontain4 smooth gas distribution a turbulent ISM develops, with dense
percentoftheirmassinthedisk,aconstant(small)fraction starformingregionsbeingcontinuallyformedbyself-gravityand
in the central supermassive black hole, 1 percent of their destroyedthroughfeedbackfromyoungstars.
mass in the bulge, and the rest of the mass in the halo (see
Table1).Thestellarandgaseousdisksaregivenexponential
surface density profiles with the same disk scale length set Table 1. The properties of the initial conditions used for the
bythehalospinparameterofλ=0.04(Rd ={1.7,2.3,3.7} three isolated disks presented in this paper. Below, Mtot is the
kpc). The relative masses in stars and gas is set by the ini- total mass of each system and mp is the baryon particle mass.
tial disk gas fraction, for which we adopt 20 percent. The Adaptivegravitationalsofteninglengthsareemployedwith(cid:15)be-
disks are initialized to be nearly in equilibrium (e.g., with ing the minimum (Plummer equivalent) force softening length
usedforeachparticletype.
Q ∼ 1). The dark matter halo and stellar bulge are setup
using a Hernquist (1990) profile initial density distribution. ID GasDisk Halo StellarDisk
AnNFW(Navarroetal.1997)equivalentconcentrationpa-
rameter of C = 10 is used for all systems. The bulge scale
lengthissettobe0.2timesthestellarandgaseousdiskscale D1 Mtot[M(cid:12)] 1.11×109 1.32×1011 4.45×109
lengths. The single supermassive black hole is included in mp[M(cid:12)] 103 8.0×104 3.0×103
(cid:15)[pc] 1 10 2
eachgalaxyasacollisionlessparticlethatisinitiallyplaced
atthecenterofeachgalaxy.Sincethemassofthesupermas- D2 Mtot[M(cid:12)] 2.72×109 3.23×1011 1.09×1010
siveblackholeismuchlargerthanthetypicalgasorstellar mp[M(cid:12)] 103 8.0×104 3.0×103
particle mass, it remains close to the potential minimum of (cid:15)[pc] 1 10 2
thehostgalaxyviaself-consistentlycaptureddynamicalfric- D3 Mtot[M(cid:12)] 1.11×1010 1.32×1012 4.45×1010
tion without any additional special treatment. Additional mp[M(cid:12)] 103 8.0×104 3.0×103
valuessuchastheforcesoftening,massresolutionandnam- (cid:15)[pc] 1 10 2
ingconventionthatweusedforthisinitialconditionsetare
summarized in Table 1.
Figure 2 shows the mass profiles, rotational velocity
curves,andradiallydependentdynamicaltimeforourfidu- Figure2.Weidentify3and40Myrlinesbecausethisroughly
ciallowmassgalaxy(D1).Thedynamicaltimesaredefined corresponds to the time when stars begin moving off the
usingtheenclosedmassprofilestocalculatethecircularve- main sequence and when the most of the energy from the
locity. Although we neglect black hole accretion and feed- newlyformedstellarpopulationhasbeendeposited,respec-
backinthesesetsofsimulations,thepresenceofthismassive tively. For this setup the dynamical time drops below ∼ 3
particleisimportantowingtoitsimpactonthedynamicsof Myrataradiusof∼100pc,whichchangesonlyslightlyfor
thenuclearregion.Thepresenceofthecentralsupermassive theotherinitialdisks.Withinthisregion,weexpectstarfor-
black hole dominates the mass enclosed and rotation curve mationandgassupplytobeunabletoachieveasteadystate
for the inner 10-30 parsecs. The rotational velocity curves solution owing to the short dynamical times. The presence
have a non-monotonic profile, but the dynamical time con- of a 107M black hole alone keeps the dynamical time be-
(cid:12)
tinually drops toward the galaxy center. The characteristic low3Myrforthecentral∼100pc,butamassive/compact
stellar evolution timescales of ∼ 3 Myr and ∼ 40 Myr are stellar bulge could play a similar role in shaping the inner
identified with horizontal dashed lines in the right panel of galactic potential.
(cid:13)c 0000RAS,MNRAS000,000–000
Unstable feedback in galactic nuclei 7
4 RESULTS However, once the feedback from this young stellar popu-
lation recedes (two rightmost panels) dense/cold molecular
4.1 Gas Distribution and Galaxy Structure
gas quickly returns to the central region.
Figure 3 shows a time sequence of the D1 galaxy evolved From Figure 5 we identify that star formation occurs
from 1 to 2 Gyr. The top panel shows a synthetic opti- in several distinct self-bound dense gas clouds within the
cal image of the galaxy while the bottom panel shows the central ∼100 pc. However, the young stars rapidly decou-
gas content. As indicated in both panels, throughout the ple from their birth clouds owing both to feedback and the
entire simulation the disk remains globally stable without shortdynamicaltime.Therapidmixingoftheyoungstellar
any strong bars developing (although there is a prominent populations drives the whole central region to react homo-
spiral pattern). This global stability is by construction. As geneously to the presence of strong feedback. Despite the
the simulation proceeds gas cools and collects into dense presence of multiple self-bound star forming clumps, the
star forming clumps which are embedded in a volume fill- feedback response of the gas in the whole nuclear region
ing warm phase of the ISM. The color of the gas distribu- acts like that of a single molecular cloud: stars form until a
tion on the bottom panel highlights this multi-phase struc- critical surface density of young stars is present (such that
ture:magentagasindicatesatomic/molecularmaterialwhile feedback overwhelms gravity) at which point they begin to
greenindicateswarmionizedgas.Thesestarformingclumps expel gas from the nuclei.
areeventuallydisruptedthroughfeedbackfromyoungstars.
The self-regulation of a multi-phase ISM through feedback
4.2 Star Formation Rates and Gas Masses Over
from stars using the FIRE feedback model has been dis-
Time
cussed extensively in Hopkins et al. (2014). Although we
onlyshowresultsherefortheD1disk,wenotethatD2,and Figure6showstheSFRwithinthecentral100pcfortheD1
D3 show qualitatively indistinguishable behavior. galaxy. We show the SFR when taken as an instantaneous
The multi-phase ISM stretches across the full radial valuefromthesimulation,aswellasaveragedoverthepast
range of each galaxy, including the central region. To high- 10 and 100 Myr by calculating SFR=M (t<t )/t where
∗ 0 0
lighttheevolutionoftheISMphasestructureinthenuclear t is the age of the star and t is the averaging timescale.
0
region,Figure4showsthegassurfacedensitywithinanar- We present smoothed definitions of the SFR for two rea-
row field of view focused on the galactic nucleus at several sons. First, the instantaneous SFR values obtained from
times. We can identify here the ISM structure extending the simulations rapidly vary and with significant magni-
down to (cid:46)100 pc scales. The gas density in the central re- tude. Smoothing the SFRs reduces the very high frequency
gioncanbeseentofluctuatethroughtheimages.Themost noiseallowingformoreclearidentificationofthemeanevo-
severe gas fluctuations are short lived and so we have se- lutionary trends. Second, observational SFR indicators are
lected times that highlight a period of central gas blowout. notsensitivetotheinstantaneousSFR,butratherthetotal
The typical blowout cycle begins with a concentration of mass of young stars formed over some recent time window.
molecular (purple) star forming gas in the central ∼ 100 The smoothed definition of the SFR is therefore more in
pc, as is found in the left most panel. This dense gas con- line with what is measured observationally through UV or
centrationformsstarsrapidlyuntilthefeedbackfromthese nebularemissionlineSFRtracers(e.g.,seeTable1ofKen-
young stars reaches the critical fcrit value. This transition nicutt & Evans 2012). In the rest of this paper, we adopt
∗,y
occurs very rapidly as the dynamical times on these scales a definition for the star formation rate as the average star
is only ∼1 Myr (see Figure 2). After gas is removed from formationrateoverthepast10Myr–whichwouldbemost
the central region, the gas remains at very low densities for similar to what would be measured though H-alpha SFR
∼20-30Myr(i.e.thethirdandfourthpanels).By∼50Myr measurements.
(therightmostpanel)moleculargasbeginstobuildupagain The episodic central star formation activity is directly
justoutsideof the centralregion.Asymmetricdisk features examinedinFigure7forthesimulatedsystems.Thepanels
are clearly present which foster the torquing of gas to the of Figure 7 show the time evolution of the gas mass, young
central 100 pc. stellar mass (i.e. M (t<40Myr)), star formation rate, and
∗
Inadditiontoconsideringthegasdistributiononthese criticalfeedbackvaluefromtoptobottom.Figure7presents
scales, we consider the distribution of young star clusters these quantities within the central 100 pc, while Figure 8
andtheirevolutionwithtime.Figure5showstheevolution shows the same quantities within the central 1 kpc.
ofthestellarpopulationsassociatedwiththissameblowout Oververylongperiodsoftime(i.e.∼1Gyrtimescales)
event.Specifically,thestellarpopulationthatwaslessthan we find the central region’s gas mass and star formation
10Myroldisidentifiedintheleftmostpanel,andthosesame rate is characterized by a relatively steady state. The time
particlesaretrackedforwardintime.Intheleftmostpanel, averagedstarformationrateorcentralgasmassevolvesonly
theyoungstars(heret<10Myr)takeonaclustereddistri- marginallyoverthefulltimeperiodshowninFigure7.This
bution that is (mostly) coincident with dense star forming isimpressive,giventhatinthenuclearregionsofthegalaxies
gas. In some cases, young stars are found in low density re- explored here are being evolved for hundreds of dynamical
gionswhentheyhavealreadybeenabletoblowaparttheir times.ThesteadynatureoftheSFRover∼1Gyrtimescales
birth cloud. Within 10 Myr (second panel from left), the is consistent with star formation being effectively regulated
stars have already randomized their locations significantly. by stellar feedback. Presumably, however, this would break
There is little correlation between these young stars and down on time scales much greater than ∼1 Gyr, where the
the dense gas distribution. By 20 Myr (central panel), the absence of IGM accretion would lead to disk gas depletion.
stars continue to deposit significant momentum and energy Onshortertimescalesthecentral100pcregionischar-
via supernovae but the gas is almost completely evacuated. acterized by an oscillatory gas mass, star formation rate,
(cid:13)c 0000RAS,MNRAS000,000–000
8 Torrey et al.
Figure 4. The projected gas surface density is shown at five snapshots during the evolution of the D1 system, focused around the
galacticnucleusduringaperiodofnucleargasblowout.Thewhitecircleindicatesacylindricalapertureof100pccenteredonthecentral
supermassiveblackhole,whichstayscloselypinnedtothepotentialminimum.Significantmoleculargasisinitiallypresentinthecentral
100 pc aperture. The concentration of cold/molecular gas (indicated with purple) in the nuclear region in the first (leftmost) snapshot
has an associated elevated level of nuclear star formation. Soon after the feedback from young stars hits the critical level, leading to a
near complete blowout of gas within the central 100 pc. The gas density remains suppressed for the next ∼30 Myr, while the feedback
fromyoungstarscontinuestobeefficientinspiteoftheabsenceofmuchongoingstarformation.
Figure 5. Same as Figure 4, but zoomed further on the galaxy nucleus and with the locations of young star particles identified with
black circles. The left panel indicates the locations of the youngest stars which are highly clustered and generally coincide with the
densegasinwhichtheywereformed.Subsequentpanelsshowthesamestarparticlesastheyage.Byafewtensofmillionsofyearsthe
distribution of young stars becomes more randomized in this region with the stellar populations being well separated from their birth
clouds.Thisispartiallydrivenbyfeedback,butaidedsignificantlybytheshortdynamicaltimesonthesescales.
and young stellar mass within the nuclear region. The be- mass of gas and young stars in the central regions shows
havior of each of these components is offset in time but no clear oscillatory behavior. The instantaneous star for-
closely related. The central gas mass increases in the ab- mationrate(notshown)showssignificantvariability,which
senceof currentstrong stellarfeedback untilstar formation becomes nearly smooth when averaged over 10 Myr peri-
becomesefficientinthecentralregionowingtothepresence ods (the plotted curves). The increased stability on 1 kpc
of sufficiently dense, self-bound gas clumps. Star formation scales is consistent with the analytic arguments outlined in
thenproceedsefficiently untilstellar feedbackis sufficiently Section2.Althoughfeedbackoperatestodisruptgasclouds
strong as defined by equation 11 to balance gas collapse. locallyfeedbackfailedtocollectivelyclearmaterialfromthe
Oncefeedbackissufficientlystrong,thegassupplyisrapidly region. The rate of gas consumption via star formation and
expelled (characterized by very sharp drops in the central gas expulsion via feedback are well balanced with the rate
gas content for any of the disks). The young stars that are ofgasreturnfromagingstellarpopulationsandinfluxfrom
responsibleforthestrongfeedbackcontinuetoinjectstrong larger radii producing stable gas mass evolution.
feedbackintothenuclearregionfor∼30-40Myrcausingthe
gasrecoveryinthisregiontoberelativelyslow.Theepisodic
sharp drops in the central gas content are correlated with 4.3 The Threshold for Burst/Quench Cycles
the periods of peak young stellar mass. The timescale for
In Section 2 we derived the surface density of young stars
this behavior is of order ∼100 Myr – though the detailed
thatwererequiredtodrivetheoutflows.Thesurfacedensity
duty cycle for this process is not very regular, and is differ-
ofyoungstarsnormalizedbythecriticalsurfacedensitythat
ent for our three disk models. The time evolution shown in
is expressed in Equation 11 is plotted as a function of time
Figure 7 can be compared back against the central gas sur-
inthebottompanelsofFigures7and8forourthreedisksas
facedensityplotsshowninFigure4.Bothfigurespresenta
evaluatedwithinthecentral100pcand1kpc,respectively.
complementarypictureofthecyclicalgasblowoutfromthe
Periods of time with fcrit <1 indicate that the system has
central region. ∗,y
low levels of feedback compared to what is required for gas
Thebehaviorofthesamequantitiesaveragedonlarger blowout and is therefore stable/collapsing while periods of
spatial scales becomes significantly more stable. Examining time with fcrit > 1 indicate that the system is feedback
∗,y
Figure 8 we find that averaged over the central 1 kpc, the dominated and is likely driving outflows. The central por-
(cid:13)c 0000RAS,MNRAS000,000–000
Unstable feedback in galactic nuclei 9
<SFR>
Inst
<SFR>
10Myrs
yr)10-1 <SFR>100Myrs
/
fl
M
(
R10-2
F
S
10-3
0.0 0.2 0.4 0.6 0.8 1.0
time (Gyr)
Figure 6.TheSFRwithinthecentral100pcasafunctionoftimefortheD1diskusingtheinstantaneousSFR(greenline),themass
of stars formed over the past 10 Myr (blue line), and the mass of stars formed over the past 100 Myr (black line). The mass of young
stars gives a smoothed estimate of the SFR and provides a better match to observational SFR probes. We exclusively use the mass of
youngstarsformedinthepast10Myr(blueline)todeterminetheSFRthroughouttherestofthispaper.
r<100 pc
D1
) 108
fl D2
M
( 107 D3

s
a
g
M 106
107
)
fl
M
106
(

y
,∗
M 105
) 100
r
y
/ 10-1
fl
M
( 10-2

R
SF 10-3
critf,y∗100
/
y
,∗
f
10-1
0.0 0.2 0.4 0.6 0.8 1.0
time (Gyr)
Figure7.Toptobottom:Thegasmass,youngstellarmass,starformationrates,andcriticalyoungstellarfractionasafunctionoftime
forthecentral100pc.Theepisodicstarformationinthecentral100pcisclear.Gascontentfirstbuildsupinthecentralregion,then
starformationrapidlybuildsapopulationofyoungstars,whichdrivesufficientlystrongfeedbacktodrivedownthecentralgasdensity.
tionofeachgalaxyspendsthemajorityofitstimeinasta- unity, there is a strong correspondence to a rapid drop in
ble/collapsingstate,bookendedwithbriefperiodsofstrong the central gas content of the system. Contrasting the pan-
feedback. We can compare the periods of time with high els within Figure 7 reveals that the fcrit values rise much
∗,y
fcrit values with the time evolution of the quantities pre- moresharplythantheyoungstellarmass.Thedrivingforce
∗,y
sented in the bottom panel of Figure 7. When fcrit exceeds behind this feature is that while the mass of young stars is
∗,y
(cid:13)c 0000RAS,MNRAS000,000–000
10 Torrey et al.
r<1 kpc
D1
) 108
fl D2
M
( 107 D3

s
a
g
M 106
107
)
fl
M
106
(

y
,∗
M 105
) 100
r
y
/ 10-1
fl
M
( 10-2

R
SF 10-3
critf,y∗100
/
y
,∗
f
10-1
0.0 0.2 0.4 0.6 0.8 1.0
time (Gyr)
Figure8.SameasFigure7butforthecentral1kpc.Theoscillatorypropertiesthatwerepresentforthecentral100pcarenotevident
forthecentral1kpc.Allmeasuredquantitiesarestableoverthislargeraperture,incontrasttothe100pccentralregion.
initially responsible for pushing fcrit above unity, once this disk dynamics, and less subject to influence from rapidly
∗,y
is achieved the rapid outflow of gas mass from the central outflowing material. The Q value hovers around unity con-
region drives up the derived turbulent velocity and drives sistentlywithtimeatlargeradii,indicatingsteady-stateself-
down the orbital frequency and gas mass. The combination regulation. The smaller radii bins have episodic and strong
ofthesethreeeffectsresultsintheverysharpriseofthefcrit increasesintheToomreQvalue,indicativeoflocaloutflows.
∗,y
values. Very large fcrit values are not long lived because (i) Theinnerregionsofthesimulateddisksexperiencebrief
∗,y
the gas fraction drops, (ii) new star formation ceases, and periodsoflargegasvelocitydispersionscoincidentwiththe
soonafter(iii)themassofyoungstellarmassdecreases(ow- presence of locally-driven outflows. Between these periods,
ing to aging). thegasvelocitydispersionsareoforder30-50km/sec.More
specifically, the D1, D2, and D3 disks have median vertical
gasvelocitydispersionsofinthecentral1kpcof35.1,46.3,
4.4 Toomre Q and Gas Velocity Dispersions and30.7km/s,respectively.Thisisconsistentwiththepre-
diction that c =f σ/2 for Q=1 (our c is slightly higher
t g t
Furtherinsightintothedynamicbehaviorofthecentralgas owing to the multi-component disk and order unity correc-
reservoir can be obtained from Figure 9 which shows the Q tions calculated in Faucher-Gigu`ere et al. 2013). The gas
parameterasafunctionoftimeforseveraldiskradii.TheQ velocity dispersion measured on smaller scales is similar in
value is evaluated via Q=2ctΩ/πGΣg using the zˆcompo- magnitude (i.e. 30-50 km/sec), but with significantly more
nentofthemassweightedgasvelocitydispersionasaproxy variability – particularly around burst episodes.
for the turbulent velocity dispersion. We adopt the verti-
cal, mass weighted, velocity dispersion because the radial
and azimuthal velocity dispersions are more subject to in-
4.5 Nuclear Kennicutt-Schmidt Relation
fluencebythecoherent,non-circularmotionofdenseclouds
in the disk plane. Moreover, the mass (rather than volume) Equilibrium analytic models have been used to explore the
weightedvelocitydispersionisbothmoreimportantforthe low efficiency of star formation that is observed in the KS
(cid:13)c 0000RAS,MNRAS000,000–000