Table Of ContentDraft version January 6, 2016
PreprinttypesetusingLATEXstyleemulateapjv.5/2/11
GALAXY OUTFLOWS WITHOUT SUPERNOVAE
Sharanya Sur1,2, Evan Scannapieco2, Eve C. Ostriker3
Draft version January 6, 2016
ABSTRACT
High surface density, rapidly star-forming galaxies are observed to have ≈ 50−100kms−1 line-of-
sightvelocitydispersions,whicharemuchhigherthanexpectedfromsupernovadrivingalone,butmay
arise from large-scale gravitational instabilities. Using three-dimensional simulations of local regions
of the interstellar medium, we explore the impact of high velocity dispersions that arise from these
disk instabilities. Parametrizing disks by their surface densities and epicyclic frequencies, we conduct
aseriesofsimulationsthatprobeabroadrangeofconditions. Turbulenceisdrivenpurelyhorizontally
6
and on large scales, neglecting any energy input from supernovae. We find that such motions lead
1
to strong global outflows in the highly-compact disks that were common at high redshifts, but weak
0
or negligible mass loss in the more diffuse disks that are prevalent today. Substantial outflows are
2
generated if the one-dimensional horizontal velocity dispersion exceeds ≈35kms−1, as occurs in the
n dense disks that have star formation rate densities above ≈ 0.1M yr−1kpc−2. These outflows are
(cid:12)
a triggered by a thermal runaway, arising from the inefficient cooling of hot material coupled with
J
successive heating from turbulent driving. Thus, even in the absence of stellar feedback, a critical
4 value of the star-formation rate density for outflow generation can arise due to a turbulent heating
instability. This suggests that in strongly self-gravitating disks, outflows may be enhanced by, but
] need not caused by, energy input from supernovae.
A
Subject headings: galaxies : evolution - galaxies : starburst - ISM : structure
G
.
h
1. INTRODUCTION tion of galaxy outflows remains a challenge.
p
One difficulty in understanding these processes stems
- Multi-wavelength observations reveal the existence of
o massive, galaxy-scale outflows of multiphase material, from the complex physics of the interstellar medium
r drivenfromrapidlystar-forminggalaxies(Heckmanetal. (ISM), where heating by UV photons, cosmic rays, and
t supernova shocks operate in combination with radiative
s 1990;Bomansetal.1997;Pettinietal.2001;Weineretal.
a cooling processes, leading to a multiphase, supersonic
2009; Martin et al. 2013; Rubin et al. 2014). Such out-
[ medium. Given the extremely short cooling times and
flows are thought to play a fundamental role in galaxy
theimportanceofsmallphysicalscales,simulationswith
1 evolution: enriching the intergalactic medium (IGM)
≈ parsec resolution are required to model the multi-
v (Songaila&Cowie1996;Simcoeetal.2002;Pichonetal.
phase gas distribution and expansion of supernova rem-
9 2003;Schayeetal.2003;Ferraraetal.2005;Scannapieco
nants directly (e.g., de Avillez & Breitschwerdt 2004;
5 et al.2006; Steidelet al. 2010;Martin etal. 2010), shap-
Hill et al. 2012; Hennebelle & Iffrig 2014; Walch et al.
6 ing the galaxy mass-metallicity relation (Dekel & Silk
0 1986; Tremonti et al. 2004; Erb et al. 2006; Kewley & 2014; Girichidis et al. 2015; Kim & Ostriker 2015). In
0 Ellison 2008; Mannucci et al. 2010), and affecting the particular, neither the correct mass of hot gas nor the
1. baryonic content and number density of dwarf galaxies correct momentum injected to the ISM by supernovae
can be captured unless both the Sedov and shell forma-
0 (Somerville & Primack 1999; Cole et al. 2000; Scanna-
tionstagesaresufficientlywellresolved(Kim&Ostriker
6 pieco et al. 2002; Benson et al. 2003). Outflows also aid
2015). Thus, galaxy-scale numerical simulations have
1 in eliminating small-scale random magnetic fields from
instead relied on a number of approximations, includ-
: galactic disks (Shukurov et al. 2006; Sur et al. 2007;
v ing: temporarily lowering the densities and cooling rates
Chamandy et al. 2014), thereby preventing the large-
i of heated particles (Gerritsen & Icke 1997; Thacker &
X scale dynamo from undergoing catastrophic quenching.
Couchman 2000; Stinson et al. 2006; Governato et al.
Yetdespitetheirimportance, understandingthedynam-
r 2007), using an empirical heating function (Mac Low
a ical and microphysical processes that control the evolu-
et al. 1989; Mac Low & Ferrara 1999), imposing an arti-
[email protected] ficial temperature floor (Suchkov et al. 1994; Tenorio-
1Indian Institute of Astrophysics, 2nd Block, Koramangala, Tagle & Munoz-Tunon 1998; Fujita et al. 2004), and
Bangalore560034,INDIA
implementing exaggerated momentum kicks (Navarro &
2School of Earth and Space Exploration, Arizona State Uni-
White1993;Mihos&Hernquist1994),andinthecaseof
versity,POBox876004,Tempe-85287,USA
3DepartmentofAstrophysicalSciences,PrincetonUniversity, feedback from active galactic nuclei, storing the energy
Princeton,NJ08544,USA
until it is sufficient to raise the temperature of the sur- sound speed of the warm medium and are greater than
rounding gas above a threshold value (Dalla Vecchia & the thermal sound speed of the cold medium. Thus any
Schaye 2012). characterization of the effective Q should depend on the
The connection to observations is further complicated total (thermal plus scale-dependent turbulent) velocity
by recent claims that instead of supernovae, outflows dispersion(Romeoetal.2010;Hoffmann&Romeo2012).
could primarily be driven by radiation pressure on dust In addition, ISM turbulence is driven by a combination
(Thompsonetal.2005;Murrayetal.2011;Hopkinsetal. of feedback from star formation and gravitational insta-
2011, 2012) or by non-thermal pressure contributed by bilities, and characterizing stability in realistic disks is
cosmicrays(Socratesetal.2008;Samuietal.2010;Uhlig therefore quite complex (Agertz et al. 2015). However,
et al. 2012; Booth et al. 2013; Hanasz et al. 2013; Salem both Milky-Way type galaxies simulations (Wada et al.
& Bryan 2014). In addition, the typically high ISM 2002; Agertz et al. 2009, 2015) and high-redshift galaxy
Reynolds numbers of Re ≈ 105 or greater, imply that a simulations(e.g.,Immelietal.2004;Ceverinoetal.2010;
completeunderstandingofthesemassiveoutflowscannot Geneletal.2012b,a)showthatatsufficientlylargescales
beobtainedbyignoringsmall-scaleturbulentstructures. gravitational instabilities promote an increase of the ve-
Interestingly, using a numerical sub-grid model for the locity dispersion, σ, until
unresolved turbulent velocities and length scales, Scan-
napieco & Bru¨ggen (2010) showed that it was possible Q ≡[σ2+c2]1/2κ/πGΣ≈1. (1)
eff s
to produce outflows of multiphase material arising from
simultaneous turbulent heating and radiative cooling in This implies that disks with high surface densities must
the disk. Notwithstanding the approximate nature of develop significant turbulent motions even if stars are
theschemeemployed,thisofferssupporttothestandard unable to stir the disks sufficiently to stabilize them, be-
picturethatmassive,galaxyscaleoutflowscanindeedbe causegravitationalinstabilitieswillleadtotheformation
produced from turbulence in the disk. of clumps moving at typical velocities σ ≈πGΣ/κ.
Onehintastothephysicsofgalaxyoutflowsmayliein Based on these ideas, Scannapieco et al. (2012) con-
the properties of the galaxies that host them. Outflows ducted simulations of turbulently stirred, radiatively
areobservedoverawiderangeofgalaxymasses,butina cooled media. These simulations modeled a local patch
smaller range of galaxy surface densities. Large outflows of the galaxy as a stratified medium in which turbulence
are ubiquitous in galaxies in which the star-formation was driven at a rate that matched the overall cooling
rate density per unit area exceeds a critical value of rate. At low velocity dispersions, such as occur in the
Σ˙cr ≈ 0.1M yr−1kpc−2 (Heckman 2002, 2003), while Milky Way, this configuration was stable for many dy-
th(cid:63)e ejection (cid:12)of material is more sporadic for Σ˙ < Σ˙cr namical times. On the other hand, the critical star-
(cid:63) (cid:63) formation rate density for galaxy outflows corresponds
values(Chenetal.2010). Recentobservationsalsoshow
to a gas surface density of Σ ≈100M pc−2, which as-
that disks with strong outflows are characterized by ve- g (cid:12)
locitydispersionsbetweenσ1D ≈50−100kms−1(Genzel suming typical values of Σ ≈ 2Σg, cs ≈ 10kms−1, and
et al. 2011; Swinbank et al.v2011). Such high velocities κ−1 ≈ 15Myr gives one-dimensional turbulent velocity
are difficult to obtain from supernovae acting alone, as dispersion of σ ≈ 35kms−1 to have Qeff ≈ 1. At these
high-resolutionsimulationsoftheISMwithawiderange highdispersions,Scannapiecoetal.(2012)discoveredthe
of SN rates show that such explosions can only drive ve- onset of a thermal runaway, where multiphase material
locities to ≈ 10−20kms−1 (Dib et al. 2006; Joung & moved upward from the disk and out of the simulation
Mac Low 2006; Joung et al. 2009; Kim et al. 2011, 2013; domain,implyingtheabsenceofastableequilibriumbe-
Shetty & Ostriker 2012; Hennebelle & Iffrig 2014; Gatto yond this critical value.
et al. 2015; Martizzi et al. 2015). Therefore, how does In this paper, we examine effects of turbulence in
one account for such high velocity dispersions in these Toomre-critical disks more closely by conducting three-
disks? dimensionalnumericalsimulationsofalocalpatchofthe
A possible solution lies in the gravitationally-driven ISM for a range of values of the total surface density
motions that occur in high surface density disks. The Σ and the epicyclic frequency, κ. Our goal is to better
Toomre stability criterion (Toomre 1964) relates the to- understand the varying environments in galaxies of dif-
tal disk surface density Σ, the epicyclic frequency κ, ferent surface densities, and the extent to which these
and the sound speed c in infinitesimally thin disks variations lead to the direct driving of galaxy-scale out-
s
that are marginally unstable to axisymmetric modes as flows and to conditions that are favorable to the driving
Q ≡ κc /πGΣ = 1, where G is the gravitational con- of outflows by supernovae. Motivated by the idea that
s
stant. Allowing for non-axisymmetric instabilities, mag- large-scale gravitational instabilities maintain a level of
netic fields, and interaction with a stellar disk increases turbulence in which Qeff ≈ 1, and that these instabili-
the critical Q value, while thick disk effects decrease it ties primarily involve in-plane motions, we apply driven
(e.g.,Romeo1992;Kim&Ostriker2001;Kimetal.2002; turbulence at a forcing level that results in a horizontal
Kim&Ostriker2007;Romeo&Falstad2013). Inthereal velocity dispersion
ISM, turbulent velocities are comparable to the thermal
σ ≈πGΣ/κ. (2)
H
2
Asweshallshow,thehorizontally-driventurbulencealso scalesandwillcascadetosmallerscalesvianonlinearin-
leadstoverticalmotions(ataloweramplitude), andthe teractions. Here, we chose a characteristic spatial scale
shocks from both horizontal and vertical motions heat for simplicity. We note that the adopted stirring scale
the gas. Thermal pressure gradients and vertical turbu- is comparable to the range where the turbulent power
lence combine to drive outflows from the disk. To fo- spectrum reaches its maximum in the shearing-box sim-
cus on the driving of turbulence purely by gravitational ulations of Kim & Ostriker (2007), who found that the
instabilities, we do not include any vertical mechanical power is flat above λ ≈8GΣ/κ2.
x
or thermal energy input from supernovae in this study, The simulations were conducted with the multidi-
such that all outflows obtained represent a lower limit mensional, grid-based (magneto)-hydrodynamic code
over which stellar processes will lead to additional con- FLASH(version4.2)(Fryxelletal.2000). WhileFLASH
tributions. is capable of incorporating dynamical grids of varying
Because our study is focused on obtaining a better resolution by virtue of the adaptive mesh-refinement
understanding of the role of gravitational instabilities (AMR) technique, we chose to perform our simulations
in changing the nature of the medium in which stellar on a uniform grid with the unsplit hydrodynamic solver
processes operate, we deliberately do not attempt to re- (Lee et al. 2009, 2012), in a box of size 3R in the x and
create a full model of the ISM, in which feedback from y directions and −3R to +3R in the vertical direction,
supernovae (Dib et al. 2006; Kim et al. 2011; Shetty & where R is defined by eq. (3). We adopted a 2562×512
Ostriker 2012), ionization fronts (Matzner 2002; Walch gridforthemajorityofoursimulations,and,forthepur-
et al. 2012), chemical transitions (Koyama & Inutsuka poseofaresolutionstudy,wealsoconductedsimulations
2000; Walch et al. 2011; de Avillez & Breitschwerdt atresolutionsof642×128,1282×256and5122×1024(see
2012), radiation pressure (Krumholz & Thompson 2012; Appendix). Furthermore, we chose periodic boundary
Sales et al. 2014), cosmic rays (Zirakashvili et al. 1996; conditionsinxandy and‘diode’boundaryconditionsin
Boettcher et al. 2013; Booth et al. 2013) and magnetic thez direction,whichallowforoutflowsofmaterialfrom
fields(Gresseletal.2008;Gentetal.2013)allplayarole. the simulation domain, but prevent inflows.
Although simulations that include this physics would be The initial conditions were characterized by two free
closer to real galaxies, they would also be much more parameters: the total matter surface density, Σ and the
difficult to interpret, as we could never be sure how to epicyclic frequency, κ. We further assumed that the gas
connect causes and effects unambiguously. In this sense, surface density was given by Σ = f Σ, where the gas
g g
simulationsthatstartfromsimplerinitialconditionsand fraction, f , was fixed at 1/2 for all our simulations. Be-
g
exploretheroleofafewfreeparametersarecomplemen- cause we do not know a priori what vertical distribution
tarytomorecomplexISMsimulationsinthesciencethey themediumwilltakeatlatetimes,weadoptedaninitial
are able to target. density distribution that maintained the desired surface
The structure of this paper is as follows. In Section density and approximated the vertical distribution as an
2 we discuss our numerical methodology. In Section 3 exponential profile:
we describe our results, focusing on outflow rates and
(cid:18) (cid:19)
phase distribution of the gas, turbulent properties, and Σ −|z| 1
ρ= exp , (4)
thecoolingandfree-falltimesinthemedia. Conclusions
2R/C R/C 1−exp(−3C)
are presented in Section 4.
where C is a ‘compression’ factor that relates R and the
2. NUMERICALMODELING initial scale height such that H(t = 0) = R/C, and the
Our simulations contain only four components: (i) the 1−exp(−3C)insuresthatthetotalsurfacedensitywithin
equations of compressible fluid dynamics; (ii) a continu- the finite domain is equal to Σ. In the gradual course of
ouslyupdatedaverageverticalgravitationalacceleration the simulation, and in the absence of any other feed-
to capture the evolution of the disk scale height in re- back processes, the scale height, H, the vertical velocity
sponse to thermal pressure and vertical turbulence; (iii) dispersion, σz, and the disk thermal structure adjusted
radiativecoolingofatomicgasintheoptically-thinlimit; automatically to the most compact and coldest distri-
and (iv) purely horizontally-driven turbulence that ap- bution available for a given choice of Σ and κ. If the
proximates the impact of gravitational instabilities in a turbulent motions were isotropic and thermal pressure
rotating disk in the absence of stellar feedback. support were minimal, then we would expect our simu-
A self-gravitating disk with velocity dispersion σ has lations to reach a state-state density distribution similar
a scale height H = σ2/(πGΣ). Taking this as a charac- to eq. (4) with C ≈1. However, if vertical turbulent mo-
teristic turbulent forcing scale λ and assuming that the tionsremainsmallerthanhorizontalones,asweshallsee
f
effectiveToomreparameterisunitysothateq.(2)holds, is the case, then we would expect much more compact
we have distributions,correspondingtoeq. (4)withlargerC val-
λ ≈R≡πGΣ/κ2, (3) ues. Forthesimulationspresentedinthispaper,wehave
f
used C = (5−7) at t = 0. In general, this is somewhat
for the turbulent stirring scale. In reality, power may more compact than the final steady-state distributions
be driven by gravitational instabilities over a range of described below, but because our gas distribution col-
3
lapses while the turbulence is developing, we find that
smallerchoicesofC leadtolongerdelaysinreachingthe
same quasi-steady state. In addition, in all our simula-
tions, we employed an ideal gas equation of state with
γ = 5/3 and an atomic mass of one, an initial constant
temperature of T =2×105K, and zero initial velocities.
Because our simulations were not conducted in a
shearing-boxinwhichlargescalegravitationalforcesare
opposed by centrifugal and Coriolis forces, and because
we do not incorporate sink particles to handle the col-
lapseofJeans-unstablegas,wecannotsimplyimplement
self-gravity in our simulations. Thus we adopted an ide-
alized approach. First, to include the impact of gravity
on the overall vertical profile of the gas, we computed a
vertical acceleration that was a time-dependent function
of z, and did not depend on x and y. The time de-
pendence ensures that the system finds its own vertical
scale-heightincourseofthesimulation,ratherthanhav-
ingtheuserguessascale-heightapriori,aswouldbethe
case if the user were to specify a gravitational profile at
runtime. Second, we added a horizontal stochastic driv-
ing term, to include the impact of large-scale horizontal
turbulent motions driven by gravity.
To compute the vertical gravity, at every time step we
Figure 1. Total atomic, ionic, and Bremsstrahlung cooling rate
calculated a vertical acceleration profile, g(z), directly (toppanel)andthecoolingtime(bottompanel)asafunctionofthe
from the average vertical mass profile as temperature for solar metallicity material. Note that the cooling
rate drops above T ≈ 2.5×105K, leading to an increase in the
(cid:90) |z| (cid:12)(cid:12)dΣ(cid:12)(cid:12) coolingtimescale.
g(|z|)=2πG (cid:12) (cid:12)dz, (5)
(cid:12)dz(cid:12) κ, the amplitude of the forcing was adjusted to get the
−|z|
desired value of σ = πGΣ/κ. Note that choosing a
H
where we approximate dΣ/dz as (1/f)dΣg/dz. Because differentvalueoftf wouldrequirere-adjustingtheforcing
this acceleration is a function of height and not of hor- amplitude so as to obtain the correct value of σ .
H
izontal position, it does not directly drive turbulence. The majority of the energy input from turbulent driv-
To implement turbulent driving as would be induced by ing in our simulations leaves the system through radia-
large-scale gravitational instability, we used a version of tive cooling, as it does in a real galaxy. Following the
the stirring module in FLASH, adapted from one of our prescription of Gray & Scannapieco (2010), we imple-
recent studies (Sur et al. 2014b), which models turbu- mentedatomic,ionic,andBremsstrahlungcoolinginthe
lent driving random motions as a stochastic Ornstein- opticallythinlimit,assuminglocalcollisionalequilibrium
Uhlenbeck (OU) process (Eswaran & Pope 1988; Benzi with,
et al. 2008). Specifically, this corresponds to a Gaussian ρΛ(T,Z)
random vector field, F, with an exponential temporal E˙cool =(1−Y)(1−Y/2) (µm )2 . (8)
correlation, t , in the momentum equation. p
f
We drove turbulence in the range of wavenumbers 2≤ Here E˙ is the radiated energy per unit mass, ρ is
cool
|k|L/2π ≤3, such that the average forcing wavenumber the density in the cell, m is the proton mass, Y is
p
was kfL/2π (cid:39)2.5. Here L=3R is the horizontal extent the helium mass fraction, µ the mean atomic mass, and
(cid:113)
of the box and the wave vector |k| = (k2+k2+k2). Λ(T,Z) is the cooling rate as a function of temperature
x y z
and metallicity. For the cooling source terms, we im-
This corresponds to a turbulent driving scale,
plemented a sub-cycling scheme (Gray & Scannapieco
L =2π/k =1.2R. (6) 2010), such that T and Λ(T,Z) were recalculated ev-
f f
ery time E /E > 0.1. This is equivalent to an inte-
cool
Since we wish to restrict driving of turbulence to a two- gralformalismthatassumesaconstantdensityovereach
dimensional plane (to simulate the effect of self-gravity hydrodynamic time step (Thomas & Couchman 1992;
withinadisk),wemodifiedourroutinestoforceonlythe Scannapieco et al. 2001). Furthermore, in order to help
horizontal component of the velocity. The correlation thesimulationsreachquasisteady-stateconditionsmore
time of turbulence, tf, was chosen to be ≈0.1ted, where quickly, cooling is disabled for the first 0.5R/σH time of
t =L /σ =1.2κ−1, (7) the simulation, to avoid excessive vertical collapse while
ed f H the disk gains an initial level of turbulence.
is the eddy-turnover time. For a given value of Σ and Thecoolingrateswerecalculatedusingthetablescom-
4
piled by Wiersma et al. (2009) making the simplifying 3. RESULTS
approximation that the metallicity is always solar and
3.1. Outflow Rates
that the abundance ratios of the metals always occurs
in solar proportions. In the upper panel of Figure 1, we As our goal is to obtain a better understanding of the
show the total radiative cooling rate as a function of the turbulentISMasafunctionofgalaxyproperties,wecon-
temperature. The lower panel shows the behavior of the ducted a suite of simulations with different values of the
local cooling timescale defined as the ratio of the energy gas surface density, Σg, and the epicyclic frequency, κ.
per unit volume to the radiative cooling rate per unit In Table 1, we present a summary of the properties of
volume, these simulations, including the run parameters (Σg, κ),
the key quantities for each model that can be directly
tcool ≡1.5nkBT/(ρE˙cool)∝ρ−1, (9) derived from these parameters - the forcing scale Lf,
box size above the midplane, z = 3R, the resolu-
max
where kB is the Boltzmann constant and n = ρ/µmp is tion dz = zmax/256, the horizontal velocity dispersion
the number density of the gas. σH, and the escape velocity ves, and the quantities that
Following the nomenclature introduced by Wolfire canonlybemeasuredfromthefullsimulations-thever-
et al. (1995, 2003), much of the interstellar medium tical velocity dispersion σz, the scale-height H, and the
can be described as comprised of a cold neutral medium mass loss rate Σ˙ . The idealized nature of our simula-
g
(CNM)withT ≈100K,awarmneutralmedium(WNM) tions enables us to probe a wide parameter space, with
withtemperaturesintherange,T ≈6×103−104Kanda Σ ∈ [50,500]M pc−2 and κ−1 ∈ [6.5,30]Myr. Note
g (cid:12)
hot ionized medium (HIM) with T ≥106K. As Figure 1 however that for runs with κ−1 = 30Myr, simulations
shows, we have only considered cooling in the tempera- with an initial Σ >75M pc−2 result in box sizes that
g (cid:12)
turerangeT ≈103−107Kanduptoamaximumdensity are larger than the typical size of entire disk galaxies,
of n =6.0×107cm−3. We note that the colder parts and we have therefore chosen to omit such simulations
max
of the ISM occupy a small fraction of the volume (only from our study.
a few percent), and while they are important for star Our motivation for using a broad range of values is to
formation, they interact much less strongly with the hot capture the properties of high surface density galaxies
medium than the warm gas that occupies most of the over a broad range of masses and redshifts. Several at-
volume. Since our goal is to capture the formation of tempts have been made to make predictions about the
galaxy outflows which are mainly comprised of hot gas evolution of the galaxy size with redshift (see Shapley
in the temperature range T ≈106−107K, we to lowest 2011; Silk & Mamon 2012; Conselice 2014, for recent re-
order neglect the inclusion of the cold dense phases of views), and observational studies have shown that dis-
the ISM in our simulations. We therefore do not include tant galaxies are more compact than those of the same
low-temperature cooling instead placing a temperature mass in the nearby Universe (e.g., Daddi et al. 2005;
floor at 1000 K throughout the simulations. This there- Trujillo et al. 2007; Buitrago et al. 2008; van Dokkum
fore amounts to adopting a simplified description of the et al. 2010). From these studies, the dependence of
multi-phase ISM, which as we shall show in later sec- galaxy size with redshift selected at a constant stellar
tions is sufficient for the physical problem at hand. In a mass can be roughly characterized as a power law of the
testcase,weloweredthetemperaturefloorto300Kand form R ∝ (1 + z)−α, although Ferguson et al. (2004)
found the mass outflow rates to remain unchanged (see find R ∝ H−1(z) from z ≈ 1 − 5, where H(z) is the
Appendix). Hubble parameter. Assuming, for simplicity, a power
At temperatures T ≈ 104 − 106.5K, cooling re- law relation between the galaxy size and redshift, it can
sults mainly from line emission, while cooling due to be shown that the epicyclic frequency for a fixed mass,
bremsstrahlung (free-free emission) becomes important scales as κ ∝ Σ1/2(1+z)α/2, while the surface density
at temperatures (cid:62) 107K. Below 104K, collisions are at fixed mass scgales as Σ ∝ (1+z)2α. These scalings
not energetic enough to excite atomic transitions, lead- g
imply that star-forming galaxies at higher redshifts are
ing to a drop in the cooling rate. The atomic cooling
both more compact (i.e., have larger values of Σ ) and
rate attains a peak value of Λ ≈ 10−21ergscm3s−1 in g
are more rapidly rotating (i.e., have lower values of κ−1)
the temperature range T ≈ (1−2.5)×105K. In this
thantheirlocalcounterparts. Asweshowbelow,thisal-
regime, the cooling time is roughly constant as seen in
lows us to make some quantitative statements about the
thelowerpaneloftheplot. Thepeakat T ≈2.5×105K
role of gravity-driven turbulence in powering outflows in
is dominated by line emission from metal ions, whose
distant galaxies. This is all the more important as the
atomic energy levels are easily excited by collisions at
integratedstar-formationratepeaksatz ≈2(seeMadau
this temperature. However, beyond T ≈ 2.5 × 105K,
& Dickinson 2014, for a review).
mostoftheatomsbecomefullyionized, theeffectiveness
In the top panel of Figure 2, we show the variation of
of the line cooling decreases, and the cooling rate drops,
the averaged gas mass loss rate as a function of Σ for
leading to a gradual increase in the cooling time. What g
different values of κ−1. The filled symbols denote runs
thisimpliesforthenatureofthemultiphasemediumwill
be discussed in the next section. where we find outflows with Σ˙g (cid:62) 0.02M(cid:12)yr−1kpc−2,
5
Simulation Σg κ−1 Lf zmax dz H σ¯H1D σ¯z ves Σ˙g
Name [M(cid:12)pc−2] [Myr] [pc] [kpc] [pc] [pc] [kms−1] [kms−1] [kms−1] [M(cid:12)yr−1kpc−2]
S500K6.5* 500 6.5 723 1.80 7.03 7.0 68 5.0 309 0.22
S250K6.5* 250 6.5 362 0.9 3.50 9.4 35 5.5 155 0.073
S150K6.5* 150 6.5 217 0.53 2.10 17 21 6.3 92 0.04
S250K10* 250 10 856 2.13 8.32 11 50 5.5 238 0.06
S150K10* 150 10 514 1.27 4.96 20 32 6.0 142 0.03
S100K10 100 10 342 0.85 3.32 27 21 6.2 95 0.017
S150K20* 150 20 2054 5.12 20.0 22 60 5.0 285 0.03
S100K20* 100 20 1369 3.40 13.2 26 53 5.0 190 0.02
S75K20 75 20 1027 2.56 10.0 35 30 6.0 143 0.009
S50K20 50 20 685 1.70 6.64 53 21 6.3 95 0.004
S75K30 75 30 2311 5.74 22.4 33 45 5.8 214 0.005
S50K30 50 30 1541 3.83 15.0 50 30 6.0 143 0.004
Table 1
Summaryofthesimulationrunsatauniformgridresolutionof2562×512. Thedifferentcolumnsfromthelefttotherightare-
1)Simulationname,2)Gassurfacedensity,3)Inverseoftheepicyclicfrequency,4)thehorizontalforcingscale,L =1.2R,5)box
f
sizeabovethemidplane,6)theresolution,7)verticalscaleheightobtainedfromequation(13),8)time-averagedmass-weighted
1Dhorizontalvelocitydispersion,9)time-averagedmass-weightedverticalvelocitydispersion,10)escapevelocityofthegasas
definedbyequation(11),and11)themasslossrate. Simulationsthatshowpersistentoutflowsformanyeddy-turnovertimesare
markedbyan*.
persistingovermanyeddy-turnovertimesindisks. Below
this threshold, outflows appear to be sporadic or absent
with negligible values of Σ˙ (denoted by open symbols),
g
asobservedinlowersurfacedensitygalaxies(Chenetal.
2010). Inpanel2(b)weshowthevariationoftheratioof
η =Σ˙ /Σ˙ , whereas the inset figure shows the variation
g (cid:63)
ofthestar-formationratealone;forthispurposewehave
adopted the Kennicutt (1998) fit to the star-formation
rate density as a function of (Σ ,κ−1)
g
Σ˙ ≈0.017Σ Ω=0.012Σ κ, (10)
(cid:63) g g
where Ω is the angular velocity and in the last equality,
√
we have assumed κ = 2Ω, as appropriate for a flat
rotation curve. The symbols (open and filled) follow the
convention of those in panel 2(a).
We recall that large-scale outflows are most likely
to occur in systems in which the star-formation rate,
Σ˙ (cid:62)0.1M yr−1kpc−2. This implies that for each pair
(cid:63) (cid:12)
ofvaluesofΣ andκ−1,massiveoutflowsareexpectedfor
g
all runs except for: S50K30, S75K30, S50K20, S75K20,
S100K20 and S150K20. Comparing with the gas mass
lossrateplottedinpanel2(a),wefindthatrunsS50K30,
S75K30, S50K20 and S75K20 show very weak outflows
(i.e. below the threshold value of 0.02M yr−1kpc−2).
(cid:12)
Figure 2. (color online). The upper panel shows the averaged Among the other two runs, the mass loss rates in
gas mass loss rate as a function of Σg for different values of κ−1. S100K20 and S150K20 are on and above the threshold
Filledsymbolsdenoterunswherethecomputedgasmasslossrate
(cid:62) 0.02M(cid:12)yr−1kpc−2 (horizontal dashed line). The lower panel value respectively. Run S100K10 lies above the ‘Heck-
showstheratioη ofthegasmasslossrateandthestarformation man limit,’ but shows only a weak outflow. The ratio
rate(wherethestarformationrateisobtainedfromequation10). (η) of the gas mass loss rate and the star formation rate
Theinsetpanelshowstheadoptedvariationofthestarformation
for runs with κ−1 = 6.5 and 10Myr stays nearly con-
rateasafunctionofthegassurfacedensitywhilethedashedline
denotes the ‘Heckman limit’. The filled (open) symbols in these stant, η ≈ 0.1. However, for κ−1 = 20Myr, η varies
panelsfollowtheconventionofthoseinpanel2(a). weaklywithΣ forrunsbothaboveandbelowtheHeck-
g
manlimit. Thetrendforκ−1 =30Myrisuncleardueto
6
our choice of simulation parameters of Σ . our initial conditions, and so we avoid making measure-
g
In Figure 3, we show a three-dimensional rendering mentsordrawingconclusionsfromthisphaseofthesim-
of the density at four different times, in a run with ulations.
Σ = 250M pc−2 and κ−1 = 10Myr. These volume On the other hand, after the initial phase has passed,
g (cid:12)
renderings show the evolution of an outflow with plumes the simulations settle into a quasi-steady state that de-
moving gradually both upwards and downwards from pends almost purely on Σ and κ. We thus evaluate the
g
the midplane, eventually leaving the simulation domain. time-averagedvaluesinTable1onlyinthisquasi-steady
Notethatthequestionofwhethertheejectedmaterialis state. Bythistime,themediumintheS250K10runset-
able to leave the simulation domain depends on whether tlesintoaroughlyconstantmass-weighted1Dhorizontal
the outflow velocity exceeds the escape velocity, which velocitydispersionofσ¯1D =50kms−1 whilethemedium
H
variesfromsimulationtosimulationinourstudy. Specifi- intheSK50K30runreachesσ¯1D =30kms−1. Whilethe
H
cally, given theverticalextentof the box abovethe mid- higher surface density run is somewhat hotter than the
plane and the gas surface density, the critical velocity moderate surface density case, in both runs the sound
requiredtoescapethegravitationalpotentialofthehost speed is well below the average velocity dispersion. This
galaxy is given by implies that the turbulent motions are supersonic in na-
ture. In fact, the mass weighted rms Mach number ≈ 6
v =(cid:112)2gz ≈(cid:18) 103 (cid:19) (cid:18) Σg (cid:19)1/2 (cid:18) z (cid:19)1/2,for the higher surface density run and ≈ 4 in the lower
es max kms−1 100M pc−2 kpc surface density run, signifying the occurrence of strong
(cid:12)
(11) shocksinbothruns. Thus, althoughthedrivinginthese
where g = 4πGΣ is the acceleration due to gravity simulationsispurelyhorizontal,italsocausessubstantial
g
(assuming Σ = 2Σ ) and z is the vertical extent of pressureanddensityperturbations,whichvarybothhor-
g max
the simulation volume. However, since the vertical ex- izontallyandvertically. Theseleadinturntoverticalve-
tent is three times the driving scale R, which in turn locityfluctuations,whichareroughlyconstantintimeat
would be related to the typical radial scale length, Rl σ¯z = 5.5kms−1 in the S250K10 case and σ¯z = 6kms−1
(as GπΣR2 ≈Ω2R , which implies R ≈πGΣ/Ω2 ≈R) in the S50K30 case. This shows that the presence of
l l l
in a more complete model disk evolution, we can rea- supernova energy input is not required for high-density
sonably expect that in a real disk as material travels to disks to develop significant turbulent motions perpen-
heights approaching z , g(z) would be falling rapidly, dicular to the plane of the galaxy, provided large-scale
max
resulting in escape velocities similar to those in our sim- self-gravitating instabilities are able to maintain hori-
ulations. zontal motions of several tens of kms−1. Quite remark-
ably, as Table 1 shows, the value of the time averaged
mass-weightedverticalvelocitydispersionisconsistently
3.2. Time Series, Vertical structure, and Thermal
within the range from 5−6.3kms−1 in our simulations,
Phases
which we discuss in more detail below.
Beforewedelvefurtherintothedependenceoftheout- A closer look at these plots also reveals that for the
flows on the disk properties, we show the time evolu- S250K10 run, a distinct outflow occurs throughout the
tion of two runs that illustrate the range of behaviors full quasi-steady state evolution (see panel 4c), persist-
seen in our simulations. Figure 4 shows the evolution ing over many eddy-turnover times until the simulation
of S250K10, a high surface density, rapidly-rotating case is terminated at t = 330Myr. This is also evident from
with Σg = 250M(cid:12)pc−2 and κ−1 = 10 Myr, and Fig- the two-dimensional x−z slices of the density, vertical
ure 5 show the evolution of S50K30, a moderate sur- component of the velocity, and temperature in Figure 6.
face density, slowly-rotating case with Σg = 50M(cid:12)pc−2 Hereweseethat,incourseoftheevolution,thecombined
and κ−1 = 30 Myr. In both figures, we plot the mass- action of turbulent heating, radiative cooling, and gravi-
weighted on√e-dimensional horizontal velocity dispersion, tationalcollapseleadstostrongdensityandtemperature
σ1D = σ / 2, the vertical velocity dispersion, σ , the contrasts throughout the simulation. Furthermore, the
H H z
mass-weighted rms Mach number, and the gas mass sur- slicesofthedensityplottedinthefirstrowofthisfigure,
face density. In both of these runs, the initial time evo- together with the vertical velocity plotted in the second
lution is marked by a transient phase, which lasts for 60 row, clearly show an outflow of material.
Myr and 200 Myr respectively (denoted by the dashed At this point, it is worth describing in greater de-
vertical lines), after which the system reaches a steady tail how such outflows are generated. To answer this
state. In this stage, the action of gravity in the absence question, we return to the total atomic, ionic, and
of turbulent support causes the material to move down- Bremsstrahlung cooling curve in Figure 1. The temper-
wards towards the mid plane and become compressed. ature slices in Figure 6 show that by t = 311Myr, the
Then as turbulence develops, pressure support also in- initialconstanttemperaturedistributionhasevolvedinto
creases, puffing the layer back up, and causing a rapid amultiphasestructure. Inparticular,materialincertain
expansionthatdrivesasmallfractionofthematerialout regions close to the mid plane is heated to T ≈ 106K
of the simulation volume. During these rearrangements, and higher, while some patches of gas are still close to a
our results are strongly dependent on the particulars of few times 105K. At the same time, Figure 1, shows that
7
Figure 3. (color online). A 3D rendering of the density in the run S250K10, showing the gradual evolution of the outflow at four
differenttimesstartingfromt=312to319Myr(lefttoright). Forclarityofthecolorcontrast,wehaverestrictedthedensityrangefrom
logρ=−30.50to−24[gcm−3]. Plumesofgasmoveawayfromthediskmidplaneandleavethesimulationdomainthroughthevertical
boundaries.
beyond T ≈ 2.5×105K, the radiative cooling rate de- and ∂lnΛ/∂lnT (cid:46) 0 throughout much of the medium.
creases. Thisimpliesthattheseregionswouldtakemany This will result in an unstable state in which rapid, run-
dynamical times to cool to the mean temperature of the away heating of the hot, low density medium will lead
medium. But by then, successive heating resulting from to the significant removal of the gas from the simulation
turbulent driving would further heat the material. domain (see also Sharma et al. 2012; Scannapieco et al.
This behavior is also expected from a linear stability 2012). Thus, from t=311Myr onwards, Figure 6 shows
analysis of a temperature perturbation in a medium in the continuing impact of this hot gas, which pushes its
which the cooling rate per unit mass is ∝ ρΛ(T), and way outward both through the top and bottom bound-
the heating rate per unit mass is roughly constant, as is aries.
the case with turbulent heating (Pan & Padoan 2009). In contrast, for a slowly-rotating, moderate surface
In this case, the the condition for such a perturbation to density disk, the temperature slices in Figure 7 show
grow exponentially (at constant pressure) is simply that the majority of the material away from the mid-
plane is at temperatures (cid:54) 2.5×105K. In this phase,
∂lnΛ
<1, (12) ∂lnΛ/∂lnT ≥ 1 and the cooling time is roughly con-
∂lnT stant or dropping as a function of temperature, leading
to a stable configuration with less significant mass loss
(Field 1965; Defouw 1970; McCourt et al. 2012). In a
through the vertical boundaries.
steadily-heatedthermally-unstablemedium,highdensity
Such multiphase outflowing material is also detected
regions condense and cool, while low density regions ex-
in observed high surface density starbursting galaxies,
pandandheat. Whenσ islarge,thetemperatureishigh
8
Figure 4. (color online). Time evolution of the mass-weighted Figure 5. (coloronline). SameasinFig.4forrunSK50K30with
one-dimensionalhorizontalvelocitydispersion(black)andthever- Σg =50M(cid:12)pc−2, κ−1 =30Myr, and ted =36Myr. In this case
tical velocity dispersion (red) in panel 4a), the rms Mach num- theseisnegligiblemasslossovera400Myrlongperiod.
ber (blue) in panel 4b), and the gas mass surface density (dark
green) in panel 4c) for run SK250K10 with Σg = 250M(cid:12)pc−2, we find the existence of a horizontal, high density tail at
κ−1=10Myr,andted=12Myr. Thedashedverticallinedepicts T ≈ 104K; a feature commonly seen in phase diagrams
thedurationoftheinitialtransientphase. ThevaluesofσH1D,σz, in ISM simulations. In the S50K30 case, corresponding
and Σ˙g presented in Table 1 are thus computed from 60 Myr on- to a slowly rotating, moderate surface density galaxy, a
wards. Noticethatthereissignificantmasslossattimeswellafter
majority of material is below the threshold temperature
theinitialtransientphase.
of2.5×105K. Thisrunhasonlynegligiblemassloss,and
which contain both 107−108K emitting material (Mar- the T −vz plots show that the vertical component of the
tin1999;Strickland&Heckman2007,2009)and≈104K velocity corresponding to the gas near the top boundary
material detectable in a variety of optical and near is 50−100kms−1, as compared to the escape velocity of
UV absorption and emission lines (Pettini et al. 2001; 143 km/s.
Tremonti et al. 2007). Certainly, in real galaxies, much On the other hand, in both S250K10 with σ¯1D = 50
H
of the hottest of these phases is infact material from su- kms−1 and S500K6.5 with σ¯1D = 68 kms−1, more of
H
pernovae. On the other hand, the presence of ≈ 106K the low density gas above the disk mid plane is at con-
gas may provide pockets in which the cooling of mate- siderably higher temperatures than the temperature at
rial shocked by supernova remnants is minimal, making whichthecoolingtimeisequaltotheeddyturnovertime.
such galaxies particularly prone to outflows. Or in other Due to the longer cooling time, this gas becomes ther-
words, gravitationally-driventurbulenceislikelytoboth mally unstable and leads to a rapid, runaway heating
drive a baseline outflow as well as provide an ISM dis- resulting in an eventual expulsion of the gas from the
tribution in which supernova driving is particularly effi- host galaxy, dragging along additional material with it.
cient. Thus, such high density, faster spinning disks are more
Figure 8 shows the phase distributions of the gas for likely to host galaxy outflows than their slow rotating
three different runs, including the temperature-density lessdensercounterparts. Themass-weightedphaseplots
distribution over the whole simulation volume and the of the temperature and the vertical component of the
temperature-velocity distribution of the gas at high lati- velocity for S250K10 and S500K6.5 show that the out-
tudes. In the T −n plots, T =2.5×105K is denoted by flowing gas is made up of a combination of phases with
the horizontal line, above which the cooling rate starts most of the mass in the hot phase at T ≈ 106 −107K
to decrease. The vertical lines in the plots denote the which also drags along some relatively low temperature
densities at which the cooling time at T = 2.5×105K gas as it escapes the galaxy. In both cases the vertical
is equal to the eddy turnover time. To the left of this velocitiesrangefrom≈100kms−1 toclosetotheescape
line and for T >2.5×105K, the cooling time is greater velocity, which is 238 kms−1 in the S250K10 case and
thant ,whiletotheright,t <t . Inallthreeplots, 309 kms−1 in the S500K6.5 case.
ed cool ed
9
Figure 6. (color online). Two-dimensional slices in the x−z plane showing the time evolution of the logarithm of the density (upper
row) and the vertical component of the velocity (middle row) and the logarithm of the temperature (lower row) in a local patch of the
turbulent ISM for a run with Σg =250M(cid:12)pc−2, κ−1 =10Myr and ted =12Myr. The total vertical extent of the simulation domain is
2z=4.26kpcwhilethehorizontalextentisx=y=2.13kpc.
3.3. Turbulent Parameters and Vertical Scale Heights mass loss rate with the one-dimensional horizontal ve-
locity dispersion, we find that with the exception of
To further explore the connection between the mass
S75K30,strongoutflowsareubiquitousinsystemswhere
loss rate and the structure of the disk, we measured
σ¯1D (cid:62) 32kms−1. This is in agreement with the earlier
how turbulent velocities and vertical scale-heights vary H
numerical (Scannapieco et al. 2012) and analytical pre-
with our input parameters, Σ and κ. In Figure 9, the
g dictions(Scannapieco2013),whicharguedthatprovided
filled (open) symbols follow the classification of disks
turbulentvelocitiescanbemaintainedabovethisthresh-
hostingoutflowsabove(below)0.02M yr−1kpc−2,sim-
(cid:12) old value (about 35kms−1 in their paper), the material
ilar to those in panel 2(a) of Figure 2. Comparing the
isshockedintothethermallyunstableregime,leadingto
10
