Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://star.arm.ac.uk/preprints/401.ps Äàòà èçìåíåíèÿ: Wed Sep 17 19:19:00 2003 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 21:34:31 2012 Êîäèðîâêà: Ïîèñêîâûå ñëîâà: molecular cloud

Mon. Not.
R. Astron. 339, 133--147 (2003)
The instability fast shocks molecular clouds
Michael Smith and Alexander Rosen
Armagh Observatory, College Hill, Armagh BT61 9DG
Accepted 2002 October
7. Received 2002 September original form October
ABSTRACT
report
on discovery moderately shocks dense molecular clouds with low
transverse magnetic fields likely
to
be unstable. The instability
is triggered promoted
cooling that results the formation
of carbon monoxide water molecules extended
warm shock section. Numerical methods employed demonstrate the absence
magnetic fields, instability regime restricted
to densities above
0
, velocities
between 30--70 km
s
, and
C abundances above #10
, that cooling from reforming
molecules dominates warm without being suppressed ultraviolet dissociation. The
result either quasi­periodic chaotic collapse re­establishment
of warm shock
a typical time­scale
of
6 with variations shorter time­scales and
changes
in period being possible. Infrared emission from unstable region, including
2 lines, exhibit orders magnitude variability. Atomic H# display constant
fluxes undergo radial velocity variations.
words: hydrodynamics instabilities molecular processes
-- shock waves ISM: lines
bands ISM: molecules.
NTRODUCT
Fast shock waves alter dynamical, physical chemical
erties dense interstellar clouds (Hollenbach &McKee 1979).
shocks driven internally externally
a variety sources,
including supernova blast waves, protostellar jets, cloud--cloud
lisions, stellar winds planetary nebula. define
a
one capable destroying molecules path.
a dense cloud,
many
of molecules subsequently reform
in cooling
layer compressed accelerated material (Hollenbach &McKee
1989; Neufeld &Dalgarno 1989a). compressed cloud may
form stellar systems fragment disperse. Hence,
behaviour shocks critical
to distributions
stars interstellar gas.
Here investigate time­dependent properties shocks
through one­dimensional (1D) numerical simulations.
ture signatures steady fast shocks was studied
in
by
Hollenbach
& McKee (1989), Neufeld
& Dalgarno (1989a,b)
Wolfire K˜onigl (1991).
a shock
is unstable, however,
reconsider only observed emission properties
of radiative
shocks also transfer turbulent energy cloud,
fragmentation destruction cloud, result­
ing shock­induced molecular abundances and accuracy
of
numerical simulations not resolve individual shocks.
Hydrodynamic shocks destroy molecules
at shock
possess speeds excessvadjust #23
-1 (Kwan
E­mail: mds@star.arm.ac.uk (MDS); rar@star.arm.ac.uk
Hollenbach
& McKee 1980; Smith 1994). dense molecular
shielded ultraviolet (UV) radiation, shocks continuous
(C­type)
in which
a fraction leads drift
of the magnetic
through neutral case, molecules largely survive
in shocks with speeds
to #30
s
if magnetic field
parallel (Smith 1993)
s
is transverse
(Draine, Roberge &Dalgarno 1983). Below these speeds, shock
cooling layer overlap. Above these speeds, however,
shock classified fast molecules abruptly disso­
ciated shock front, followed
a wide zone radiative
cooling. demonstrate moderately shocks
in dense molec­
clouds are thermally dynamically unstable. Thermal insta­
bilities small scales expected cooling layer follows
shock front provided cooling unit mass increases
as temperature (McCray, Kafatos Stein 1975). This
is
cause
a small region within the with
a lower temperature than
surroundings faster, resulting
a growing density contrast
Smith 1989).
To dynamically unstable, cooling must increase
sufficiently rapidly temperature falls (Langer, Chanmugam
Shaviv 1981; Chevalier
& Imamura 1982).
a certain time­scale,
a cooling shocked layer collapse simultaneously
as shock
shrinks back, reducing shock strength. shock
velocity promotes faster cooling, which leads even faster
collapse
of cooling layer. Thus,
a resonance occurs and entire
collapse `catastrophically'. The lowered shock velocity
leads
a lower pressure. Hence,
a resisting medium confin­
can decelerate retreat
of shock
C 2003134 Smith and Rosen
period increasing shock velocity
a cooling layer increas­
ing temperature shock layer inflates. higher temperature
gas efficiently,
so expands outwards, driving
shock
to even higher speeds. The shock perform
ing periodic oscillations
to high­amplitude quasi­periodic
oscillations Imamura, Wolff Durisen 1984; Strickland
&
Blondin 1995).
This dynamical overstability relevant atomic shocks.
For velocities greater than approximately
s
, strong atomic
cooling near
K predicted initiate overstability
Innes, Giddings
& 1987; Gaetz, Edgar Chevalier 1988).
instability expected, however, molecular shocks
all
individual cooling functions positive functions
of tempera­
ture, shown Fig.
1.
therefore surprised
to find
a dynamical instability
in
numerical simulations. version
a three­dimensional
(3D) molecular hydrodynamic code, weperformedone­dimensional
high­resolution code improves Suttner
et
al.
(1997) including
H
2 chemistry cooling
Appendices). previous code not signs dynami­
instability. fact,
it production
2 Omolecules
cooling layer
to
of catastrophic cooling.
The formation
of molecules causes instability
by altering
cooling function
is time­dependent
as moves down­
stream. Note (1983) recognized that rapid formation
water molecules could
to thermal instabilities molecular
clouds. The fundamental condition dynamical instability
derived through numerical solutions specific cases, especially
for
monatomic with
a cooling proportional square
of
density, with atomic shocks being considered (Chevalier
& Imamura
1982; Chanmugam, Langer
& Shaviv 1985; Dgani Soker 1994).
Numerical simulations demonstrate instability leads
amplitude oscillations (Strickland Blondin 1995). Section
2,
present hydrodynamic equations discuss criteria
for
which molecular shocks should
be unstable, unstable modes
and oscillation periods. also present cooling
tion temperature steady­state shocks equilibrium
C
and
O chemistry employed
in numerical calculations
for
which time­dependent oxygen carbon chemistry
is practical
Figure cooling line) cooling functions
a
with indicated fixed density
2 fraction. cooling components
are:
1, gas--grain;
2,
2 ro­vibrational;
3, atomic;
4,
H
2 rotational;
2
O
vibrational
2
;
6,
2 vibrational
H dissociation;
8,
H
2
reformation heating (thin double line);
9, rotational; vibrational
with
H
; COvibrational with H;12,
I] structure; rotational.
to follow. compare cooling functions steady­
shocks non­equilibrium chemistry conclude
instability will probably present, simulations
time­dependent chemistry ionization needed verify this
conjecture determine minimum abundances necessary. Note
molecular hydrogen formation rate remains
certain influence instability regime.
In Section present numerical results demonstrate
nature instability yield time­scales. strong influ­
instability observational properties, emphasis
infrared emission
is then discussed (Section 4.1).
Several conditions must fulfilled instability oper­
First, sufficient
2 gas below
#8000 Hollenbach McKee (1989) have shown and
water cooling indeed dominate steady shocks
in high­density
regions. consider carbon oxygen depletions that
stabilize shock. Ameliorating factors also include thermal con­
duction and depletion. Furthermore, sufficiently strong trans­
magnetic cushions shock inhibits dynamical
instability (T’oth Draine 1993).
that radiation will influence location which
molecules reform (Hollenbach McKee 1989). shocks with
speed
v
#
80
s
, radiation will inhibit formation
of
2 molecules. However, calculation Neufeld
Dalgarno (1989a) shows cooling then replaces cooling,
providing
a similar steep increase
in cooling temperature falls.
Nevertheless, restrict results shocks below speed.
shock instabilities operating. These mainly
distort fragment dense gas accumulates either
between shocks
or between
a single shock medium
high­thermal pressure. The instabilities include Vishniac thin­
linear instability (Vishniac 1983; Mac Low
& Norman 1993),
non­linear thin­shell instability (Vishniac 1994) trans­
acceleration instability (Dgani, Buren
& Noriega­Crespo
1996). They require two­dimensional studies, whereas
cillatory instability present even one­dimensional simulations.
Moreover, instability occurs
in cooling layer rather than
accumulated cold possessing distinct signatures.
2 THE FLOW EQUAT
I ONS AND
ANALYT
I CAL PRED
I ON
The equations
model basic radiative shocks: hydrodynamic flows
dimension. Numerically, will solve time­dependent flow
equations: (#v)
x
(#v)
t (#v
2
0
#(ev) -#(T,
f
)
f
f nv)
n,
f
)
- D(T,
),
where
n
is hydrogen nuclei density,
e internal energy
density and
f
is molecular hydrogen abundance
2
f take
a helium abundance n(He)
= Therefore, total
particle density
n (1.1
-
f and temperature
C
# 2003 RAS, MNRAS 133--147

Shocks molecular clouds 135
p/(kn
tot internal energy loss through radiation chemistry
per volume
is dissociation reformation rates
of
molecular hydrogen
by
D and
R, respectively.
The internal energy related pressure p/(#
1)
where effective specific ratio taken
as
5.
5
f
3.
3
-
f
(5)
which assumes
2 possesses two rotational degrees
of freedom.
We define
#
tot
p
=
1.
4
1. f.
(6)
2.2 Power­law cooling predictions
Shock stability studies have been mainly restricted flows
power­law cooling
of
0
#
#
(7)
(Chevalier Imamura 1982; Wolff, Gardner Wood
Strickland Blondin 1995) superpositions power­law
ing Saxton
et 1998 further references). Thermal instabil­
ities interstellar medium first analysed detail
(1965). Smith (1989) investigated small­wavelength isobaric
condensation modes within shock waves. The instability condition
(equation
of Smith 1989) reduces
to
# 1)/(3
-
#
)
shock and approaches condition
<
# downstream.
Note that
#
=
2, instability condition remains close
<
whereas
1 condition more stringent.
The oscillatory instability condition
is even more stringent,
pending upon cooling processes boundary conditions.
For
# oscillations unstable fundamental
mode
#
0.
4 and
in higher modes
c where values
c exceed (Chevalier Imamura 1982). Note, however,
apparently
a range values
#
< which non­linear
effects critical (Strickland Blondin 1995).
In this range,
ondary shocks form near wall restrict, not
collapse. fundamental
is thought most signifi­
cant since long­wavelength perturbations would least damped
motions transverse direction. Numerical simulations
demonstrate, however, higher­order modes generate
lations
#
0. (Strickland
& Blondin 1995).
applications
to shock propagation through molecular
we
should consider
#
=
7
5 (fully molecular)
# (fully molec­
ular cent helium). Furthermore, unity several
cooling processes high­density CO rotational
2
ra­
diative cooling. dissociative shocks, however,
hydrogen almost completely atomic during cooling
down below that instability
is caused
2 vibrational cooling (with molecules excited collisions
with which
# these respects, existing analyses
provide relevant criteria present conditions.
2.3 steady­state model
The cooling function
a uniform constant density
and chemistry displayed Fig.
1. component cooling
tions listed Appendix
A chemical reactions Appendix
They were selected problems associated shocks
in
dense molecular clouds. assumes equilibrium chemistry
described Appendix conclude cooling functions
Figure cooling functions corresponding
to dissociative shock waves
prove
to unstable, assuming equilibrium chemistry
for
C with
standard conditions
as shown abundances
(
5
â
=
-4 panel)
#
=
â
2
(lower panel) Alfv’en speed
s
. panel corresponds
to unstable depleted
O case. individual cooling/heating components
labelled
positive functions temperature with power­law index
being above oscillatory instability critical value each
section. This holds only indicated parameters under
much wider conditions.
unstable shocks that motivated this analysis involve com­
multicomponent time­dependent cooling function Im­
mediately following shock front,
is zone
of steep cooling,
dominated
2 radiative and dissociative cooling atomic cool­
Therefore, expect with
T
> 8000
be stable.
that thin section maintained throughout, moving
accommodate oscillating cool layer, signals trav­
elling quickly through signal propagates slowly through
following cool layer, which refer `infrared layer', and
determines dynamical instability properties.
steady­state model solves hydrodynamic equations
post­shock (with partial derivatives being dropped).
shock
is replaced discontinuity satisfying
Rankine--Hugoniot conditions, attention being spe­
heat ratio. pre­shock parameters represented
p
0
,
v
,
0
,
#
0
,
#
0
,
f
0
2
=
(#
v
2
/#
p
0 Since molecular
dissociation requires
a finite time, immediate post­shock values
f
f
,
1
0
,
C 2003 MNRAS 339, 133--147

136 Smith and Rosen
S
1
0
0
1
0 1)M
2
0 1)M
2
0
(8)
p
1
0
=
1
1
M
(9)
and
1
0
=
p
0.
The radiative
on eliminating
e, described
=
1
1
p
2
1
1
v
and, substituting,
#[p/(#
-
1)]
+
#
p
= -#(n,
f
),
which, using equation substitute from
##
#x
-4. (3.
3
)
2
#
and
#
f
1.
p
#
1
1
leave single first­order equation that solved numerically
without difficulty.
Before running steady shock model simulations, added
equilibrium oxygen carbon chemistry addition hydro­
gen chemistry. basic reactions form
OH molecule were chosen (Appendix Previously, had
abundances follow molecular hydrogen fraction
not included
2 simplest approximation that allowed
us
model molecular
in dimensions (Suttner 1997).
Here, however, assume conversion oxygen
and and
O described equilibrium
dances through reactions
2
C. This especially critical
large­scale simulations since advantageous
to avoid
introduction
of variables.
H
2
form ahead
2 reformation since small fraction
2 sufficient instigate complete conversion
oxygen. Note OH oxygen structure cooling included
steady­state slab models, excluded from
lations. find cooling contributions neglected
present shock simulations densities exceeding
3
4
,
as demonstrated hydrogen nucleon density
6
-3
Fig.
2.
Time­dependent chemistry included steady­state
calculations
to test equilibrium assumption. Equilibrium chem­
istry found
to provide reasonable estimate cooling function
high densities. demonstrated
by comparing panels
of Figs
2
and
3, cooling generated
H
O formation
is weaker
equilibrium This would influence oscillatory structure
do believe that sufficient remove instability.
instability present when either
C
or depleted
an
order magnitude, yielding more gradual increase
in
ing (e.g. lower
of Fig.
2 where
2 bump
is
Figure cooling functions dissociative shock waves non­
equilibrium chemistry standard conditions described
individual cooling/heating components labelled
greatly reduced). We note predict detailed structure,
dynamical processes
as critical chemical effects. Two­
dimensional simulations full time­dependent chemistry will
needed.
a density
4
, however, non­equilibrium chemistry
is certainly required
to reproduce important contributions
oxygen fine structure water cooling below 8000 shown
in non­equilibrium chemistry, increase
in
vibrational cooling cools below
K
is longer
mirrored
in cooling.
Dynamical instabilities may occur when molecular hydrogen
fraction inversely related temperature, such possible
in cooling layer shock,
as now discussed. cooling
is presented volume,
#,
in since cooling
functions simple power­law functions
of density. Hence,
to predict which these states would unstable, need inter­
critical temperature indices discussed Section terms
of index
#
by
#
f
(T
#
T
#
T
, which follows
a fluid parcel through cooling
at almost constant pressure.
interpretation contains uncertainty, however, since chem­
(and, hence, cooling)
is time dependent. Nevertheless,
that
a high negative value
#
is present temper­
regime between 3000 8000 standard non­depleted
O abundances. Moreover, width shock
is mainly
determined
by temperature which rapid cooling about
to
K confirmed Fig. This loca­
adjacent cooling zone possesses
a inverse
C
# 2003 RAS, MNRAS 133--147

Shocks molecular clouds 137
Figure comparison
of equilibrium non­equilibrium chemistry
cooling functions within dissociative shock waves density
of
4
standard conditions described individual
cooling/heating components labelled
as
dependence temperature.
us
to conclude
shock
is,
at least, linearly unstable.
example
C strongly depleted confirms
cause instability. Low oxygen and carbon abundances permit
little
2 CO OH form. shown
in
should stable since
#
# -1.39 except narrow temperature
range,
as confirmed simulations (lower panel
of
Section 3.6).
HYDROCODE
I MULAT
I ONS
3.1 numerical model: ZEUS­3D
We employ ZEUS­3D one­dimensional mode (Stone
&
Norman 1992) provide basic hydrodynamics. This
is second­
order Eulerian finite­difference code. Here study compressible
hydrodynamics without gravity, self­gravity thermal conduction.
physical viscosity
is modelled, numerical viscosity remains
present,
a Neumann artificial viscosity determines
sipation shock front.
module molecular chemistry molecular atomic
cooling added. functions rates listed
in
Appendices. ultimate goal develop
a reliable code
which tackle three­dimensional molecular dynamics,
adding gravity, magnetic fields, ambipolar diffusion radiation.
Figure Profiles temperature, density, pressure molecular fractions
(H
,
H demonstrating collapse shocked region.
display profiles from simulation covering stages
of cycle.
shows collapse
in temperatures
of #500--
initiated
by pressures (seen
as
the
in density
collapse) region.
restrict cooling chemistry
to just those items
essential
to dynamics. employed simultaneous
method discussed Suttner (1997)
in which time­
adjusted
in
to limit change
in internal energy
in any
present primarily two­shock simulations. addition, have
examined single shocks, reflection boundary condition simu­
lating
a wall. Symmetric inflow boundary conditions chosen
ends grid
in two­shock simulations. Hence, dense
cold material accumulates centre grid. then
behaviour each shock independent accu­
mulated material behaves exactly wall. Initially, however,
is effective wall shock.
a direct result, twin shock
pattern might display oscillatory instability since waves
brought resonance through reflection. Thermal instabil­
present, however, and rapidly generate some structure and
asymmetry.
C 2003 MNRAS 339, 133--147

138 Smith and Rosen
Figure cooling functions within
a stable dissociative shock
with abundances
#
2
â
#
=
10
ditions
as individual cooling/heating components labelled
3.2 Resolution convergence
Once central dense shocks usually
begin behave independently dynamical instability
in.
Table
1 includes data detailed resolution study
exploration parameter space standard two­shock case,
shock--wall example. define shock resolution
as
shock cooling length zone length.
Fig.
7 displays locations
of shock defined
positions
of maximum temperature) one shock
a
function time. This demonstrates presence
of instability
at
resolutions. however,
at resolutions above
4
10
variation time­scale amplitude approach
a consistent
behaviour. remark that oscillations change period
at
irregular intervals. propose these changes, encountered
Table
1. Shock parameters.
a logn
v
s 1­/2­shock unstable Period
b
8(10)
6
60
2(10)
6
60
3.2(11)
6
40
1.6(11)
6
40
8(10)
6
40
6D40­C 8(10)
6
40
t
< 500
S: 8.7,
<
t
< 1000 35.6, 9.4--16.6,
1000
<
t
yr:
M: 49.8,
<
t
< 750
S:
6D40­O 8(10)
6
40
M: 0.27,
6E40­CO 8(10)
6
40
4(10)
6
40
2(10)
6
40
6F40­1 2(10)
6
40 0.42, 1.3,
8(10)
6
20
1(10)
6
20
3.2(11)
5
40
1.6(11)
5
40
2.56(12)
4
40
1.28(12)
4
40
a designation simulation order, number density,
a representation resolution, pre­shock
velocity
-1 letters resolution
are associated number zones
in cooling length,
as
follows: 8--16; 16--32; 32--64; 64--128;
E, 128--256; and
F, 256--512.
b period years determined using indicating
a weak, moderate strong peak, respectively.
in atomic shocks, caused nature cooling, which
possesses
a time dependence through molecular abundances.
typically require shock resolution
of #180
to obtain physical
instability.
It significant the instability generates some chaotic
haviour, different cooling components able
to produce dif­
ferent resonances. thus find, similarly studies
of turbulent
that behaviour converges structure,
as
resolution increased.
notable points
as follows.
(i) When
a low resolution employed, shock
is unstable
unphysical long­term high­amplitude variations found.
increased period
at which shock width recorded
100 Along with change, more short­period variations
become apparent. Hence, spiky structure that appears during
shock evolution curves caused display
resolution rather than numerical resolution.
behaviour does converge
a single pattern
passes through various cyclic patterns. may caused
by
presence
of different harmonic modes interaction through
complex cooling function.
The collapse re­expansion
display excerpt
in from recorded
a time
resolution displayed short variations
time­scale superimposed variations
a time­scale
of Although most collapses expansions exhibit smooth
haviour, collapses place high speed. This
is displayed
explicitly
9 where plot distribution shock
velocity,
as calculated between each consecutive pair.
asymmetry demonstrates that collapses often,
C
# 2003 RAS, MNRAS 133--147

Shocks molecular clouds 139
Figure dependence instability spatial resolution. position
of maximum temperature clarity,
of shock)
a function
time simulations (from down) 6C40, 6D40, 6E40 6F40 span different spatial resolutions.
,
is displayed
in each
panel. The frequency
of displayed identical each panel, show datum
for
the thereafter one
yr.
Figure Expanded view short­term oscillations
in high­resolution
case. position maximum temperature
a function
shocks, where
is shown
a dotted high­resolution
simulation 6F40. positions plotted shows
the
similar,
if somewhat time, behaviour independent shock.
exclusively,
at high speed,
in range
s
probable shock front velocity
is larger average
shock speed
-1
. Note these speeds
an
order magnitude smaller shock speed, they
large­amplitude variations shock length.
The distances shock fronts the centre
displayed Since investigating instability,
reason shock location should coincide.
Figure probability distribution shock velocities high­
resolution relative fraction velocities, which deter­
mined consecutive times positions, shocks
high­resolution simulation 6F40, portion
of simulation
0.1­yr intervals. The
of velocity
-1
Shock profiles
displays profiles temperature, density, pressure and
molecular fractions
at different phases during collapse. The
upper panel demonstrates collapse indeed coincides with
shrinking warm infrared layer. This layer determines
C 2003 MNRAS 339, 133--147

140 Smith and Rosen
length
of shock cools efficiently temperatures
below until some COmolecules formed. end
of
collapse, remaining infrared layer consists #8000
K
gas. other hand, both occupy
regions. Note pressure trough appears when shock wide.
low pressure results shock collapse.
In this produces
a
weak secondary shock returns upstream.
The logarithmic profile hydrogen molecular fraction
tains simple shape since
2 reforms infrared
Fast reformation does occur cooled below
300 which close
to interface between shocks.
shape
of
H
2
O fraction follows
2 quite closely.
contrast, fully reforms infrared layer. Note generation
secondary shock
x =-3â
12
in figure) during
collapse that
is sufficiently strong
to destroy most reforming
COmolecules. Although shock
in comparison
mary shock,
it occur regime where
2 and chemistry
most sensitive temperature.
3.5 Density velocity
Fig. shows that instability
is present densities consid­
ered here. code assumes equilibrium chemistry remains
a
valid approximation other cooling mechanism
neglected. density
of
10
, however, assumption
of
equilibrium cooling invalid slow formation CO
bined enhanced
O structure cooling, leads
a
flatter cooling function.
standard
-1 that highest ampli­
tude greatest variation time­scale occur lowest density.
The time­scale dependence expected since lower densities
have longer cooling times. apparent from difference
in
shock width shown Fig. Quasi­periodic oscillations
shock found some stages simulations.
Figure The dependence dynamical instability
on density. position maximum temperature
as function time different densities
plotted simulations 6E40, 5F40, 4F40
the down). Fig.
7 displays two­shock simulations, show position
of only
shocks. Here, interval
yr
in 6E40, every
in every 0.25
yr
for
all
t
in simulation 5F40, every
4F40.
the shock
is inversely proportional density, vertical different. time­scale oscillations
increases density
is decreased.
employed Fourier transform (FFT) highest­resolution
and present
in
1 some periods that appear.
oscillatory pattern changes periods that show
in analysis very weak (<3#
, denoted
W
in
Table). Still, these values show inverse relationship between
period the pre­shock density. Specifically, periods
approximately
6 where
6 n/10
6
, range
densities simulations shown
in Fig.
A study dependence shock width
on shock speed
is presented density
6
. The
instability present speeds above km
s
, where
drogen molecules destroyed immediately behind shock
rapid reformation leads
to dynamical instability. have
investigated shocks
-1 since,
at speeds,
radiation from gas delays formation
of molecules
infrared layer. Nevertheless, instability could
deeper infrared layer. includes radiative transfer
effects
be necessary confirm
that altered time­step display
at 100
In two­shock simulations, allows shocks progress
point where
is independent one another and enables
us
to display structure that would recorded through discrete
observations made every years.
manner which shock structure changes depends
strongly
on upstream speed imposed.
s example
occasionally displays narrow shock regions, roughly 30­yr
interval. Cyclic patterns uncovered
s
.
At
20
s shock stable.
Chemistry oscillatory instability found
in previous simulations
of shocks (Suttner
et because, owing lack
computing power, oxygen and carbon chemistry included.
these prove vital nature
of
a fast shock. This
C
# 2003 RAS, MNRAS 133--147

Shocks molecular clouds 141
Figure The dependence instability pre­shock velocity. position maximum temperature
as function time
3 different initial
(pre­shock) velocities simulations 6D60, 6E40, 6D20 (from down). with previous position figures two­shock simulations, show
position one
of
the shocks. Here
are shown year
t
<
yr
yr
in
every simulation 6D20. Note different vertical scales,
the
40
s simulation largest region.
confirmed which displays evolution
width different
O abundance combinations.
that shock stable only when
C depleted.
2 formation generate instability,
as expected
shape
of cooling function
in steady shock analysis.
evident when molecule responsible
instability (middle panels),
a regular oscillation
is generated.
suggests that changes pattern other simulations mainly
caused multipeaked cooling function. small­amplitude
oscillations uncovered oxygen­depleted shock related
to
inefficiency
of
to continue cool below #2000K.
Note when instability
is instigated through molecule
formation cooling, high­amplitude oscillation
is apparent.
In­
teresting changes temporal behaviour found carbon­
depleted shock simulation shown detail
A dramatic
change occurs
t
= seems plausible shocks
prior
t
= interact through the dense, wall,
suspect change period, effectively
of previous
one, signifies become enough allow
shocks behave independently.
3.7 Single shock impacting
We also verified
a model produces same unstable
structures
in equivalent two­shock model, shown
Hence, conclusions independent chosen shock
set
up, given limitations one­dimensional shock waves.
dimensions, that forms between shocks well
and fragment owing induced transverse motions.
MPL
I CAT ONS
4.1 Observational aspects
We calculate the effect
of instability observable quanti­
ties. First, chose molecular hydrogen emission
to represent
signatures from infrared layer. The upper energy level asso­
ciated
H
2
1 transition corresponds
to
an excitation
temperature
of nearly 7000 Emission arising from cascade fol­
lowing molecule reformation accounted here assume
energy released reformation efficiently redistributed
thermal component collisions).
displays
in grey­scale showing location
of emission
as
a function integrated emission
is shown
in panels
right. also show
2
1 emission, which arises
more excited
2 spatial distribution, however,
is very similar
to and
is illustrated.
panels demonstrate characteristics
a steady shock,
produced taking abundances. central wall
slowly builds
in with constant velocity, producing
expected triangular shape.
2 emission
is spread large
region, since shock dominated spatially
by slowly cooling
between 4000 8000 The integrated emission
of
H
constant, with diagnostic
0 (1)]/F[2
1 This consistent with range values from
= standard steady case
=
25 reduced abun­
dances example
= 1.0(-4) and
= 2.3(-4)] tabulated
Hollenbach
& McKee (1989).
The bottom row
of shows example large variations
integrated flux
of
H
a simulation where instability
is active. Thevariations 1#0S(1) order magnitude.
caused
by rapid appearances disappearances
of
whole 2000--4000
K section. This causes
R vary
in
range #1.5--30. time­scale variation than
some significant variations
in time­scale
of months,
in this
high­density example. from
a low­lying rotational (dashed line)
overwhelmingly generated
in downstream wall,
is simply related accumulated mass (Fig.
Variations atomic flux displayed) present
relatively small since dissociated cooling layer
is preserved throughout. shock front speed does vary, however,
C 2003 MNRAS 339, 133--147

142 Smith and Rosen
Figure dependence
of instability abundances carbon oxygen. position maximum temperature
a function four
combinations
of nominal depleted
C and
O abundances
is plotted. panel, simulation 6E40
is presented, default abundances
of #(C)
2.0(-4)
= 5.0(-4). second panel bottom), simulation 6D40­O
is displayed, which abundances
of #(C)
= 2.0(-4) #(O)
2.1(-4) thus almost
of oxygen end
up
in
2
In panel simulation 6D40­C abundances 1.0(-5)
#(O) 5.0(-4). final
is simulation 6E40­CO with abundances #(C)
= 1.0(-5) #(O)
= 2.0(-5). Note simulations
in lower three
panels with similar sizes,
the decreasing cooling larger cooling length­scales change resolution case.
The panels display every bottom panels every
Figure Change period amplitude instability carbon­depleted shock simulation. display position
of shock
in simulation
6D40­C times, intervals between data. dramatic change
in
the character
of instability
at
= prior
is
series paired oscillations, shallow followed
a deeper During period
<
< other shock
is completely
phase with shown, with shallow collapse.
and generates moderate changes more recognizable
signature
of instability
is pattern variation
radial velocity when viewed shock normal.
The radial velocity through sporadic variations
, accompanied
by dips flux.
grated
H emission above, large variations occur
a
months
a density
6
-3 30--50
yr
at
a density
of
.
Physical interpretation
complexity time dependence total molecular cooling
function generate complex unstable patterns. Cooling (atomic and
H dissociative) following shock front extremely rapid, and
treated the shock front itself. These processes
produce
a narrow
in temperature profile (see
thus carried out
a linear analysis modelling shock front
C
# 2003 RAS, MNRAS 133--147

Shocks molecular clouds 143
Figure Variation shock position
a single­shock simulation. display position
of maximum temperature simulation 6F40­1, which
shock reflecting boundary
at
a every thereafter. instability
is present
and variation
of shock position
a amplitude similar two­shock simulation 6F40.
Figure Spatial temporal variations molecular emission lines. display
the results simulations, every
2
yr simulation
6E40­CO and
yr 6F40 bottom
In each displays spatial variation emission
2
#
0 S(1),
second variation
0
0 emission. two shading indicates emission, most
of emission
2
cooling closest
to central Similarly, bulk
of emission, contributions include
a maximum cut­off
at
-1
from central
In rightmost panel row, temporal variation integrated emission
the mentioned
(H
dashed)
2
2
1 (short dashed).
and rapidly cooling
by conditions with
a range
of
effective specific heat ratios shock
of 10--40. Hence,
modelled front treats shock heating, atomic cooling
dissociative cooling discontinuity. subsequent cooling
layer
5
3
. linear instability analysis solves equations
presented Chevalier Imamura (1982) with revised
shock boundary conditions.
that rapid cooling stabilizing
influence shock (Smith Rosen 2002),
as expected
previous linear analysis included
a second steep component
(Saxton 1998). indeed that negative values
power­law index required produce linear instability (Smith
Rosen 2002). example, with
a compression factor
S
corresponding 21/19, first overtone grows
#
-0.
4 and fundamental mode only We expect
find significant regimes which higher harmonics dominate.
The many different manifestations instability have been
explored classified
by Walder Folini (1996). They propose
different oscillation types. particular, first
overtone exists
a strong form results
a high­amplitude
quasi­periodic collapse re­expansion
pearance
of subshocks (M­type). also uncover
a smooth
sinusoidal mode (S­type). Some similarity
of long­term simula­
their simulations suggests overtone may
responsible much behaviour recorded simula­
however, generally demonstrate complex behaviour sev­
superimposed non­stationary usually being present.
The simpleoscillations encountered oxygen­depleted shock
simulation have been examined harmonic signatures. Here,
a simple wave form present reveals harmonic mode.
When shock
at
its maximum width, high­amplitude velocity
grows strongly wave represents infalling
C 2003 MNRAS 339, 133--147

144 Smith and Rosen
material and
is accompanied
a strong mass influx
to
This strong rarefaction wave propagates outwards, until
it meets
collapsing shock front, which collapse
is halted.
shock re­expands, accompanied much smoother veloc­
profile. This sequence repeats throughout simulation.
behaviour
a characteristic overtone, caused
velocity (Chevalier Imamura 1982) and inflow (Walder
&
Folini 1996) being phase with shock front position.
The oscillating modes have described standing waves,
with zero velocity
a reflecting dense cold layer)
shaken shock front (Chevalier
& Imamura 1982). The
pattern analogous
in half­open organ pipe (Saxton
et
al.
1998). yields
a series resonances with angular frequen­
cies
n
= 2#(2n
+
0
,
where
t
0 period
of fundamental mode. Calculated periods
#20.6--22.6 #6.7--7.1 #4.2, respectively, units
of
0
, derived from linear analysis (Chevalier Imamura
# range
1
2 and shocks with
a compression factor
of
S derive relationship
t
0 2#(S
1)/1.
6 from
a
analysis (Smith
& Rosen 2002); (linear) instability period
is
approximately increased compression factor.
Now wish determine whether numerically derived
re­
sults consistent theoretical expectations.
A very detailed
comparison attempted because total cooling function
complex. The cooling beginning
of unstable
section #6000 estimated
ir
#
2
0 (v/40
s
)
2
erg
-3
s Fig. density location
n
#
n
0
m
p
2
0 /(kT) 32.2
n
0
s
). Hence cooling time­scale
6000
n kT/L
ir 120
4
0
) agrees
in
order magnitude variation time­scales
in simulations.
The oscillation dynamical time­scale
is defined
t
d
L
c
ir
,
where
v
ir
is speed 6000
c
â
10
4
0 measured from simulations. yields
t
=
200
4
-3
0 Hence, relevant cooling, dynamical
measured variation time­scales comparable.
The theoretical oscillation period described above
for
harmonic,
c
0 190
4
0
yr
on substitu­
tion
c
. Hence linear analysis consistent
numerical results.
SUMMARY AND CONCLUS
I ONS
We shown cooling layers hydrodynamic molecular
shocks speeds between #30 #70
s subject
to
cyclic, potentially large­scale collapse reformation.
ing
to
a form catastrophic cooling produces overstability
caused rapid formation
of CO
2 molecules. When
C
and
O depleted, instability absent.
Quasi­periodic oscillations dominate sections
simulations. time­scale variations may change.
riodic oscillations present when either
or
H
2 dominate
cooling, controlled abundances.
Large variations predicted those lines produced
in
gas temperatures range 1000--8000 shown
this explicitly near­infrared
H
it also applies
to high­J
CO lines. Quite large radial velocity variations expected
emission produced
in moving with shock front, such
as
those associated with recombination common tracers
of
hydrodynamic shocks optical range, S[II].
application
of results limited
by following assump­
magnetic field must either parallel
to the shock
direction
or have
a transverse component Alfv’en speed
-1 cooling function remain unchanged above
K. Shock flows transverse magnetic fields characterized
bypre­shock Alfv’en speeds
of
s (shock speed 30km
-1
to
s (shock speed
-1 maintain gradient
in cooling initiated instability oxygen­depleted
simulation. suspect that these flows also exhibit
instability. This needs verified with magnetohydrodynamic
(MHD) simulations. Secondly, one­dimensional analysis does
allow other instabilities that require transverse motions
be present. These instabilities depend specific problem and
could make signatures complicated
to interpret.
In this
respect, specific signatures results should allow
determine whether observed variations caused instability.
Thirdly, time­dependent abundances
2 formation dissoci­
included equilibrium
2
is assumed
in
simulations. have shown here cooling function derived
non­equilibrium chemistry
is smoother than, shape
approximately similar that assuming equilibrium pre­shock
densities
n excess
of
4
.
Low­resolution numerical simulations should treated cau­
When sharp gradients abundances smeared
parently stable shock fronts may result. This
in addition
to
coarse mesh problems noted Walder
& Folini (1996), which
eliminate particular modes.
instability occurs
in high­density low­speed dissociative
shocks. For shocks with speed
0
, radiation
inhibit the reformation
of
2 molecules (Hollenbach
&McKee 1989). Hence, restrict our results shocks below
speed, although note unstable conditions arise
through OH cooling
v
#
80
s (see Neufeld
Dalgarno 1989a). variations shock­generated emission lines well
documented. Recently, Chrysostomou (2000) found variability
of
1
#
0 S(1)
2 years several locations within
Herbig--Haro (HH) objects. Herbig Jones (1981) documented
of optical structural variations
1. future, flux
variations
be increasingly documented higher spatial reso­
lution permit observations coherent shock sections. What
would make
it possible associate these variations radiative
cooling instability? series simultaneous spectroscopic obser­
vations emission lines optical infrared would provide
a comprehensive
In the absence
of information,
shown interpretation shocks
as steady flows
lead misleading results.
ACKNOWLEDGMENTS
numerical calculations were FORGE super­
computer, acquired through PPARC JREI initiative SGI
participation.
is funded PPARC. thank referee
many suggestions helped
us improve paper.
REFERENCES
Biham Furman Pirronello Vidali 2001, ApJ,
Chanmugam Langer Shaviv 1985, 299,
Chevalier Imamura J.N.,
Chrysostomou Hobson Davis Smith Berndsen 2000,
MNRAS,
C
# 2003 RAS, MNRAS 133--147

Shocks molecular clouds 145
Dgani Soker 1994, ApJ,
Dgani van Buren Noriega­Crespo
Draine Roberge Dalgarno ApJ,
Field G.B., 1965, 142,
Gaetz T.J., Edgar Chevalier R.A., 1988, 329,
Herbig Jones B.F., 1981,
Hollenbach McKee C.F., 1979, ApJS,
Hollenbach McKee C.F., 1980, 241,
Hollenbach McKee C.F., 1989, 342,
Imamura J.N., Wolff Durisen R.H.,
Innes D.E., Giddings Falle S.A.E.G., 1987, MNRAS,
Jacquet Staemmler Smith M.D., Flower 1992,
B,
25,
Kwan 1977,
Langer Chanmugam Shaviv
Bourlot Pineau For“ets Flower D.R., 1999, MNRAS,
Lepp Shull ApJ,
Mac M.­M., Norman M.L., 1993,
McCray Kafatos R.F., 1975, 196,
McKee Storey J.W.V., Watson D.M., Green 1982, 259,
Neufeld D.A., Dalgarno 1989a,
Neufeld D.A., Dalgarno 1989b,
Neufeld D.A., Kaufmann 1993,
Pineau For“ets Flower D.R., Dalgarno 1988, MNRAS,
Saxton Pongracic Shaviv 1998, MNRAS,
Shapiro Kang 1987, ApJ,
Silk
J., MNRAS, 205,
Smith M.D., 1989, MNRAS,
Smith M.D., 1993, ApJ,
Smith M.D., 1994, MNRAS,
Smith M.D., Rosen 2002, Ap&SS, submitted
Stone J.M., Norman 1992, ApJS,
Strickland Blondin J.M., 1995, 449,
Sutherland Dopita ApJS, 253
Suttner Smith M.D., Zinnecker
T’oth Draine B.T., 1993, 413,
Vishniac 1983,
Vishniac 1994,
Walder Folini
Whitworth Clarke 1997, MNRAS,
Wolff M.T., Gardner
J., K.S., 1989, 346,
Wolfire M.G., K˜onigl 1991, 383,
APPEND
X
: THE COOL
I FUNCT
Here describe cooling function used numerical
(equation
3) present version steady model outlined
Section The function
is composed
of
11 separate parts (one
of
which further components were introduced
steady shock models.
11
#
=1
#
i.
The components follows:
#
is gas--grain (dust) cooling,
#
is
collisional cooling associated with vibrational rotational modes
molecular hydrogen,
#
is collisional cooling atoms,
#
is
cooling through rotational modes water induced collisions
with both atomic and molecular hydrogen,
5 cooling through
vibrational modes
of water induced collisions with molecular
hydrogen
6 vibrational modes water induced
sions with atomic hydrogen,
7 cooling dissociation
of
molecular hydrogen,
#
is heating resulting from reformation
molecular hydrogen,
9
is cooling through rotational modes
of
carbon monoxide induced collisions with atomic molec­
ular hydrogen, cooling through vibrational modes carbon
monoxide induced collisions with molecular hydrogen
is cooling through vibrational modes carbon monoxide induced
collisions with atomic hydrogen. Therefore, excluded
cooling
13 structure cooling
12
present simulations. form dust gas--grain) cooling
is taken from equa­
(2.15) Hollenbach McKee (1989):
#
n
2
#
1 (A2)
where, assuming standard dust properties,
#
1
=
3.
T
T
â
[1. 0. exp(-75/T
)]
s
3. (A3)
In present calculations,
=
20 hence assume
dust cools rapidly after shock front, justified
Whitworth Clarke (1997).
H
2 vibrational and rotational cooling based equations (7)--
Lepp Shull (1983). have shock tested these formula
against detailed subroutines presented Bourlot, Pineau
For“ets Flower (1999) found significant differences
to
detailed shock properties. take
#
n
2
1
vL
+
L
1
rL
)
#
, (A4)
where vibrational cooling coefficients
at and density
first formula corrects
a print error):
L
1.
â
10 exp(-6744/T
s
, (A5)
L
8.
â
10 exp(-6840/T
H
k
+
n
2
k
2
-1 (A6)
terms
H
1) and
k
2
0
1 collisional
excitation rates (-6840/T
) term converts these
to
excitation rates], follows:
k
H
(0,
1)
=
1.
4 exp[(T/125)
- (T/577)
2
T
v
1.
0
T
1/2 exp(-1000/T
T
v
, (A7)
where
T 1635
k
H
=
1.
T
1/2 exp[-28728/(T
1190)]. (A8)
rotational cooling coefficient density
L
dex[-19.
+
0. 474x
1.
]
if
T
r
3. exp(-6118/T
)
if
T
r (A9)
where
T
r 1087 density,
coefficient
is
L
dex[-22. 0. 553x
1.
]
if
T
1. exp(-9243/T
)
if
T
, (A10)
where
T
K,
=
#
H
2
0. 1.
2
H
0. 77# . (A11)
Atomic cooling takes expected form
#
n
2
#
2 (A12)
where have Sutherland Dopita (1993) (with
= -0.5)
of
#
2 where additional thermal
bremsstrahlung term equal
to
is added
to
#
2
T 000
C 2003 MNRAS 339, 133--147

146 Smith and Rosen
Water rotational cooling takes
4
=
(n
2
1.
H )n(H
2 O)#
3
3
1. (T/1000)
#
,
where
#
= 1.35--0.3 log(T /1000) values tabulated
by
Neufeld Kaufmann (1993).
Cooling water vibrational modes
in collision
2 and
H
given Hollenbach McKee (1989):
5
=
1.
n
2
exp(-2325/T
)
exp(-47. 5/T
),
6
=
7.
n
#
2
exp(-2325/T
)
exp(-34. 5/T
). The cooling from dissociation molecular hydrogen follows
Shapiro &Kang (1987) modifications Lepp&Shull (1983)
form:
7
=
7.
#
2
H
k
2
H
H
2
#
where 4.48 dissociation energy,
1.
2 exp(-52400/T
)
â[0. exp(-17950/T
)]
#
3
-1
2
1.
3
â exp(-53300/T
â[0. 0908 exp(-16200/T
)]
3
s
with
1.
n
2
f
1
n
1
n
#
+
1
n
1 ##
-1
and finally, critical densities dissociation collisions
of
molecular hydrogen with atomic hydrogen,
n
1
, itself,
2
,
by following approximations:
1
dex(4. 0. 416x
0.
2
)
2
dex(4. 845
1. 1.
2
-3.
Heating from reformation hydrogen including
a
ing' term, following form:
8
=
7
where
7
R given appendix
H
#)7.
â
-12.
Hence employ parametrize fraction
of released energy
that
is thermalized rather radiated.
Cooling rotational modes through collisions with molecu­
or atomic hydrogen have described McKee (1982).
We our cooling their equations (5.2)--(5.5),
9
=
n
n kT#v
T
n
a
/n
1.
a
cr
1/2
,
where mean speed molecules
v
T
= 8kT/(#m
H
2
cr
=
3.
â
10
T
0.
75
3
-3
3.
0
,
where
T
=
T /1000
K. have average density,
n
=
0.
#
2
).
vibrational cooling from collisions
H
,
(Neufeld
& Kaufmann 1993)
#
1.
â
10
H
2
n
T exp(-3080/T exp(-68/T
). (A28)
vibrational cooling through collisions with atomic
drogen,
#
1.
â
10
H
T
1/2 exp(-3080/T exp[-(2000/T
)
3. ]. (A29)
Oxygen cooling through 63­µm structure line been
included
in steady shock model intend
to include
further simulations
a convenient large­scale dynamical
computations molecular turbulence. omit calculation
upper energy level generates much weaker 146­µm line.
cooling
is taken
#
2.
â
10
O
1
H
10
L
, (A30)
where
A
10 8.95
-5
is spontaneous transition rate,
f
=
0.
6 exp(-228/T
1
0.
6 exp(-228/T
)
0.
2 exp(-326/T
) (A31)
is fractional occupation
of
3
1 level thermodynamic
equilibrium
r
r
r
H
r
H
=
[4.
37
â
0. 0.
6 exp(-228/T
+
1.
T
0.
80
0.
2 exp(-326/T
+
0. 48n
H
) (A32)
r
H
[2.
â
10
T
0.
35
0.
6 exp(-228/T
)
6. 0. 0. exp(-326/T
2 (A33)
collisional (values kindly provided
by Flower),
cluding rates calculated Jacquet
et
al. (1992)
2 collisions
which
a single weighted index
to combine
para values expediency, with ortho ratio
sufficiently accurate
for present purpose.
cooling
is taken (see Hollenbach
& McKee 1989)
#
2.
2
T
3/2. (A34)
APPEND
I
: THE CHEM
I STRY
consider limited network
of chemical reactions, which
these one­dimensional simulations which
is
cumbersome included
in three­dimensional MHD self­gravity
simulations. These reactions
(1)
H
H#
2
(2)
H
+
H#
(3)
H
+
H
2#
2
+
H
2# H# C#
+
H#
+
2#
H
2
2
+
H#
H
.
reactions included
a semi­implicit cooling­
chemistry step calculate temperature
2 fraction.
C
# 2003 RAS, MNRAS 133--147

Shocks molecular clouds 147
Reaction
on surfaces rate taken
Hollenbach McKee (1979):
R
=
(3
s
)
â
f
0. 04(T
g
)
T
with factor
f
a exp(-600/T
g
)]
,
which means simulations,
f
a
is quite close
this formula
is subject ongoing debate (e.g. Biham 2001),
one also
f
a
to test sensitivity.
Dissociation given Appendix
A.
fix oxygen carbon abundances,
#
0 (O)
0
then calculate equilibrium abundances, assuming
that reactions much faster
H
2 reactions.
remaining free oxygen
is distributed
O,
2
cording equilibrium abundances. manner, calculate
approximately chemistry within shocks without overloading
hydrocode.
At temperatures below #200
K, equilibrium
dances unlikely reached within cooling layers. The
error will
2 abundance: cooling below
a
hundred Kelvin, code predict oxygen will
in
O,
whereas
H
O forms would remain
as
considerably longer. gas--grain cooling dominate, however,
between 1000
K shape cooling function
this region approximately correct.
Predicted strengths
O
2
O however,
will error. differences published the
tions (4)--(6) Hollenbach
& McKee 1989; Pineau For“ets,
Flower Dalgarno 1988). From Hollenbach McKee (1989),
two
=
2.
T
1.
93 exp(-3940/T
3
-1
,
where
T
K, simplicity, assume
2
in
ground vibrational state,
# (4b)
=
6. 1. exp(-2970/T
)
3
s
, (B4)
which yield equilibrium
of
#(O) #(OH)
0. -0.
4 exp(970/T
2
). (B5)
In comparison, formula
of Pineau For“ets
et (1988) yields
#(O) #(OH)
0. exp(1030/T
) n(H) n(H
2
). (B6)
abundance
of relative determined,
abundance then given through
# (5a)
1.
1/2
300
3
-1 (B7)
# (5b)
=
1.
1/2 exp(-77700/T
)
s
, (B8)
which combine equilibrium CO abundance
#(CO)
# ##(OH)
+ ##(OH) (B9)
where
#
# 700
#
n(H). (B10)
CO abundance determined, equilibrium
H
O and
abundances given equation (B6)
# (6a)
8. 1.
300 exp(-1429/T
)
3
s
, (B11)
# (6b)
=
7. 1. exp(-9140/T
)
3
s
, (B12)
which combine equilibrium:
#(OH)
2
=
8. -0. -7711
#
2
). (B13)
paper been typeset
a
A
T
X prepared author.
C 2003 MNRAS 339, 133--147