MNRAS
511,
859–876
(2022)
https://doi.org/10.1093/mnras/stab3351
Advance
Access
publication
2021
No
v
ember
23
Sur
vi
v
al
and
mass
growth
of
cold
gas
in
a
turbulent,
multiphase
medium
Max
Gronke
,
1
,
2
‹
†
S.
Peng
Oh,
3
Suoqing
Ji
4
,
5
and
Colin
Norman
1
,
6
1
Department
of
Physics
&
Astronomy,
Johns
Hopkins
University,
Bloomberg
Center,
3400
N.
Charles
St.,
Baltimore,
MD
21218,
USA
2
Max
Planck
Institut
fur
Astrophysik,
Karl-Sc
hwarzsc
hild-Str
aße
1,
D-85748
Garching
bei
M
̈
unchen,
Germany
3
Department
of
Physics,
University
of
California,
Santa
Barbara,
CA
93106,
USA
4
Shanghai
Astronomical
Observatory,
Chinese
Academy
of
Sciences,
Shanghai
200030,
China
5
TAPIR
&
Walter
Burke
Institute
for
Theoretical
Physics,
Caltech,
Pasadena
CA
91125,
USA
6
Space
Telescope
Science
Institute,
Baltimore,
MD
21218,
USA
Accepted
2021
No
v
ember
15.
Received
2021
October
28;
in
original
form
2021
July
19
A
B
S
T
R
A
C
T
Astroph
ysical
g
ases
are
commonly
multiphase
and
highly
turbulent.
In
this
work,
we
investig
ate
the
survi
v
al
and
gro
wth
of
cold
gas
in
such
a
turbulent,
multiphase
medium
using
three-dimensional
hydrodynamical
simulations.
Similar
to
previous
work
simulating
coherent
flow
(winds),
we
find
that
cold
gas
survives
if
the
cooling
time
of
the
mixed
gas
is
shorter
than
the
Kelvin–
Helmholtz
time
of
the
cold
gas
clump
(with
some
weak
additional
Mach
number
dependence).
Ho
we
ver,
there
are
important
differences.
Near
the
survival
threshold,
the
long-term
evolution
is
highly
stochastic,
and
subject
to
the
existence
of
sufficiently
large
clumps.
In
a
turbulent
flow,
the
cold
gas
continuously
fragments,
enhancing
its
surface
area.
This
leads
to
exponential
mass
growth,
with
a
growth
time
given
by
the
geometric
mean
of
the
cooling
and
the
mixing
time.
The
fragmentation
process
leads
to
a
large
number
of
small
droplets
which
follow
a
scale-free
d
N
/d
m
∝
m
−
2
mass
distribution,
and
dominate
the
area
co
v
ering
fraction.
Thus,
whilst
survi
v
al
depends
on
the
presence
of
large
‘clouds’,
these
in
turn
produce
a
‘fog’
of
smaller
droplets
tightly
coupled
to
the
hot
phase
which
are
probed
by
absorption
line
spectroscopy.
We
show
with
the
aid
of
Monte
Carlo
simulations
that
the
simulated
mass
distribution
emerges
naturally
due
to
the
proportional
mass
growth
and
the
coagulation
of
droplets.
We
discuss
the
implications
of
our
results
for
convergence
criteria
of
larger
scale
simulations
and
observations
of
the
circumgalactic
medium.
Key
words:
hydrodynamics
–ISM:
clouds
–ISM:
structure
– galaxy:
kinematics
and
dynamics
– galaxies:
evolution
– galaxies:
haloes.
1
INTRODUCTION
Turbulent,
multiphase
gases
are
extremely
common
in
astrophysics.
We
find
them
in
the
interstellar
medium
(ISM),
circumgalactic
medium
(CGM),
intracluster
medium,
and
even
intergalactic
medium
(for
re
vie
ws,
see
e.g.
Meiksin
2009
;
Kravtsov
&
Borgani
2012
;
Klessen
&
Glo
v
er
2016
;
Tumlinson,
Peeples
&
We r k
2017
).
Un-
derstanding
the
dynamics
of
such
multiphase
gases
is
crucial.
It
is
essential
to
understand
the
phase
structure
of
such
gases
and
the
mass,
momentum,
and
energy
exchange
between
phases,
which
can
be
mediated
by
physical
mechanisms
such
as
radiative
cooling,
heating,
and
thermal
conduction.
In
a
static
medium,
this
exchange
is
relatively
well
understood
(Field
1965
;
Wa t e r s
&
Proga
2019
;
Das,
Choudhury
&
Sharma
2021
).
Ho
we
ver,
in
a
dynamic
situation,
hydrodynamical
instabilities
and
turbulence
can
mix
the
phases,
and
thus
can
quickly
dominate
the
mass,
momentum,
and
energy
transfer
rates
between
hot
and
cold
medium.
As
astrophysical
systems
are
rarely
static
–and
in
fact
often
highly
turbulent
–this
regime
is
of
great
interest.
Sev
eral
studies
hav
e
been
carried
out
with
the
broad
goal
of
characterizing
the
mass
transfer
rate
between
different
phases,
E-mail:
maxbg@mpa-garching.mpg.de
†
Hubble
fellow.
for
instance,
focusing
on
thermal
instabilities
(Sharma,
Parrish
&
Quataert
2010
;
McCourt
et
al.
2012
;
Sharma
et
al.
2012
;
Vo i t
et
al.
2015
),
turbulent
mixing
layers
(Begelman
&
Fabian
1990
;
Ji,
Oh
&
Masterson
2019
;
Fielding
et
al.
2020
;
Ta n ,
Oh
&
Gronke
2021
),
the
interaction
of
a
hot
wind
with
a
cold
cloud
(e.g.
Klein,
Mckee
&
Colella
1994
;
Mellema,
Kurk
&
R
̈
ottgering
2002
;
Scannapieco
&
Br
̈
uggen
2015
;
Armillotta,
Fraternali
&
Marinacci
2016
;
Br
̈
uggen
&
Scannapieco
2016
;
Schneider
&
Robertson
2017
;
Gronke
&
Oh
2018
;
Sparre,
Pfrommer
&
Vogelsberger
2019
;
Gronke
&
Oh
2020b
;
Kanjilal,
Dutta
&
Sharma
2020
;
Li
et
al.
2020
;
Abruzzo,
Bryan
&
Fielding
2021
;
Farber
&
Gronke
2021
),
or
more
complex
multiphase
geometries
(Kim
&
Ostriker
2015
;
Banda-Barrag
́
an
et
al.
2020
,
2021
;
Rathjen
et
al.
2021
).
In
this
paper,
we
focus
on
the
dynamics
of
a
cold
cloud
embedded
in
a
turbulent
hot
medium,
as
can
be
found
in
essentially
all
the
astrophysical
systems
mentioned
abo
v
e.
We
focus
on
the
regime
where
the
cooling
time
of
the
hot
medium
is
long.
We
do
not
study
thermal
instability
in
a
turbulent
medium,
which
can
also
lead
to
multiphase
structure
and
has
been
the
focus
of
several
previous
studies
(Hennebelle
&
P
́
erault
1999
;
Gazol
et
al.
2001
;
Kritsuk
&
Norman
2002
;
Saury
et
al.
2014
;
Kobayashi
et
al.
2020
).
The
paper
is
structured
as
follows.
We
will
first
set
the
stage
in
Section
2
where
we
discuss
the
analytic
expectations.
Then,
we
will
describe
in
Section
3
our
numerical
setup
before
we
present
its
results
in
Section
4.
We
discuss
our
findings
in
Section
5
before
we
© 2021
The
Author(s)
Published
by
Oxford
University
Press
on
behalf
of
Royal
Astronomical
Society
Downloaded from https://academic.oup.com/mnras/article/511/1/859/6433663 by California Institute of Technology user on 08 April 2022
860
M.
Gronke
et
al.
conclude
in
Section
6.
Videos
visualizing
our
results
can
be
obtained
at
http://
max.lyman-alpha.com/
multiphase-turbulence
.
2
ANALYTIC
CONSIDERATIONS
In
this
section,
we
re
vie
w
characteristic
length
scales,
and
discuss
expectations
for
cold
gas
growth
rates.
We
consider
a
cloud
of
size
r
cl
,
o
v
erdensity
χ
,
and
temperature
T
cold
embedded
in
a
hot
medium
with
temperature
T
hot
which
has
an
rms
turbulent
velocity
v
turb
=
M
c
s
,
hot
.
The
‘cloud
crushing
problem’
where
the
cloud
is
subject
to
a
uniform
wind
has
been
heavily
studied.
There,
hydrodynamical
instabilities
will
destroy
the
cold
gas
on
a
time-scale
of
the
order
of
the
Kelvin–Helmholtz,
Rayleigh–Taylor,
or
shock-crossing
time
of
the
cloud
t
cc
∼
χ
1/2
r
cl
/
v
turb
(Klein
et
al.
1994
).
From
our
previous
work
(Gronke
&
Oh
2018
),
we
expect
the
cloud
to
survive
if
t
cool
,
mix
<
αt
cc
,
(1)
that
is,
when
the
cooling
time
of
the
mixed
gas
t
cool
(
T
mix
,
n
mix
)
with
T
mix
∼
√
T
hot
T
cold
,
n
mix
∼
√
n
hot
n
cold
is
smaller
than
the
cloud
crushing
time.
In
this
‘wind-tunnel’
setup,
we
found
equation
(1)
to
hold
with
α
∼
1
(Gronke
&
Oh
2018
;
see
also
discussion
of
this
criterion
in
Kanjilal
et
al.
2020
and
references
therein)
but
as
we
will
see
the
dynamics
is
more
complex
in
a
turbulent
multiphase
medium
and
thus
we
leave
α
as
a
free
fudge
parameter
for
now.
Equation
(1)
can
be
rewritten
as
a
geometrical
criterion
that
states
that
clouds
larger
than
a
characteristic
size
survive
the
ram
pressure
acceleration
process
(Gronke
&
Oh
2018
,
2020a
),
R
>
r
crit
,
w
∼
v
wind
t
cool
,
mix
χ
1
/
2
α
−
1
≈
2
pc
T
5
/
2
cl
,
4
M
wind
P
3
mix
,
−
21
.
4
χ
100
α
−
1
,
(2)
where
T
cl
,
4
≡
(
T
cl
/
10
4
K),
P
3
≡
nT
/
(10
3
cm
−
3
K
),
mix
,
−
21
.
4
≡
(
T
mix
)
/
(10
−
21
.
4
erg
cm
3
s
−
1
),
M
wind
is
the
Mach
number
of
the
wind,
and
we
write
v
wind
=
c
s
,
wind
M
wind
∼
c
s
,
cl
M
wind
χ
1
/
2
.
This
implies
that
for
typical
galactic
conditions
clouds
larger
than
parsec
size
usually
fulfill
the
requirement.
The
scale
r
crit,
w
can
be
compared
to
the
characteristic
length
scale
of
cooling
induced
fragmentation
shatter
∼
c
s,
cl
t
cool,
cl
(McCourt
et
al.
2018
;
Gronke
&
Oh
2020b
)
to
give
r
crit
,
w
shatter
≈
50
M
wind
χ
100
(
(
T
cl
)
/
(
T
mix
)
0
.
5
)
α
−
1
,
(3)
where
quantities
are
e
v
aluated
at
v
alues
appropriate
for
our
fiducial
simulation.
Thus,
a
cloud
should
be
significantly
larger
than
the
characteristic
fragmentation
scale
to
exhibit
this
behaviour.
Note
that
we
use
r
crit,
w
for
the
‘wind
tunnel’
solution
to
differ
from
the
survi
v
al
scale
r
crit
in
a
turbulent
medium
investigated
in
this
work.
If
a
cloud
survives,
we
expect
the
cold
gas
to
grow
in
mass
due
to
continuous
cooling
which
we
characterize
by
̇
m
∼
v
mix
A
cl
ρ
hot
(4)
with
an
ef
fecti
ve
cold
gas
surface
area
A
cl
(see
below),
and
a
surrounding
hot
gas
density
ρ
hot
.
The
dynamics
of
the
mixing
layer
is
dominated
by
turbulence.
The
mixing
velocity
follows
(in
the
‘fast
cooling
regime’,
when
the
cooling
time
is
shorter
than
the
eddy
turno
v
er
time)
the
scaling
(Fielding
et
al.
2020
;
Gronke
&
Oh
2020a
;
Ta n
et
al.
2021
)
v
mix
∝
(
u
)
3
/
4
(
L
cold
t
cool
,
c
)
1
/
4
,
(5)
where
u
is
the
turbulent
velocity
in
the
cold
medium
and
depends
on
the
specific
setup.
1
Here,
t
cool
and
the
integral
scale
of
turbulence
L
cold
are
also
e
v
aluated
for
the
cold
medium.
In
summary,
generally
the
mass
transfer
between
the
phases
depends
on
the
cold
gas
surface
area
and
the
mixing
rate
per
unit
area.
The
former
dependence
is
intuitive.
In
combustion,
the
reaction
rate
is
proportional
to
the
reactant
surface
area;
likewise
here,
mass
growth
(where
hot
‘fuel’
is
converted
via
a
cold
‘reactant’
to
cold
‘ashes’)
is
proportional
to
the
cold
gas
surface
area.
In
our
wind-
tunnel
simulations,
the
‘ef
fecti
ve’
surface
area
A
cl
of
a
monolithic
cloud
in
equation
(4)
followed
a
simple
scaling
A
cl
∼
(
m
/
ρ
cl
)
2/3
.
Ho
we
ver,
in
turbulent
setup
where
cold
gas
is
constantly
fragmenting
to
smaller
scales,
the
surface
area
can
grow
faster
with
mass.
The
cold-hot
gas
interface
has
a
fractal
geometry
(Fielding
et
al.
2020
),
and
it
is
well-known
that
fractals
have
surface
areas
which
grow
faster
than
the
Euclidean
expectation
of
A
∝
V
2/3
,
as
inferred
from
e.g.
fractal
respiratory
organs.
For
instance,
one
possible
ansatz
is
that
A
∝
V
D
/3
,
where
2
<
D
<
3
is
the
fractal
dimension
(Barenblatt
&
Monin
1983
).
Another
limiting
case
is
when
large
clouds
continuously
fragment
down
to
some
scale
̃
r
.
If
so,
the
surface
area
is
simply
proportional
to
the
number
of
such
small
clouds,
so
that
A
∝
m
.
These
differences
are
important:
for
A
cl
∝
m
2/3
one
obtains
a
power
law
solution
of
the
form
m
(
t
)
=
(1
+
at
)
3
(with
a
being
a
combination
of
the
model
parameters)
whereas
for
A
cl
∝
m
it
follows
m
∼
m
0
exp
(
t
/
t
grow
)
with
t
grow
∼
χ
̃
r
/v
mix
where
̃
r
is
the
fragmentation
size.
We
shall
later
see
that
in
a
turbulent
scenario
the
mass
growth
rate
is
indeed
exponential.
How
the
turbulent
surface
area
A
T
scales
with
physical
parameters
is
a
thorny
problem
which
is
a
central
focus
of
the
turbulent
combustion
literature.
Many
different
scalings
have
been
proposed,
each
with
experimental
support
in
different
regimes
(see
discussion
in
Ta n
et
al.
2021
and
references
therein).
Ta n
et
al.
(
2021
)
side-
stepped
this
problem
by
arguing
that
a
turbulent
medium
should
have
a
net
cooling
time
given
by
the
geometric
mean
of
the
eddy
turno
v
er
and
cooling
time
of
the
cold
medium
(cf.
equation
(25)
of
Ta n
et
al.
2021
),
̃
t
cool
∼
(
l
u
t
cool
,
c
)
1
/
2
,
(6)
where
l
and
u
are
the
characteristic
length
scale
and
turbulent
velocity
in
the
cold
medium
,
respectively.
This
can
be
thought
of
as
the
geometric
mean
of
the
elastic
and
inelastic
time-scales
for
a
fluid
element,
similar
to
random
walk
time-scales
for
a
photon
in
the
presence
of
both
absorption
and
scattering,
which
gives
rise
to
an
ef
fecti
ve
optical
depth
τ
eff
∼
√
τ
abs
τ
scatter
or
an
electron
in
the
presence
of
both
Coulomb
(elastic)
and
atomic
(inelastic)
scattering,
which
gives
rise
to
the
Field
length
λ
F
∼
√
λ
e
λ
cool
.
This
ansatz
was
supported
by
detailed
simulations
of
turbulent
mixing
layers.
This
implies
a
mass
growth
time:
t
grow
≡
m
̇
m
∼
χ
̃
t
cool
∼
χ
M
−
1
/
2
(
l
shatter
)
1
/
2
(
l
L
box
)
−
1
/
6
t
cool
,
c
.
(7)
We
have
assumed
that
–as
seen
in
hydrodynamic
simula-
tions
–the
turbulent
pressure
is
continuous
across
the
interface
ρ
c
(
u
)
2
≈
ρ
h
v
2
turb
.
Thus,
the
cold
gas
turbulent
velocity
is
u
∼
v
turb
,
hot
(
l
)
/χ
1
/
2
=
M
c
s
,
c
(
l
/L
box
)
1
/
3
at
the
scale
l
.
The
scale
l
is
the
length
scale
characterizing
the
cold-hot
gas
interface.
Initially,
this
is
given
by
the
cloud-size;
later
on,
a
combination
of
mass
growth,
frag-
mentation
and
coagulation
renders
this
length
scale
more
ambiguous.
1
For
a
plane-parallel
shearing
layer
(Tan
et
al.
2021
)
found
u
∝
M
4
/
5
.
MNRAS
511,
859–876
(2022)
Downloaded from https://academic.oup.com/mnras/article/511/1/859/6433663 by California Institute of Technology user on 08 April 2022
Evolution
of
turbulent,
multiphase
gas
861
Figure
1.
Density
projections
of
individual
droplets
of
different
sizes
(from
r
d
∼
50
shatter
to
r
d
∼
500
shatter
from
top
to
bottom
row)
and
o
v
erdensity
χ
∼
100
in
a
turbulent
medium
with
M
∼
1.
The
small
droplet
is
dispersed
quickly
whereas
the
larger
one
manages
to
survive.
The
white
circle
indicated
the
initial
droplet
size.
Videos
showing
these
simulations
can
be
seen
at
http://
max.lyman-alpha.com/
multiphase-turbulence
.
In
principle,
at
late
times,
l
approaches
the
driving
scale
of
turbulence
L
box
.
This
is
true
of
mixing
layer
simulations,
where
cold
gas
fills
the
box.
Fortunately,
equation
(7)
implies
that
t
grow
only
has
a
weak
t
grow
∝
l
1/3
dependence
on
the
length
scales
l
.
We
have
found
that
with
our
limited
simulation
domain
and
the
boundary
conditions
we
use,
l
∼
r
cl
is
a
good
approximation
throughout.
Furthermore,
since
t
grow
=
m/
̇
m
∝
l
1
/
3
∝
m
1
/
9
is
roughly
constant
2
(and,
equi
v
alently
v
mix
),
we
expect
exponential
mass
growth
m
∝
exp(
t
/
t
grow
).
We
will
explore
the
validity
of
both
the
survival
estimate
and
the
simple
models
for
mass
growth
using
numerical
simulations.
2
The
cooling
time
of
gas
in
the
cold
medium
t
cool,
c
depends
only
on
pressure
and
metallicity,
which
we
assume
to
be
roughly
constant.
3
NUMERICAL
METHODS
For
our
hydrodynamical
simulation,
we
use
the
AT H E N A
++
code
(Stone
et
al.
2020
).
We
use
the
HLLC
Riemann
solver,
second-
order
reconstruction
with
slope
limiters
in
the
primiti
ve
v
ariables,
and
the
van
Leer
unsplit
integrator
(Gardiner
&
Stone
2008
).
We
implemented
the
Tow n s e n d
(
2009
)
cooling
algorithm
that
allows
for
fast
and
accurate
computations
of
the
radiative
losses.
We
adopt
a
solar
metallicity
cooling
curve
to
which
we
fitted
broken
a
power
law,
and
a
temperature
floor
T
floor
=
4
×
10
4
K.
3
Note
that
shatter
is
e
v
aluated
at
this
temperature
floor.
3
We
use
this
value
to
be
comparable
with
earlier
studies
(McCourt
et
al.
2015
;
Gronke
&
Oh
2018
)
and
to
ensure
the
gas
is
fully
ionized
since
we
assume
collisional
ionization
equilibrium
cooling
rates.
MNRAS
511,
859–876
(2022)
Downloaded from https://academic.oup.com/mnras/article/511/1/859/6433663 by California Institute of Technology user on 08 April 2022