Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://star.arm.ac.uk/preprints/381.ps Äàòà èçìåíåíèÿ: Fri Aug 2 12:50:59 2002 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 21:26:59 2012 Êîäèðîâêà: Ïîèñêîâûå ñëîâà: m 5

Mon. Not. R. Astron. Soc. 000, 000{000 (0000) Printed 1 August 2002 (MN L A T E X style le v2.2)
Hydrodynamical simulations of the decay of high-speed
molecular turbulence. I. Dense molecular regions
Georgi Pavlovski, 1? Michael D. Smith, 1 Mordecai-Mark Mac Low, 2
& Alexander Rosen 1
1 Armagh Observatory, College Hill, Armagh BT61 9DG, Northern Ireland, U.K.
2 American Museum of Natural History, Department of Astrophysics, 79th St. at Central Park West, New York, NY 10024-5192, USA
Accepted .... . Received .... ; in original form ....
ABSTRACT
We present the results from three dimensional hydrodynamical simulations of decay-
ing high-speed turbulence in dense molecular clouds. We compare our results, which
include a detailed cooling function, molecular hydrogen chemistry and a limited C and
O chemistry, to those previously obtained for decaying isothermal turbulence.
After an initial phase of shock formation, power-law decay regimes are uncov-
ered, as in the isothermal case. We nd that the turbulence decays faster than in the
isothermal case because the average Mach number remains higher, due to the radiative
cooling. The total thermal energy, initially raised by the introduction of turbulence,
decays only a little slower than the kinetic energy.
We discover that molecule reformation, as the fast turbulence decays, is several
times faster than that predicted for a non-turbulent medium. This is caused by moder-
ate speed shocks which sweep through a large fraction of the volume, compressing the
gas and dust. Through reformation, the molecular density and molecular column ap-
pear as complex patterns of laments, clumps and some diuse structure. In contrast,
the molecular fraction has a wider distribution of highly distorted clumps and copious
diuse structure, so that density and molecular density are almost identically dis-
tributed during the reformation phase. We conclude that molecules form in swept-up
clumps but eectively mix throughout via subsequent expansions and compressions.
Key words: hydrodynamics { shock waves { turbulence { molecular processes { ISM
{ clouds { ISM: kinematics and dynamics
1 INTRODUCTION
An understanding of turbulence is of fundamental impor-
tance to many research areas in astrophysics. In particular,
compressible turbulence plays a central role in the process
by which stars are formed (reviewed by Vazquez-Semadeni
et al. (2000); Padoan et al. (2001)). Stars form out of dense
clouds of molecular gas, which have formed out of diuse
clouds of atomic gas. The turbulent energy deduced from
observations of the molecular gas is suôcient to delay the
gravitational collapse, making molecular turbulence of great
interest. The purpose of this study is to provide insight
through the rst three-dimensional simulations of molecu-
lar turbulence.
Isothermal or polytropic equations of state have been
the central premise of previous investigations of turbulent
? email: gbp@star.arm.ac.uk; mds@star.arm.ac.uk;
mordecai@amnh.org; rar@star.arm.ac.uk
models for molecular clouds (e.g. Vazquez-Semadeni et al.
(1996); Mac Low et al. (1998); Padoan et al. (1998); Stone
et al. (1998)). This simplication has the advantage that pa-
rameter space can be explored with numerical simulations,
and the in uence of magnetic elds and gravity can be de-
termined in some detail (e.g. Balsara et al. (2001); Passot
et al. (1995); Klessen (2000); Heitsch et al. (2001) ). The
isothermal approximation is indeed valid in many specic
models where molecular cooling is eôcient and low gas tem-
peratures are maintained. Provided the temperature is low,
one can argue that an isothermal equation of state is as
good as any other, with little feedback from the the thermal
pressure on the dynamics.
However, for cloud formation and evolution molecular
chemistry and cooling is critical (Langer et al. 2000; Lim
2001; Lim et al. 1999). Molecular hydrogen forms most ef-
ciently where the gas is warm but the grains are cool (H2
forms mainly when atoms combine after colliding and stick-
ing to dust grains). Simple molecules like OH, CO and H2O
c
0000 RAS

2 Georgi Pavlovski et al.
form in the gas phase with H2 as the reactive agent. These
molecules are not only important coolants, but associated
emission lines provide a means of measuring the cloud prop-
erties. Molecules are dissociated as a consequence of fast
shocks, UV radiation, X-rays and cosmic rays (Herbst 2000).
We thus need to study molecular turbulence to determine
the distribution and abundances of molecular species.
Even in cool optically thin regions, as assumed here,
molecular pressures may in uence the dynamics. Strong
shock waves may be driven from hot atomic regions, formed
by shock waves, into the proximity of cold molecular gas.
Areas devoid of dust, or in which molecule formation is in-
eôcient (e.g. because atoms will evaporate o warm grains
rather than recombine) may attain and maintain pressures
as high as turbulent pressures.
The rate of decay of supersonic turbulence is important
to the theory of molecular clouds. A possible consequence
of the rapid decay of kinetic energy is that the turbulent
clouds we observe have short lives. The short lives, however,
may not provide suôcient time for relaxation into thermo-
dynamic and chemical equilibrium. Hence, cloud structure
and content evolve simultaneously as a cloud evolves dy-
namically, rather than proceeding through a series of equi-
librium states. Therefore, conclusions based on assuming a
molecular fraction to be a function of density alone (e.g.
Ballesteros-Paredes et al. (1999)) might not always be ac-
curate. We note, however, that the conclusion reached by
Ballesteros-Paredes et al. (1999) of rapid cloud formation is
supported by this study.
Decaying supersonic turbulence is the subject of this
initial study. The work is a direct extension, in terms of code
and initial conditions, of isothermal simulations presented
by Mac Low et al. (1998). Our goals here are as follows:
 To directly compare the decay of the kinetic energy
Ekin with those derived in the papers of (Mac Low et al.
1998) and (Smith, Mac Low & Zuev 2000; Smith, Mac Low
& Heitsch 2000).
 To determine the distribution of the molecular column
density as opposed to the total column density (Vazquez-
Semadeni & Garca 2001)
 To study the rate of formation of molecules.
We omit in this study magnetic eld, self-gravity and pho-
todissociating radiation. We begin with a fully molecular
cube of dense gas and apply an initial turbulent eld of ve-
locity perturbations. The turbulence was chosen to be suf-
cient to dissociate the gas in one case (root mean square
(rms) speed of 60 km s 1 ) and to leave the gas molecular in
another case (rms speed of 15 km s 1 ).
2 METHODS
2.1 The equations
Numerically we solve the time-dependent ow equations:
@
@t +r
( v) = 0; (1)
@ ( v)
@t + (v
r)v = 1
 rp; (2)
@e
@t + v
r e = pr
v +  (T; n; f) ; (3)
@ (fn)
@t +r
(fnv) = R (T; n; f) D (T; n; f) ; (4)
where n is the hydrogen nuclei density, e is the internal en-
ergy density and f is the molecular hydrogen abundance (i.e.
nH 2
= fn). We consider the gas as a mixture of atomic and
molecular hydrogen with 10% of helium (i.e. nHe = 0:1 n).
Therefore, the total particle density is n tot = (1:1 f) n,
mass density is  = 1:4nmH and the temperature is T =
p = (k n tot ). An internal energy loss term has been added to
the r.h.s. of the energy ow equation:  is the loss of en-
ergy through radiation and chemistry per unit volume. The
function consists of 11 distinct parts, summarised below. R
and D are reformation and dissociation rates of molecular
hydrogen respectively. A detailed description of the chem-
istry can be found in Smith & Rosen (2002). Formulae for
, R and D are given in the Appendix A
2.2 The numerical model: ZEUS-3D
As a basis, we employ the ZEUS-3D code (Stone & Norman
1992a,b). This is a second-order, Eulerian, nite-dierence
code. We study here compressible hydrodynamics without
gravity, self-gravity or thermal conduction. No physical vis-
cosity is modelled, but numerical viscosity remains present,
and a von Neumann articial viscosity determines the dissi-
pation in the shock front.
Further coding for molecular chemistry and molecular
and atomic cooling has been added. Our ultimate goal is
to develop a reliable code with which we can tackle three
dimensional molecular dynamics, later adding self-gravity,
magnetic elds, ambipolar diusion and radiation. We thus
restrict the cooling and chemistry lists to just those items
essential to the dynamics. We have employed the simulta-
neous implicit method discussed by Suttner et al. (1997) in
which the time step is adjusted so as to limit the change in
internal energy in any zone to 30%.
The cooling is appropriate for dense cloud material
of any atomic-molecular hydrogen mixture. We include H2
ro-vibrational and dissociative cooling, CO and H2O ro-
vibrational cooling, gas-grain, thermal bremsstrahlung and
a steady-state approximation to atomic cooling (see Smith
& Rosen 2002, Appendix A) and Appendix A of the article.
We take a very basic network of chemical reactions.
Time-dependent hydrogen chemistry is included, but C and
O chemistry is limited to the reactions with H and H2 which
generate OH, CO and H2O (see Smith & Rosen 2002, Ap-
pendix B) and Appendix A of the article. Equilibrium abun-
dances are calculated, which are generally accurate within
the shocks where molecules are rapidly formed and de-
stroyed. In cold molecular gas, however, the available oxygen
will probably remain primarily in the form of water even if
c
0000 RAS, MNRAS 000, 000{000

Simulations of molecular turbulence 3
Table 1. Power law indices  of kinetic energy decay. The 
coeôcients are derived from linear tting on the time interval
[90; 600] years. Average error 0:006.
resolution initial rms velocity for the run
15 km s 1 30 km s 1 60 km s 1
32 3 1:19 1:33 1:33
64 3 1:32 1:38 1:36
128 3 1:45 1:48 1:52
256 3 1:34 ? 1:51 1:60
? { linear tting on time interval [90; 300] years; the run
terminated half way due to numerical diôculties
it has not been shock-processed. Therefore, the code is con-
structed to follow the shock-enhanced chemistry and cooling
rather than the cold molecular gas. This will generally not
be a problem at the high densities where CO cooling and
gas-grain heating and cooling determine the properties of
the low-temperature gas.
The formula for H2 formation is critical to the results.
Reformation takes place on grain surfaces with the rate
taken from (Hollenbach & McKee 1979):
kR = 3  10 18 [cm 3 s 1 ]


T 0:5 fa
1 + 0:04 (T + Tg ) 0:5 + 2  10 3 T + 8  10 6 T 2 (5)
with a factor
fa = 1 + 10000 exp( 600=Tg )

1
; (6)
which means, with Tg = 20 K assumed, that fa is quite
close to unity in our simulations. Recent experimental re-
sults, summarised by Katz et al. (1999), indicate that the
H2 formation rate could be considerably lower. The inher-
ent grain mass and uniform space distribution adopted here
also in uence the results. This would imply that the abso-
lute reformation times derived here may need revision; we
intend to explore a range of possibilities in a following study.
The computations were performed on a Cartesian grid
with uniform initial density and periodic boundary condi-
tions in every direction, to simulate a region internal to a
molecular cloud. The hydrogen was xed to be fully molecu-
lar at the beginning (fraction: f = 0:5). Initial stress was in-
troduced by Gaussian perturbations applied to model veloc-
ities with a at spectrum extending over a narrow range of
wave numbers given by 3  jkj  4. We have performed runs
with root mean square velocities of 15 km s 1 , 30 km s 1
and 60 km s 1 , and computational domain sizes of 32 3 , 64 3 ,
128 3 and 256 3 zones. The number density has been taken
to be n = 10 6 cm 3 , physical box size L = 10 16 cm, and
the initial temperature has a homogeneous distribution of
T i = 100 K. Abundances for helium, free oxygen and car-
bon were taken as 0.1, 5  10 4 and 2  10 4 .
3 DATA ANALYSIS
3.1 Decay rates
The decay in total kinetic energy is displayed in Fig. 1. An
initial phase of evolution during which waves steepen and
velocity gradients transform into shocks is apparent. The
Table 2. Power law indices  of thermal energy decay. The 
coeôcients are derived from linear tting performed on the time
interval [90; 600] years. Average error 0:008.
resolution initial rms velocity for the run
15 km s 1 30 km s 1 60 km s 1
32 3 0:78 0:82 0:47 y
64 3 0:97 0:99 0:34 y
128 3 0:82 0:81 0:63 y
256 3 0:85 ? 0:75 0:62 y
? { linear tting on time interval [90; 300] years; the run termi-
nated half way due to numerical diôculties
y { linear tting on the time interval [140; 600] years, present nal
decay behaviour after initial strong cooling regime (see Fig. 2)
Table 3. Power law indices  of Mach number decay. The 
coeôcients are derived from linear tting performed on the time
interval [90; 600] years. Average error 0:005.
resolution initial rms velocity for the run
15 km s 1 30 km s 1 60 km s 1
32 3 0:08 0:13 0:46 y
64 3 0:16 0:15 0:55 y
128 3 0:21 0:16 0:49 y
256 3 0:23 ? 0:24 0:54 y
? { linear tting on time interval [90; 300] years; the run termi-
nated half way due to numerical diôculties
y { linear tting on the time interval [140; 600] years, present nal
decay behaviour after initial strong cooling regime (see Fig 3)
shock formation time is short, with a timescale
t i  L
2 km vrms
= 45:3  L
10 16 [cm]
10 [km s 1 ]
vrms
3:5
km
[yr]: (7)
where km is dened as the initial mean wavenumber (Mac
Low et al. 1998; Stone et al. 1998; Smith et al. 2000). Hence,
shocks develop faster in the 60 km s 1 turbulent eld, as
evidenced by the immediate power-law decay.
Figure 4. Average molecular fraction as a function of time. There
is virtually no dissociation during the run with an initial rms of
15 km s 1 { all four resolution curves merge in one straight line
to the 0.5 level. 30 km s 1 { intermediate case; 60 km s 1 {
practically all molecules are temporarily destroyed.
c
0000 RAS, MNRAS 000, 000{000

4 Georgi Pavlovski et al.
Figure 1. Average kinetic energy
as a function of time; log{log plot.
Figure 2. Average thermal en-
ergy as a function of time; log{log
plot.
Figure 3. Average Mach number
as a function of time; log{log plot.
c
0000 RAS, MNRAS 000, 000{000

Simulations of molecular turbulence 5
A rapid power-law decay of kinetic energy with time
was uncovered in the isothermal case (Stone et al. 1998; Mac
Low et al. 1998; Smith et al. 2000). This actually followed
an initial phase of evolution during which velocity gradients
increase and shocks appear. The time scale for decay is of the
order of a wave crossing time, which implied that turbulent
molecular clouds must be short lived if the turbulence is not
regenerated.
The rst result here is the conrmation that molecu-
lar turbulence behaves similar to isothermal turbulence, but
with higher rates of decay. The curves of average kinetic en-
ergy decay are shown in Fig. 1. Power laws deduced from
linear tting (log Ekin /  log t) of the curves are shown in
Table 1.
We have found that the formula
Ekin = E0(1 + t
t i
)  ; (8)
ts the data quite well. However, values of jj derived from
from tting the curve (8) are slightly larger due to nonzero
values of t i . This functional dependence explains why values
of the average kinetic energy at the end are compatible,
while at the beginning they are dierent by the factor of 4.
The power-law index  increases with vrms . This is con-
sistent with the high power-law index found for high Mach
number, M = 50, isothermal turbulence (Smith et al. 2000).
The second result is that the total thermal energy
(E th = nkT (3:3 f)=2) also decays with time , as opposed to
the imposed constant value for the isothermal case. The av-
erage thermal energy density behaviour is presented in Fig. 2
and the corresponding decay rates : E th / t  are shown
in Table 2. The thermal energy aids our description of the
physics. We have found that it also has a power-law decay
for the low speed case in which molecules are not destroyed.
The thermal energy decays, although not as fast as the ki-
netic energy. Thus the Mach number remains high relative to
the isothermal case, which may well account for the steeper
kinetic energy decay rate. Eventually, however, the gas cools
down to just below the grain temperature and the turbulent
heating input falls o. At this temperature range, however,
the numerical code is not accurate enough to resolve the
minimum.
The fact that the thermal energy exhibits a dependence
on numerical resolution (to a higher degree than kinetic en-
ergy does) is due to the importance of cooling layer res-
olution for the energy balance and smoothening eects on
coarser grids. This implies lower values of thermal energy on
ner grids, although a smoothening eect for the 64 3 run led
to higher (than the 32 3 run) values of the average thermal
energy.
With an initial rms speed of 60 km s 1 , the molecules
predominantly dissociate within 40 years. Cooling associ-
ated with the dissociation and shock heating are competing
processes (Fig. 2). After 40 years, there is a steep decay in
thermal energy as trace CO and H2O molecules eôciently
cool gas between 100 K and 8000 K. As conrmation, we
note that the number of zones within which the tempera-
ture is very high (> 10 3 K) falls from a percentage  45%
to  1% during the rst 150 yr. The gas which had become
almost fully atomic (see plots of average molecular fraction
in Fig. 4 and histograms on Fig. 5), is again over 60% molec-
ular hydrogen after 300 yr. As the molecules reform, the gas
is heated since the molecule energy released on reformation
is channelled collisionally into the gas (rather than being ra-
diated) at the high density of these simulations. Therefore,
the thermal energy subsequently decays slower than in all
other cases, approaching a shallow power-law value. These
changes are not directly re ected in the kinetic energy; how-
ever, the rapid gas cooling maintains a high Mach number,
which leads to the fast decay of kinetic energy.
For reference, the energies are plotted in erg cm 3 , with
the initial kinetic energy given by
E i = 1
2 v 2
rms = 1:2  10 6
 no
10 6 [cm 3 ]


 vrms
10 [km s 1 ]
 2
[erg cm 3 ] (9)
and the nal (thermal) energy, after complete molecule ref-
ormation ( = 10=7), is
Ef = 1
1 n k T = 1:6  10 9 no
10 6 [cm 3 ] 
Tf
10 [K] [erg cm 3 ]: (10)
Plots of the average Mach number are shown in Fig. 3, corre-
sponding decay rates, : Ma / t  are in Table 3. Note that
the display is logarithmic, with an average Mach number
exceeding 10 maintained for the rst 100 yr in the low-speed
case. Here we dene the average Mach number as
Ma =
s
hv 2 i
hc 2
s i ; (11)
where
c 2
s = (3:5 3f)(2:2 2f)
(3:3 2f) 2
e

is an adiabatic speed of sound, and the bracket notation
h: : : i means averaging over the domain: L 3
P
i;j;k . The
Mach number decaying behaviour can be qualitatively ex-
plained by noticing that it should be roughly proportional
to the fraction (Ekin=E th ) 1=2 , although the equality  =
( )=2 is not generally true.
The increase of the Mach number during the period
 [30; 120] yr in the runs with rms velocity of 60 km s 1 is
due to the strong dissociative cooling at that time, as shown
by the thermal energy decay curves, i.e. where jj > jj. Note
that this case, with the highest initial Mach number, pos-
sesses the lowest Mach number once the shocks have formed.
This is due to the dissociation, generating a warm atomic
gas. The question then is: why do the turbulent motions
decay so fast in this case, despite the low Mach number?
The suggestion is that the shocks are preferentially running
into the denser molecular gas, which can more eôciently
dissipate the energy because it is cooler. In isothermal gas,
however, this eect would not be apparent.
3.2 Molecular hydrogen evolution
The molecular fraction `f ' is displayed in Fig. 4. The three
initial states correspond to three distinct physical regimes.
With a rms velocity of 15 km s 1 , very few dissociative
shocks develop but some localised dissociation still occurs.
c
0000 RAS, MNRAS 000, 000{000

6 Georgi Pavlovski et al.
Figure 5. Histograms of the molecular fraction distribution. The black curve corresponds to the left y-axis and represents the volume
distribution of molecules: the number of zones with value of molecular fraction within the interval [f; f +
f ],
f = 10 2 . The grey
curve corresponds to the right y-axis and represents the mass distribution of molecules. The data were taken from the 128 3 run with
rms velocity of 60 km s 1 .
Figure 6. Histograms of the average molecular fraction distri-
bution with density at dierent times. Molecules reform quickly
even in regions with low densities. Data from 128 3 resolution run
with initial rms velocity of 60 km s 1 .
With 30 km s 1 , a few per cent of the molecules are dis-
sociated whereas at 60 km s 1 , the gas becomes over 80%
atomic.
Note that the minimum value of f is reached after
 4 t i . The dynamics alone determines the period of the
dissociation phase, which does not overlap with the refor-
mation phase for the chosen parameters.
Reformation of molecular hydrogen is unexpectedly
rapid. The expected H2 reformation time at 20 K and
10 6 cm 3 is tR = (kR no) 1 = 3,200 yr (Eq. 1) a factor of 5
larger than the simulation time. At 100 yr, the temperature
is  80 K, predicting a reformation time of tR = (kR no) 1 =
2,000 yr. Yet, reformation is occurring over  400 yr. This
speed up is caused by the turbulence itself: the molecules
preferentially reform in the denser and cooler locations. As
weak shocks propagate through the gas, dierent regions
are compressed and expanded. Hence the reformation time
is not only controlled by the `average' reformation time, but
also by the strength of the turbulence. Given a turbulent
dynamical timescale shorter than the average reformation
timescale, then we can expect reformation to be accelerated.
We have checked the degree of convergence of both the
dissociation and reformation processes. Molecule reforma-
tion is a gradual process and is, therefore, accurately rep-
resented in the simulations. High resolution is, however,
paramount to correctly follow the degree of dissociation.
This is critical to the 30 km s 1 turbulence since there are
numerous intermediate speed shocks, partially dissociating
the gas. For the low and high speed examples, however, the
dissociation is basically zero and complete, respectively. Ex-
haustive resolution studies of one-dimensional shocks with
this code are presented elsewhere (Smith & Rosen 2002).
The remarkable manner in which the molecular fraction
changes is illustrated in Fig. 5 for the high speed case. This
gure displays the number of zones Nz (f) with a molecular
fraction in each interval [f; f +
f ] where
f is taken as
10 2 . After 10 yr, there is an even distribution, suggesting
that dissociative shocks have in uenced about half the gas.
Up to 50 years, the majority of zones have low values of f .
Thereafter, the reformation has the following properties.
 A single maximum in Nz (f) shifts with time to higher
values of f .
 The maximum remains at Nz  10 5 zones.
 The width of the peak is roughly constant (full width
at 80% maximum is 0.16  0.02)
c
0000 RAS, MNRAS 000, 000{000

Simulations of molecular turbulence 7
Figure 9. The correlation length changes with time in runs with
dierent initial rms velocity. The data taken from the 128 3 run.
Thick lines represent correlation length, thin line represent inverse
average Mach number change.
 The distribution of f by mass displays a single wide
peak (grey line in Fig. 5).
 Even in regions with a very low density, reformed
molecules are found after 100 yr (Fig. 6).
These vital results imply that f is dominated by a small
range in values at any instant, at times after 100 yr. For
this reason, three dimensional structure and maps based on
column density of either all the gas or just the molecules
appear almost identical. In this sense, isothermal simulations
can be utilised to determine the structure of molecular clouds
or cores but not the underlying fraction of molecules.
3.3 Statistical analysis
In recent years, several studies of probability density func-
tions in numerical simulations have been advanced as a
step in their full statistical characterisation (e.g. Vazquez-
Semadeni & Garca (2001); Burkert & Mac Low (2001)). The
goal is to quantitatively compare observed maps of the inter-
stellar medium with the predictions of simulations. Here, we
present the rst analysis of the results for the column den-
sity of molecules from turbulent simulations which include
the molecular chemistry.
The autocorrelation function is dened as follows:
S2(x; y) =
L
X
i;j
(i; j) hi

(i + x; j + y) hi


L
X
i;j
(i; j) hi

2
; (12)
where (x; y) is the quantity (integrated along the LOS):
(x; y) = 1
L
L
X
z
(x; y; z); (13)
where (x; y; z) is the density (or molecular density), and
hi is an average:
hi = 1
L 2
L
X
x;y
(x; y): (14)
As a qualitative measure, however, it is more convenient to
use an averaged correlation function, i.e. polar angle inte-
grated:

S2 (r) = 2

Z =2
0
S2 (r cos ; r sin )r d; r 2 [0; L]: (15)
Having integrated the correlation function, we have lost all
information about anisotropic properties of the molecular
density map. However, in our studies we have found that
considerable deviations from polar symmetry of a correlation
function occur usually only in regions of weak correlation,
i.e. separated by a distance which is large in comparison with
the correlation length. This suggests that large scale struc-
ture which might seem to appear in molecular density maps
(e.g. see the map on Fig. 7) is not coherent, and therefore
turbulence in the simulations can be considered as isotropic.
We nd that the correlation length for the molecular
column density grows as the turbulence decays (see Fig. 9
also compare Fig. 7 and Fig. 8) as has been predicted for gas
in general by Mac Low (1999) and Mac Low & Ossenkopf
(2000). Dening the correlation length as the length over
which the average autocorrelation function has decayed by
a factor of 10 ( 
S2(`cr ) = 0:1), we nd that in the run with
resolution of 128 3 , the correlation length is related to the
average sound speed. The growth of `cr shows a good corre-
lation with the inverse average Mach number behaviour (see
Fig. 9). An explanation for this is easily seen from the fol-
lowing analysis. Dening the crossing time as tcross = L=v,
where  v is the average speed in the box, we nd that the
`information' during this time has travelled as much as
`cr = 
cs  tcross . From this, one immediately obtains the
following relation:
`cr
L =  cs
 v  Ma 1 : (16)
The autocorrelation functions for density and molecular
density maps indeed look very similar, as expected given
the above results for the evolution of the molecular fraction
within a relatively narrow range of values.
An analysis of the probability density functions (PDF)
reveals no particular dierence between PDFs for density
and molecular density maps and PDFs from dierent mo-
ments of time. Linear t coeôcients  of the histograms of
 in log-linear coordinates (N / (=max)  ) are distributed
c
0000 RAS, MNRAS 000, 000{000

8 Georgi Pavlovski et al.
Figure 7. Autocorrelation functions for molecular density and snapshot of molecular density map (from 256 3 data) at 100 yr. The
molecular density map is the result of the integration of the data cube (fraction  density) along z axis.
Figure 8. Autocorrelation functions for the molecular density and snapshot of the molecular density map (from 256 3 data) at 600 yr.
The molecular density map is the result of the integration of the data cube (fraction  density) along z axis.
in a quite wide range:   [4:0; 6:8] with a mean value of
5.13 and dispersion  0:66 (see Fig. 10).
It has been proposed that observational measurements
are often conned to within a correlation length (Burkert &
Mac Low (2001)). Linear t coeôcients of the average his-
tograms (mean from a hundred histograms) of regions which
have size of about a correlation length, yield a distribution
with a smaller dispersion ( 0:22) and have a mean value of
2.90 (Fig. 10). This is roughly consistent with the observed
value for molecular clouds found by Williams et al. (2000);
Blitz & Williams (1997). In their studies they have used
two-dimensional column density maps of clouds observed in
dierent radial velocity bins for an optically thin molecular
species. Each pixel in the data cube (galactic latitude, longi-
tude, radial velocity) has an associated antenna temperature
Ta , proportional to the column density of gas in that pixel.
The PDFs was then dened as the total number N of pix-
els with a certain column density =max , (where max is
the maximum column density found in the data cube). For
molecular clouds in Taurus the following asymptotic 1 value
of  is reached:  = 2:7  0:4. Although our simulations
are not directly comparable with the observations of Giant
Molecular Clouds, the consistency is quite remarkable.
3.4 Spatial Distribution
In Figs. 11 and 12, we present representations of the three di-
mensional distributions. These gures demonstrate how the
molecular density and total density distributions converge
between 50 and 100 yr. Initially, of course, while thermal
dissociation is occurring in the shock waves, there is less
correlation.
1  is found to be sensitive to the spatial smoothing, for smooth-
ing beam & 20 0 the PDF reaches an asymptotic value
c
0000 RAS, MNRAS 000, 000{000

Simulations of molecular turbulence 9
Figure 10. Probability distribution functions. The left gure presents histograms from the 60 km s 1 run at 100 yr. The right gure
presents histograms from the 30 km s 1 run at 100 yr. The thicker line corresponds to the averaged PDF from areas with a size of about
or less then one correlation length ( 0:23L), with the associated linear t coeôcient cr . The data are taken from the 256 3 run.
Figure 11. Molecular density distribution snapshots. Data taken
from the 128 3 resolution run with rms velocity 60 km s 1 . At
early stages (when average molecular fraction is low) the dis-
tribution of molecular density is dierent from mass density (see
Fig. 12); later, when hydrogen reforms, the distributions look very
similar
The general structure of the clouds is hard to dene.
There are some clumps and lament-like structures but little
visual evidence for sheets. Some clumps and laments are
surrounded by diuse envelopes, which are adjacent to other
laments and clumps. Self-gravity or sweeping by large-scale
shocks could generate more coherent structures.
Figure 12. Mass density distribution snapshots. Data taken from
128 3 resolution run with rms velocity 60 km s 1 . Compare the
distributions with the molecular density distribution on Fig. 11.
4 CONCLUSIONS
4.1 Summary
We have presented the properties of a specic model for
molecular turbulence. We carried out three dimensional hy-
drodynamical simulations of decaying supersonic turbulence
in molecular gas. We included a detailed cooling function,
molecular hydrogen chemistry and equilibrium C and O
c
0000 RAS, MNRAS 000, 000{000

10 Georgi Pavlovski et al.
chemistry. We studied three cases in which the applied ve-
locity eld straddles the value for which wholesale dissocia-
tion of molecules occurs. The parameters chosen ensure that
for the high-speed turbulence, the molecules are initially de-
stroyed in shocks and gradually reform in a distinct phase.
We nd the following.
 An extended phase of power-law kinetic energy decay,
as in the isothermal case, after an initial phase of slow dis-
sipation and shock formation.
 The thermal energy, initially raised by the introduction
of turbulence, decays only a little slower than the kinetic
energy.
 The reformation of hydrogen molecules, as the fast
turbulence decays, is several times faster than expected
from the average density. This is expected in a non-uniform
medium, as a consequence of the volume reformation rate
being proportional to the density squared.
 The molecules reform into a pattern of laments and
small clumps, enveloped in diuse structure.
 During reformation, the remaining turbulence redis-
tributes the gas so that the fraction of molecules is dis-
tributed relatively evenly. Hence, the density and molecular
density are almost identically distributed at any one time
after 100 yr.
 The correlation length in our simulations grows with
time as has been predicted for a gas in general.
 The probability density functions sampled from regions
with size of about one correlation length are consistent with
observations.
4.2 Discussion
We mainly wish here to emphasise the insight these simula-
tions provide into how molecular chemistry and supersonic
dynamics combine. We have found that isothermal simula-
tions are indeed very useful, not only for the rate of energy
decay but also to trace the molecules.
A simple reason for the fast decay is that a suôcient
number of strong shocks survive. As shown by Smith et al.
(2000), the rate of energy decay in decaying turbulence is
dominated by the vast number of weak shocks. These shocks
are less eôcient at energy dissipation. A second possible rea-
son is that the curved shock structures create small scale
vorticity, which leads to enhanced dissipation of kinetic en-
ergy. It is not clear, however, why relatively more vorticity
should be created when stronger shocks are present. Simu-
lations of isolated curved shocks would help clarify the dis-
sipation paths.
It is plausible that the high Mach number turbulence
could create more small-scale structure, leading to a faster
decay rates as argued by Mac Low (1999).
The simulations presented here may directly aid our
interpretation of fast molecular shocks in dense regions,
such as Herbig-Haro objects. For example, supersonic decay
would appear to occur in the wake of bow shocks such as
in DR 21 (Davis & Smith 1996). These simulations suggest
that the turbulence decays rapidly and bow shock wakes will
be bright but short.
The fast reformation of molecules indicates that
molecules may form on shorter time scales than previously
envisaged. Hence molecular clouds may appear out of atomic
clouds several times faster than anticipated in non-turbulent
scenarios. Recently, some evidence has been found for rapid
cloud formation and dissipation (Ballesteros-Paredes et al.
1999; Hartmann et al. 2001). This implies that much of the
interstellar medium may be undergoing both rapid dynam-
ical and chemical changes, driven by sources of supersonic
turbulence.
The simulations analysed here provide a basis for much
further work. For example, we can now investigate the prop-
erties of hydrodynamic clouds in which dust or chemical
abundances are non-uniform, clouds in which the turbulence
is driven uniformly or non-uniformly, and clouds which are
initially atomic.
5 ACKNOWLEDGMENTS
The computations reported here were performed using the
UK Astrophysical Fluids Facility (UKAFF) and FORGE
(Armagh), funded by the PPARC JREI scheme, in col-
laboration with SGI. M-MML was partially funded by the
NASA Astrophysical Theory Program under grant number
NAG5-10103 and NFS CAREER grant AST99-85392. AR
was funded by the PPARC. Armagh Observatory receives
funding from the Northern Ireland Department of Culture,
Arts and Leisure. ZEUS-3D was used by courtesy of the
Laboratory of Computational Astrophysics at the NCSA.
REFERENCES
Ballesteros-Paredes J., Hartmann L., Vazquez-Semadeni
E., 1999, ApJ, 527, 285
Balsara D. S., Crutcher R. M., Pouquet A., 2001, ApJ, 557,
451
Blitz L., Williams J. P., 1997, ApJ, 488, L145
Burkert A., Mac Low M.-M., 2001, Column Density Prob-
ability Distribution Functions in Turbulent Molecular
Clouds: A comparison between Theory and Observations,
arXiv:astro-ph/0109447
Davis C. J., Smith M. D., 1996, A&A, 310, 961
Hartmann L., Ballesteros-Paredes J., Bergin E. A., 2001,
ApJ, 562, 852
Heitsch F., Mac Low M.-M., Klessen R. S., 2001, ApJ, 547,
280
Herbst E., 2000, in Combes F., Pineau des For^ets G., eds,
Molecular Hydrogen in Space. The formation of H2 and
other simple molecules on interstellar grains.. Cambridge
University Press, pp 85{88
Hollenbach D., McKee C. F., 1989, ApJ, 342, 306
Hollenbach D., McKee C. F., 1979, ApJ, 41, 555
Katz N., Furman I., Biham O., Pirronello V., Vidali G.,
1999, ApJ, 522, 305
Klessen R. S., 2000, ApJ, 535, 869
Langer W. D., van Dishoeck E. F., Bergin E. A., Blake
G. A., Telens A. G. G. M., Whittet D. C. B., 2000, in
Mannings V., Boss A. P., Russel S. S., eds, The University
of Arisona Space Science Series, Protostars and Planets
IV. The University of Arizona Press, Tuscon
Lepp S., Shull J. M., 1983, ApJ, 270, 578
Lim A. J., 2001, MNRAS, 321, 306
c
0000 RAS, MNRAS 000, 000{000

Simulations of molecular turbulence 11
Lim A. J., Rawlings J. M. C., Williams D. A., 1999, MN-
RAS, 308, 1126
Mac Low M.-M., 1999, ApJ, 524, 169
Mac Low M.-M., Klessen R. S., Burkert A., Smith M. D.,
1998, Phys. Rev. Lett., 80, 2754
Mac Low M.-M., Ossenkopf V., 2000, A&A, 353, 339
McKee C., Storey J., Watson D., Green S., 1982, ApJ, 259,
647
Neufeld D., Kaufman M., 1993, ApJ, 418, 263
Padoan P., Juvela M., Bally J., Nordlund A., 1998, ApJ,
504, 300
Padoan P., Juvela M., Goodman A. A., Nordlund  A., 2001,
ApJ, 553, 227
Passot T., Vazquez-Semadeni E., Pouquet A., 1995, ApJ,
455, 536
Smith M. D., Mac Low M.-M., Heitsch F., 2000, A&A, 362,
333
Smith M. D., Mac Low M.-M., Zuev J. M., 2000, A&A,
356, 287
Smith M. D., Rosen A., 2002, The instability of Fast
Shocks in Molecular Clouds, Submitted to MNRAS;
http://star.arm.ac.uk/~rar/catinst/
Stone J. M., Norman M. L., 1992a, ApJS, 80, 753
Stone J. M., Norman M. L., 1992b, ApJS, 80, 791
Stone J. M., Ostriker E. C., Gammie C. F. F., 1998, ApJ,
508, L99
Sutherland R., Dopita M., 1993, ApJS, 88, 253
Suttner G., Smith M. D., Yorke H. W., Zinnecker H., 1997,
A&A, 318, 595
Vazquez-Semadeni E., Garca N., 2001, ApJ, 557, 727
Vazquez-Semadeni E., Ostriker E. C., Passot T., Stone
C. F., 2000, in Mannings V., Boss A. P., Russel S. S.,
eds, The University of Arisona Space Science Series, Pro-
tostars and Planets IV. The University of Arizona Press,
Tuscon
Vazquez-Semadeni E., Passot T., Pouquet A., 1996, ApJ,
473, 881
Williams J. P., Blitz L., McKee C. F., 2000, in Mannings
V., Boss A. P., Russel S. S., eds, The University of Arisona
Space Science Series, Protostars and Planets IV. The Uni-
versity of Arizona Press, Tuscon
APPENDIX A: THE CHEMISITRY AND THE
COOLING FUNCTION
We consider only a limited network of chemical reac-
tions, which has been tested in one-dimensional simulations
(Smith & Rosen 2002) and which is critical for molecular
cloud evolution. These reactions are:
H+H 7! H2 (A1)
H2 +H 7! 3H (A2)
H2 +H2 7! 2H +H2 (A3)
O+H2 7! OH+H (A4)
OH+H 7! O+H2 (A5)
OH+ C 7! CO +C (A6)
CO +H 7! OH+ C (A7)
OH+H2 7! H20 +H (A8)
H2O +H 7! OH+H2 (A9)
The rst three reactions (A1, A2, A3) are included in a semi-
implicit cooling-chemistry step to calculate the temperature
and H2 fraction. Reaction (A1) takes place on grain surfaces
with the rate given by equation (5) (see section 2.2).
We x free oxygen and carbon abundances, 0(O) and
0(C). We then calculate the equlibrium O, OH, and CO
abundances, assuming that the CO reactions are much faster
than the H2O reactions. Then, the remaning free oxygen is
distributed into O, OH, and H2O, according to equlibrium
abundances. In this manner we can calculate quite accu-
rately the chemistry within fast shocks without overloading
the hydrocode. For further details (see Smith & Rosen 2002,
Appendix B).
The cooling fuction used in the numerical code (see
equations 1 { 4 ) is composed of 11 separate parts (one
of which heats the gas):
 =
11
X
i=1
 i (A10)
The components are summarised in Table A1 below.
c
0000 RAS, MNRAS 000, 000{000

12 Georgi Pavlovski et al.
Table A1. Components of the cooling function
Formulae Description
 1 =  1  n 2 , where, on assuming standard dust properties,
 1 = 3:8  10 33:0 T 0:5 (T T dust (1:0 0:8 exp( 75=T )) [erg s 1 ]
gas-grain (dust) cooling (Hollenbach & McKee 1989). In
present calculations we x T dust = 20K.
2 = nH 2

L (h)
v
1+L (h)
v =L (l)
v
+ L (h)
r
1+L (h)
r =L (l)
r

, where the vibrational and
rotational coeôcients at high and low density are:
L (h)
v = 1:10  10 18 exp( 6744=T ) [erg s 1 ];
L (l)
v = 8:1810 13 exp(6840=T ) nH k H(0;1) + nH 2
k H 2 (0;1)


[erg s 1 ];
the terms k H(0;1) and k H 2 (0;1) are v : 0 ! 1 collisional excitation
rates,
k H(0;1) =
(
1:4  10 13 exp(T=125 T 2 =577 2 ); if T < Tv ;
1:0  10 12 T 0:5 exp( 1000=T ); if T > Tv ;
where Tv = 1635 K,
k H 2 (0;1) = 1:45  10 12 T 0:5 exp( 28728=(T + 1190));
L (h)
r =
(
dexf19:24 + 0:474x 1:247x 2 g; if T < Tr ;
3:9  10 19 exp( 6118=T ); if T > Tr ;
where Tr = 1087 K, and x = log(T=10000),
L (l)
r
Q(n)
=
(
dexf22:90 0:553 1:148x 2 g; if T < T l ;
1:38  10 22 exp( 9243=T ); if T > T l ;
where T l = 4031 K, and Q(n) = nH 2

0:77 + 1:2 (n H ) 0:77
collisional cooling associated with vibrational and
rotational modes of molecular hydrogen. These for-
mulae based on equations 7 { 12 in Lepp & Shull (1983)
 3 = (n H ) 2   2 collisional cooling of atoms. We have used Table 10 of
Sutherland & Dopita (1993) (with Fe = 0:5) for the form
of 2 and we added an extra thermal bremsstrahlung term
equal to 1:42  10 27 T 0:5 for T > 10 4 K
4 = (n H 2
+ 1:39n H )  nH 2 O  3 ; 3 = 1:32  10 23 (T=1000)  ,
where  = 1:35 log(T=1000).
cooling through rotational modes of water induced
by collisions with both atomic and molecular hydro-
gen. Given values of  ts the values tabulated by Neufeld
& Kaufman (1993)
5 = 1:03  10 26 nH 2 nH 2 OT exp( 2325=T ) exp( 47:5=T 1=3 ) cooling through vibrational modes of water induced
by collisions with molecular hydrogen (Hollenbach &
McKee 1989)
6 = 7:40  10 27 nH3nH 2 OT exp( 2325=T ) exp( 34:5=T 1=3 ) cooling through vibrational modes of water induced
by collisions with atomic hydrogen (Hollenbach & Mc-
Kee 1989)
7 = 7:18  10 12 (n H 2
) 2 kD;H 2
+ nHnH 2
kD;H


; where dissociation
coeôcients are taken to be,
kD;H = 1:210 9 exp( 52400=T ) (0:0933 exp( 17950=T ))  [cm 3 s 1 ];
kD;H 2
= 1:310 9 exp( 53300=T ) (0:0908 exp( 16200=T ))  [cm 3 s 1 ];
 =
h
1:0 + n

2f

n 1
2 n 1
1

+ n 1
1
i 1
critical densities are t by the following approximation, n 1 =
dexf4:0 0:416x 0:327x 2 g [cm 3 ] and n 2 = dexf4:845 1:3x +
1:62x 2 g [cm 3 ].
cooling from the dissociation of molecular hydrogen.
Factor 7:18  10 12 [erg] is the 4.18 eV dissociation energy;
n 1 { is the density critical for dissociation by collision of
molecular hydrogen with atomic hydrogen, n 2 { with itself.
8 = 7 , where  7 = kR (see Eq. 5) and  = nnH (1 )7:1810 12 heating resulting from the reformation of molecu-
lar hydrogen. We employ  to parametrise the fraction of
released energy which is thermalised rather than radiated.
9 = nCOnkTv T = 1 + na=ncr + 1:5(na=ncr ) 0:5


, where the mean
speed of the molecules v T =
p
8kT=(mH 2
) and ncr = 3:3 
10 6 (T=1000) 0:75 [cm 2 ],  = 3:0  10 16 (T=1000) 0:25 [cm 2 ]. Av-
erage density, na = 0:5(n H + nH 2
p
2).
cooling through rotational modes of carbon monox-
ide induced by collision with both atomic and
molecular hydrogen. We based our equations on
Eqs. 5.2{5.3 in McKee et al. (1982)
 10 = 1:83  10 26 nH 2
nCOT exp( 3080=T ) exp( 68=T 1=3 ) cooling through vibrational modes of carbon
monoxide induced by collisions with molecular hy-
drogen (see Neufeld & Kaufman 1993)
11 = 1:28  10 24 nHnCOT 0:5 exp( 3080=T ) exp( (2000=T ) 3:43 ) cooling through vibrational modes of carbon
monoxide induced by collisions with atomic hydro-
gen
c
0000 RAS, MNRAS 000, 000{000