Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://star.arm.ac.uk/preprints/401.pdf Äàòà èçìåíåíèÿ: Wed Sep 17 19:19:16 2003 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 22:06:37 2012 Êîäèðîâêà: Ïîèñêîâûå ñëîâà: m 63

Mon. Not. R. Astron. Soc. 339, 133­147 (2003)

The instability of fast shocks in molecular clouds
Michael D. Smith and Alexander Rosen
Armagh Observatory, College Hill, Armagh BT61 9DG

Accepted 2002 October 7. Received 2002 September 2; in original form 2001 October 19

ABSTRACT

We report on the discovery that moderately fast shocks in dense molecular clouds with low transverse magnetic fields are likely to be unstable. The instability is triggered by the promoted cooling that results from the formation of carbon monoxide and water molecules in an extended warm shock section. Numerical methods are employed to demonstrate that, in the absence of magnetic fields, the instability regime is restricted to densities above n 0 = 104 cm-3 , velocities between 30­70 km s-1 , and O or C abundances above 10-4 , so that cooling from reforming molecules dominates in the warm gas without being suppressed by ultraviolet dissociation. The result is either a quasi-periodic or chaotic collapse and re-establishment of the warm shock layer on a typical time-scale of 106 cm-3 /n 0 yr with variations on shorter time-scales and changes in period being possible. Infrared emission lines from the unstable region, including the H2 lines, exhibit orders of magnitude variability. Atomic lines such as H display constant fluxes but undergo rapid radial velocity variations. Key words: hydrodynamics ­ instabilities ­ molecular processes ­ shock waves ­ ISM: lines and bands ­ ISM: molecules.

1 INTR ODUCTION Fast shock waves alter the dynamical, physical and chemical properties of dense interstellar clouds (Hollenbach & McKee 1979). The shocks may be driven internally or externally by a variety of sources, including supernova blast waves, protostellar jets, cloud­cloud collisions, stellar winds and planetary nebula. We define a fast shock as one capable of destroying molecules in its path. In a dense cloud, many of the molecules may subsequently reform in the cooling layer of compressed and accelerated material (Hollenbach & McKee 1989; Neufeld & Dalgarno 1989a). The compressed cloud may then go on to form stellar systems or may fragment and disperse. Hence, the behaviour of fast shocks is critical to the distributions of both stars and interstellar gas. Here we investigate the time-dependent properties of fast shocks through one-dimensional (1D) numerical simulations. The structure and signatures of steady fast shocks was studied in depth by Hollenbach & McKee (1989), Neufeld & Dalgarno (1989a,b) and Wolfire & Konigl (1991). If a shock is unstable, however, we have ¨ to reconsider not only the observed emission properties of radiative shocks but also (i) the transfer into turbulent energy of the cloud, (ii) the fragmentation or destruction of the cloud, (iii) the resulting shock-induced molecular abundances and (iv) the accuracy of numerical simulations that do not fully resolve individual shocks. Hydrodynamic shocks that destroy molecules at the shock front possess speeds in excessvadjust of 23 km s-1 (Kwan 1977;
E-mail: mds@star.arm.ac.uk (MDS); rar@star.arm.ac.uk (AR)
C

Hollenbach & McKee 1980; Smith 1994). In dense molecular gas shielded from ultraviolet (UV) radiation, shocks can be continuous (C-type) in which a low ion fraction leads to the drift of the magnetic field through the neutral gas. In this case, molecules largely survive in shocks with speeds up to 30 km s-1 if the magnetic field is parallel (Smith 1993) but up to 50 km s-1 if the field is transverse (Draine, Roberge & Dalgarno 1983). Below these speeds, the shock front and the cooling layer overlap. Above these speeds, however, the shock is classified as fast since molecules are abruptly dissociated after the shock front, followed by a wide zone of radiative cooling. We demonstrate here that moderately fast shocks in dense molecular clouds are thermally and dynamically unstable. Thermal instabilities on small scales are expected in the cooling layer that follows the shock front provided that the cooling rate per unit mass increases as the temperature falls (McCray, Kafatos & Stein 1975). This is because a small region within the flow with a lower temperature than its surroundings will cool faster, resulting in a growing density contrast (see Smith 1989). To also be dynamically unstable, the cooling rate must increase sufficiently rapidly as the temperature falls (Langer, Chanmugam & Shaviv 1981; Chevalier & Imamura 1982). On a certain time-scale, a cooling shocked layer may collapse simultaneously as the shock front shrinks back, reducing the shock strength. The lower shock velocity then promotes faster cooling, which leads to an even faster collapse of the cooling layer. Thus, a resonance occurs and the entire layer may collapse `catastrophically'. The lowered shock velocity also leads to a lower pressure. Hence, a resisting medium or confining wall can decelerate the retreat of the shock front. This can lead to

2003 RAS


134

M. D. Smith and A. Rosen
to follow. We compare this with the cooling functions for steadystate shocks with non-equilibrium chemistry and conclude that the instability will probably be present, but 1D simulations with full time-dependent chemistry and ionization are needed to verify this conjecture and determine the minimum abundances necessary. Note also that the molecular hydrogen formation rate remains quite uncertain and may also influence the instability regime. In Section 3, we present numerical results that demonstrate the nature of the instability and yield the time-scales. The strong influence of the instability on observational properties, with the emphasis on infrared emission lines, is then discussed (Section 4.1). Several conditions must be fulfilled for the instability to operate. First, sufficient CO and H2 O must form as the gas cools below 8000 K. Hollenbach & McKee (1989) have shown that CO and water cooling do indeed dominate for steady shocks in high-density regions. We consider here the carbon and oxygen depletions that stabilize the shock. Ameliorating factors also include thermal conduction and dust depletion. Furthermore, a sufficiently strong transverse magnetic field cushions the shock and inhibits the dynamical instability (Toth & Draine 1993). ´ Note that UV radiation will influence the location in which the molecules reform (Hollenbach & McKee 1989). For shocks with speed v 0 80 km s-1 , the UV radiation will inhibit the formation of CO and H2 O molecules. However, the calculation of Neufeld & Dalgarno (1989a) shows that OH cooling then replaces CO cooling, providing a similar steep increase in cooling as the temperature falls. Nevertheless, we restrict our results to fast shocks below this speed. Other shock instabilities may also be operating. These mainly distort or fragment the dense slab of gas that accumulates either between two shocks or between a single shock and a medium of high-thermal pressure. The instabilities include the Vishniac thinshell linear instability (Vishniac 1983; Mac Low & Norman 1993), the non-linear thin-shell instability (Vishniac 1994) and the transverse acceleration instability (Dgani, van Buren & Noriega-Crespo 1996). They all require two-dimensional studies, whereas the oscillatory instability is present even in one-dimensional simulations. Moreover, this instability occurs in the cooling layer rather than in the accumulated cold slab, thus possessing distinct signatures. 2 THE FLO W E Q U A TIONS AND AN AL YTICAL PREDICTION 2.1 The equations We model the most basic radiative shocks: hydrodynamic flows in one dimension. Numerically, we will solve the time-dependent flow equations: ( v ) + =0 t x ( v ) ( v 2 + p ) + =0 t x e (ev ) v + = -p - t x x (T , n , f ) (1) (2) (3)

a period of increasing shock velocity and a cooling layer of increasing temperature as the shock layer inflates. The higher temperature gas does not cool efficiently, and so expands outwards, driving the shock front to even higher speeds. The shock may perform growing periodic oscillations that lead to high-amplitude quasi-periodic oscillations (e.g. Imamura, Wolff & Durisen 1984; Strickland & Blondin 1995). This dynamical overstability is relevant to fast atomic shocks. For velocities greater than approximately 140 km s-1 , strong atomic cooling near 100 000 K is predicted to initiate the overstability (e.g. Innes, Giddings & Falle 1987; Gaetz, Edgar & Chevalier 1988). An instability was not expected, however, in molecular shocks since all the individual cooling functions are positive functions of temperature, as shown in Fig. 1. We were therefore surprised to find a dynamical instability in our numerical simulations. To test a new version of a three-dimensional (3D) molecular hydrodynamic code, we performed one-dimensional high-resolution tests. This code improves the code of Suttner et al. (1997) also including CO and H2 O chemistry and cooling (see the Appendices). The previous code did not show any signs of dynamical instability. In fact, it is the production of CO and H2 O molecules in the cooling layer that leads to the onset of catastrophic cooling. The formation of molecules causes the instability by altering the cooling function that is now time-dependent as gas moves downstream. Note that Silk (1983) recognized that the rapid formation of water molecules could lead to thermal instabilities in molecular clouds. The fundamental condition for the dynamical instability has been derived through numerical solutions for specific cases, especially for a monatomic gas with a cooling rate proportional to the square of the density, with atomic shocks being considered (Chevalier & Imamura 1982; Chanmugam, Langer & Shaviv 1985; Dgani & Soker 1994). Numerical simulations demonstrate that the instability leads to highamplitude oscillations (Strickland & Blondin 1995). In Section 2, we present the hydrodynamic equations and discuss the criteria for which molecular shocks should be unstable, the unstable modes and their oscillation periods. We also present the cooling as a function of temperature for steady-state shocks with the equilibrium C and O chemistry as employed in the 3D numerical calculations for which time-dependent oxygen and carbon chemistry is not practical

Figure 1. The total cooling (thick line) and cooling functions for a slab with the indicated fixed density and H2 fraction. The cooling components are: 1, gas­grain; 2, H2 ro-vibrational; 3, atomic; 4, H2 O rotational; 5, H2 O vibrational with H2 ; 6, H2 O vibrational with H; 7, H2 dissociation; 8, H2 reformation heating (thin double line); 9, CO rotational; 10, CO vibrational with H2 ; 11, CO vibrational with H; 12, [O I] fine structure; 13, OH rotational.

( fn ) ( fn v ) + = R (T , n , f ) - D (T , n , f ), (4) t x where n is the hydrogen nuclei density, e is the internal energy density and f is the molecular hydrogen abundance i.e. n (H2 ) = fn . We take a helium abundance n (He) = 0.1n . Therefore, the total particle density is n tot = (1.1 - f )n and the temperature is T =
C

2003 RAS, MNRAS 339, 133­147


Shocks in molecular clouds
p /(kn tot ). The internal energy loss through radiation and chemistry per unit volume is , and the dissociation and reformation rates of molecular hydrogen are given by D and R, respectively. The internal energy is related to the pressure by e = p /( - 1) where the effective specific heat ratio is taken as = 5.5 - 3 f , 3.3 - f (5)

135

which assumes that H2 possesses two rotational degrees of freedom. We define 1.4 µ= = . (6) n tot m p 1.1 - f

2.2 Power-law cooling predictions Shock stability studies have been mainly restricted to flows with = 5 and power-law cooling of the form 3 =
0

n T



(7)

(Chevalier & Imamura 1982; Wolff, Gardner & Wood 1989; Strickland & Blondin 1995) and superpositions of power-law cooling (see Saxton et al. 1998 for further references). Thermal instabilities in the interstellar medium were first analysed in detail by Field (1965). Smith (1989) investigated the small-wavelength isobaric condensation modes within shock waves. The instability condition (equation 24 of Smith 1989) reduces to < 2( - + 1)/(3 - ) at the shock front and approaches the condition < downstream. Note that for = 2, the instability condition remains close to < 2, whereas for = 1 the condition is more stringent. The oscillatory instability condition is even more stringent, depending upon the cooling processes and the boundary conditions. For the case = 2, the oscillations are unstable in the fundamental mode for 0.4 and in higher modes for c where the values of c exceed 0.8 (Chevalier & Imamura 1982). Note, however, there is apparently a range of values 0.29 < < 0.5 for which non-linear effects are critical (Strickland & Blondin 1995). In this range, secondary shocks form near the wall and restrict, but do not halt, the collapse. The fundamental mode is thought to be the most significant since the long-wavelength perturbations would be least damped by motions transverse to the flow direction. Numerical simulations demonstrate, however, that higher-order modes still generate oscillations for 0.7 (Strickland & Blondin 1995). For applications to shock propagation through molecular gas, we should consider = 7 (fully molecular) and = 10 (fully molec5 7 ular with 10 per cent helium). Furthermore, is unity for several cooling processes in high-density gas, e.g. CO rotational and H2 radiative cooling. For fast dissociative shocks, however, we find that the hydrogen is still almost completely atomic during the cooling down to below 500 K, and that the instability is caused by CO and H2 O vibrational cooling (with the molecules excited by collisions with H) for which = 2. In these respects, existing analyses may provide relevant criteria for the present conditions. 2.3 The steady-state model The cooling function for a uniform slab of gas with constant density and chemistry is displayed in Fig. 1. The component cooling functions are listed in Appendix A and chemical reactions in Appendix B. They were selected for the problems associated with shocks in dense molecular clouds. Fig. 1 assumes the equilibrium chemistry described in Appendix B. We conclude that the cooling functions
C

Figure 2. The cooling functions corresponding to dissociative shock waves that prove to be unstable, assuming equilibrium chemistry for C and O, with standard conditions as shown and abundances ( O) = 5 â 10-4 and (C) = 2 â 10-4 (upper panel) and (O) = 2.1 â 10-4 and (C) = 2 â 10-4 (lower panel) and Alfven speed 0.01 km s-1 . The lower panel corresponds ´ to the unstable depleted O case. The individual cooling/heating components are labelled as in Fig. 1.

are positive functions of temperature with the power-law index being above the oscillatory instability critical value of 0.4 in each section. This holds not only for the indicated parameters but under much wider conditions. The unstable shocks that motivated this analysis involve a complex multicomponent time-dependent cooling function (Fig. 2). Immediately following the shock front, is a zone of very steep cooling, dominated by H2 radiative and dissociative cooling and atomic cooling. Therefore, we expect this hot zone with T > 8000 K to be stable. We find that this thin section is maintained throughout, moving en bloc to accommodate the oscillating cool layer, with signals travelling quickly through it. The signal propagates slowly through the following cool layer, which we refer to as the `infrared layer', and this determines the dynamical instability properties. The steady-state model solves the hydrodynamic equations in the post-shock flow (with the partial time derivatives being dropped). The shock itself is replaced by a discontinuity satisfying the Rankine­Hugoniot conditions, with attention being paid to the specific heat ratio. The pre-shock gas parameters are represented by p 0 , v 0 , T 0 , 0 , 0 , f 0 and M 2 = ( 0 v 2 0 / p 0 ). Since molecular dissociation requires a finite time, the immediate post-shock values are f 1 = f 0 , 1 = 0 ,

2003 RAS, MNRAS 339, 133­147


136
S=
1

M. D. Smith and A. Rosen
/
0

= v0 /v1 = 1- 1 S

(0 + 1) M 2 , (0 - 1) M 2 + 20 M
2

(8)

p1 / p0 = 1 + and T1 / T0 = p1 . Sp0

(9)

(10)

The flow in the radiative layer, on eliminating e, is described by v=
1 v1

(11)
2 1 v1

p + v 2 = p1 +

(12)

and, on substituting, v [ p /( - 1)] v + p = - (n , T , f ), x - 1 x (13)

in which, on using equation (11), we can substitute from -4.4 f = x (3.3 - f )2 x and f 1.4m p ( R - D ) = x 1 v1 (15) (14)

to leave a single first-order equation that can be solved numerically without difficulty. Before running the steady shock model and simulations, we added equilibrium oxygen and carbon chemistry in addition to the hydrogen chemistry. The basic reactions that form and react with the OH molecule were chosen (Appendix B). Previously, we had taken the CO abundances to follow the molecular hydrogen fraction and not included H2 O, as the simplest approximation that allowed us to model molecular jets in three dimensions (Suttner et al. 1997). Here, however, we assume that the conversion of oxygen into OH and then into CO and H2 O can be described by equilibrium abundances through reactions with H, H2 and C. This is especially critical to large-scale 3D simulations since it is advantageous to avoid the introduction of new variables. We thus find that CO and H2 O can form ahead of the main H2 reformation zone since a small fraction of H2 O is sufficient to instigate the complete conversion of the free oxygen. Note that OH and oxygen fine structure cooling are also included in the steady-state and slab models, but excluded from the simulations. We find that these cooling contributions can be neglected in the present shock simulations for densities exceeding 3 â 104 cm-3 , as demonstrated for a hydrogen nucleon density of 106 cm-3 in Fig. 2. Time-dependent chemistry has been included in the steady-state calculations to test the equilibrium assumption. Equilibrium chemistry is found to provide a reasonable estimate of the cooling function at high densities. As demonstrated by comparing the panels of Figs 2 and 3, the cooling generated by H2 O formation is weaker in the nonequilibrium case. This would influence the oscillatory structure but we do not believe that this is sufficient to remove the instability. The instability is still present when either C or O is depleted by over an order of magnitude, yielding a more gradual increase in the cooling rate (e.g. see the lower panel of Fig. 2 where the H2 O bump is

Figure 3. The cooling functions within dissociative shock waves with nonequilibrium chemistry and with standard conditions as described in Fig. 2. The individual cooling/heating components are labelled as in Fig. 1.

greatly reduced). We also note that to predict the detailed structure, dynamical processes are just as critical as chemical effects. Twodimensional simulations with full time-dependent chemistry will be needed. For a density of 104 cm-3 , however, non-equilibrium chemistry is certainly required to reproduce the important contributions of oxygen fine structure and water cooling below 8000 K. As shown in Fig. 4, for non-equilibrium chemistry, the rapid increase in CO vibrational cooling as the gas cools below 8000 K is no longer mirrored in the total cooling. Dynamical instabilities may occur when the molecular hydrogen fraction is inversely related to the temperature, such as is possible in the cooling layer of a fast shock, as now discussed. The cooling is presented per unit volume, , in Figs 2 and 3 since the cooling functions are not simple power-law functions of the density. Hence, to predict which of these states would be unstable, we need to interpret the critical temperature indices discussed in Section 2 in terms of the index given by f (T ) = 0 T = 0 T - , which follows a fluid parcel through the cooling layer at almost constant pressure. The interpretation contains uncertainty, however, since the chemistry (and, hence, the cooling) is time dependent. Nevertheless, it is clear that a very high negative value for is present in the temperature regime between 3000 and 8000 K for standard non-depleted O and C abundances. Moreover, the width of the shock is mainly determined by the temperature at which the rapid cooling is about to turn on, i.e. near 8000 K (as confirmed in Fig. 5). This location is adjacent to the cooling zone that possesses a high inverse
C

2003 RAS, MNRAS 339, 133­147


Shocks in molecular clouds

137

Figure 4. A comparison of equilibrium and non-equilibrium chemistry cooling functions within dissociative shock waves at the low density of 104 cm-3 and with standard conditions as described in Fig. 2. The individual cooling/heating components are labelled as in Fig. 1.

dependence with temperature. This leads us to conclude that the shock is, at least, linearly unstable. An example with O and C both strongly depleted confirms the cause of the instability. Low oxygen and carbon abundances permit little H2 O, CO and OH to form. As shown in Fig. 6, the shock should be stable since -1.39 except for a narrow temperature range, as confirmed by the simulations (lower panel of Fig. 12; see Section 3.6). 3 HYDR OCODE SIMULA TIONS 3.1 The numerical model:
ZEUS

Figure 5. Profiles of temperature, density, pressure and molecular fractions (H2 , CO, and H2 O) demonstrating the collapse of the shocked region. We display profiles from the simulation 6F40 covering two stages of a cycle. This shows that the collapse in the shock covers temperatures of 500­ 5000 K and is initiated by lower pressures (seen here as the dip in density prior to the collapse) in this region.

-3D

We employ the ZEUS-3D code in a one-dimensional mode (Stone & Norman 1992) to provide the basic hydrodynamics. This is a secondorder Eulerian finite-difference code. Here we study compressible hydrodynamics without gravity, self-gravity or thermal conduction. No physical viscosity is modelled, but numerical viscosity remains present, and a von Neumann artificial viscosity determines the dissipation in the shock front. A module for molecular chemistry and molecular and atomic cooling has been added. The functions and rates are listed in the Appendices. Our ultimate goal is to develop a reliable code with which we can tackle three-dimensional molecular dynamics, later adding gravity, magnetic fields, ambipolar diffusion and radiation.
C

We thus restrict the cooling and chemistry lists to just those items essential to the dynamics. We have employed the simultaneous implicit method discussed by Suttner et al. (1997) in which the timestep is adjusted in order to limit the change in internal energy in any zone to 30 per cent. We present primarily two-shock simulations. In addition, we have examined single shocks, with a reflection boundary condition simulating a wall. Symmetric inflow boundary conditions were chosen at both ends of the grid in the two-shock simulations. Hence, a dense slab of cold material accumulates at the centre of the grid. We then find that the behaviour of each shock is independent and the accumulated slab material behaves exactly as a wall. Initially, however, there is no effective wall or shock. As a direct result, the twin shock flow pattern might not display the oscillatory instability since waves are not brought into resonance through reflection. Thermal instabilities are present, however, and rapidly generate some structure and asymmetry.

2003 RAS, MNRAS 339, 133­147


138

M. D. Smith and A. Rosen
in atomic shocks, are caused by the nature of the cooling, which possesses a time dependence through the molecular abundances. We typically require a shock resolution of 180 to obtain the physical instability. It is significant that the instability generates some chaotic behaviour, with the different cooling components able to produce different resonances. We thus find, similarly to studies of turbulent flow, that the behaviour converges but not the exact structure, as the resolution is increased. Other notable points are as follows. (i) When a low resolution is employed, the shock is still unstable but unphysical long-term high-amplitude variations are found. (ii) We increased the period at which the shock width is recorded after 100 yr. Along with this change, more short-period variations become apparent. Hence, the spiky structure that appears during the first 100 yr of the shock evolution curves are caused by a low display resolution rather than numerical resolution. (iii) The behaviour does not converge to a single pattern but passes through various cyclic patterns. This may be caused by the presence of different harmonic modes and their interaction through the complex cooling function.

Figure 6. The cooling functions within a stable dissociative shock wave, with low abundances (O) = 2 â 10-5 and (C) = 1 â 10-5 and other conditions as in Fig. 2. The individual cooling/heating components are labelled as in Fig. 1.

3.2 Resolution and convergence Once the central dense slab has built up, the two shocks usually begin to behave independently and the dynamical instability sets in. Table 1 includes the data for our detailed resolution study and our exploration of parameter space for the standard two-shock case, and a shock­wall example. We define the shock resolution as the ratio of the mean shock total cooling length to the zone length. Fig. 7 displays the locations of the shock front (as defined by the positions of the maximum temperature) for one side of the shock as a function of time. This demonstrates the presence of an instability at all resolutions. It is, however, only at resolutions at or above 4 â 1010 cm that the variation time-scale and amplitude approach a consistent behaviour. We remark that the oscillations can change period at irregular intervals. We propose that these changes, not encountered
Table 1. Shock inlet parameters. Sim.a 6B60 6D60 6B40 6C40 6D40 6D40-C N 500 1000 500 500 1000 1000 x 8(10) 2(10) 3.2(11) 1.6(11) 8(10) 8(10) logn 6 6 6 6 6 6 v
s

3.3 The collapse and re-expansion We display an excerpt in Fig. 8 from data recorded with a time resolution of 0.1 yr. In the displayed era, short variations on a 1yr time-scale are superimposed on large variations with a time-scale of 4 yr. Although most collapses and expansions exhibit a smooth behaviour, some collapses take place at high speed. This is displayed explicitly in Fig. 9 where we plot the distribution of the shock front velocity, as calculated between each consecutive data pair. The asymmetry demonstrates that (i) the collapses are often, but not

1-/2-shock 2 2 2 2 2 2

Cat. unstable Yes Yes Yes Yes Yes Yes

Period M: 2.1, 2.7, W: 0.3, 3.6, 6.9

b

60 60 40 40 40 40

6D40-O 6E40-CO 6E40 6F40 6F40-1 6A20 6D20 5E40 5F40 4E40 4F40
a

1000 1000 2000 5000 5000 500 1000 1000 2000 1000 5000

8(10) 8(10) 4(10) 2(10) 2(10) 8(10) 1(10) 3.2(11) 1.6(11) 2.56(12) 1.28(12)

6 6 6 6 6 6 6 5 5 4 4

40 40 40 40 40 20 20 40 40 40 40

2 2 2 2 1 2 2 2 2 2 2

Yes No Yes Yes Yes No No Yes Yes Yes Yes

t < 500 yr: S: 21.7, M: 4.5, 8.7, 49.8 500 < t < 1000 yr: M: 35.6, W: 5.0, 9.4­16.6, 11.3 1000 < t < 1500 yr: M: 21.7, 49.8, W: 11.3 550 < t < 750 yr: S: 10.4 S: 1.4, 0.7, M: 0.47, W: 0.27, 0.34, 2.43 S: 16 W: 1.3, 3.0, 3.6 M: 0.42, 1.1, 1.2, 1.3, 1.8, 2.6, 3.2 W: 0.3, 0.6, 0.7, 2.6 W: 1.0, 1.1, 1.2, 1.4, 1.6 M: 5.8, 6.6, 34.5, W: 9.0, 9.8, 13.6 S: 370

The designation for each simulation lists, in order, the log of the number density, a representation of the resolution, and the pre-shock velocity in km s-1 . The letters for the resolution are associated with the number of zones in the mean shock total cooling length, as follows: A, 8­16; B, 16­32; C, 32­64; D, 64­128; E, 128­256; and F, 256­512. b The period in years determined using an FFT, with W, M, S indicating a weak, moderate and strong peak, respectively.
C

2003 RAS, MNRAS 339, 133­147


Shocks in molecular clouds

139

Figure 7. The dependence of the instability on spatial resolution. We plot the position of maximum temperature (for clarity, of one shock) as a function of time for simulations (from the top down) 6B40, 6C40, 6D40, 6E40 and 6F40 that span five different spatial resolutions. The zone size, x , is displayed in each panel. The frequency of data displayed is identical for each panel, we show one datum per year for the first 100 yr, and thereafter one every 0.1 yr.

Figure 8. Expanded view of short-term oscillations in a high-resolution case. We plot the position of maximum temperature as a function of time for both shocks, where one is shown as a dotted line, in the high-resolution simulation 6F40. The positions are plotted every 0.1 yr. This shows the similar, if somewhat offset in time, behaviour of each independent shock.

exclusively, at high speed, in the range -3 to -8 km s-1 and (ii) the most probable shock front velocity is larger than the average shock front speed by just +1kms-1 . Note that these speeds are an order of magnitude smaller than the shock speed, yet they are part of large-amplitude variations in shock length. The distances of both shock fronts from the centre of the grid are displayed in Fig. 8. Since we are investigating an instability, there is no reason why the shock front location should coincide.
C

Figure 9. The probability distribution of shock velocities in the highresolution case. We plot the relative fraction of velocities, which are determined from consecutive times and positions, of each of the two shocks for the high-resolution simulation 6F40, and for the portion of the simulation with 0.1-yr data intervals. The width of the velocity bins is 0.2 km s-1 .

3.4 Shock profiles Fig. 5 displays the profiles of temperature, density, pressure and molecular fractions at different phases during one collapse. The upper panel demonstrates that the collapse indeed coincides with the shrinking of the warm infrared layer. This layer determines the

2003 RAS, MNRAS 339, 133­147


140

M. D. Smith and A. Rosen
employed a fast Fourier transform (FFT) on the highest-resolution runs and present in Table 1 some periods that appear. Since the oscillatory pattern changes over time, many periods that show up in the FFT analysis are very weak (<3 , denoted by W in the Table). Still, these values do show an inverse relationship between the period and the pre-shock density. Specifically, we see periods of - approximately 3.5n 6 1 yr, where n 6 is n /106 cm-3 , for the range of densities and simulations shown in Fig. 10. A study of the dependence of the shock width on the shock speed is presented in Fig. 11 for the case with density of 106 cm-3 . The instability is present at all speeds above 25 km s-1 , where the hydrogen molecules are first destroyed immediately behind the shock front, and rapid reformation leads to dynamical instability. We have not investigated shocks faster than 60 km s-1 since, at such speeds, UV radiation from the hot gas delays the formation of molecules in the infrared layer. Nevertheless, the instability could still occur but deeper in the infrared layer. A code that includes radiative transfer effects will be necessary to confirm this. Note that we have altered the time-step for the display at 100 yr. In these two-shock simulations, this allows the shocks to progress to the point where each is independent of one another and also enables us to display the structure that would be recorded through discrete observations made every few years. The manner in which the shock structure changes depends strongly on the upstream speed imposed. The 60 km s-1 example occasionally displays narrow shock regions, roughly with a 30-yr time interval. Cyclic patterns are uncovered at 40 km s-1 .At20km s-1 the shock is stable. 3.6 Chemistry The oscillatory instability was not found in previous simulations of fast shocks (Suttner et al. 1997) because, owing to the lack of computing power, oxygen and carbon chemistry were not included. Yet these prove to be vital to the nature of a fast shock. This is

length of the shock since the gas cools less efficiently at temperatures below 8000 K, until some CO molecules have formed. At the end of the collapse, the remaining infrared layer consists of only 8000 K gas. On the other hand, both the hot and cold gas occupy very thin regions. Note that a pressure trough appears when the shock is wide. This low pressure results in the shock collapse. In turn, this produces a weak secondary shock that returns upstream. The logarithmic profile of the hydrogen molecular fraction maintains a simple shape since very little H2 reforms in the infrared layer. Fast reformation does not occur until the gas has cooled down below 300 K, which is close to the interface between the two shocks. The shape of the H2 O fraction follows that of H2 quite closely. The CO, in contrast, fully reforms in the infrared layer. Note the generation of a secondary shock (at x =-3 â 1012 cm in the figure) during the collapse that is sufficiently strong to destroy most of the reforming CO molecules. Although the shock is weak in comparison to the primary shock, it occur in the regime where the H2 and CO chemistry is most sensitive to the temperature. 3.5 Density and velocity Fig. 10 shows that the instability is present at all densities considered here. The code assumes that equilibrium chemistry remains a valid approximation and that no other cooling mechanism has been neglected. At a density of 104 cm-3 , however, the assumption of equilibrium cooling is invalid and the slow formation of CO combined with the enhanced O fine structure cooling, leads to a much flatter cooling function. For our standard case of 40 km s-1 , we note that the highest amplitude and greatest variation time-scale occur for the lowest density. The time-scale dependence is expected since the lower densities have longer cooling times. This is apparent from the difference in shock width shown in Fig. 10. Quasi-periodic oscillations of the shock front are found in some stages of the simulations. We have

Figure 10. The dependence of the dynamical instability on density. The position of maximum temperature as a function of time for three different densities is plotted for simulations 6E40, 5F40, and 4F40 (from the top down). As with Fig. 7 that displays two-shock simulations, we show the position of only one of the two shocks. Here, the data time interval is 2 yr in 6E40, every 2.5 yr in the first 250 yr and every 0.25 yr for all t > 250 yr in simulation 5F40, and every 10 yr in 4F40. Since the size of the shock is inversely proportional to the density, the vertical scale of each panel is different. The time-scale of the oscillations increases as the density is decreased.

C

2003 RAS, MNRAS 339, 133­147


Shocks in molecular clouds

141

Figure 11. (pre-shock) the position every 0.5 yr

The dependence of the instability on pre-shock velocity. We plot the position of maximum temperature as a function of time for 3 different initial velocities in the simulations 6D60, 6E40, and 6D20 (from the top down). As with the previous position figures of two-shock simulations, we show of only one of the two shocks. Here the data are shown once per year for t < 100 yr and once per 0.1 yr for t > 100 in 6D60 and in 6E40, and in simulation 6D20. Note the different vertical scales, and that the 40 km s-1 simulation has the largest shock region.

confirmed in Fig. 12, which displays the evolution of the shock width for four different O and C abundance combinations. Note that the shock is stable only when both O and C are depleted. Both H2 O and CO formation generate the instability, as expected from the shape of the cooling function in the steady shock analysis. It is evident that when just one molecule is responsible for the instability (middle panels), a regular oscillation is generated. This suggests that the changes in pattern in other simulations are mainly caused by the multipeaked cooling function. The small-amplitude oscillations uncovered for the oxygen-depleted shock are related to the inefficiency of the CO to continue to cool the gas below 2000 K. Note that when the instability is instigated through water molecule formation and cooling, a high-amplitude oscillation is apparent. Interesting changes in temporal behaviour are found for this carbondepleted shock simulation as shown in detail in Fig. 13. A dramatic change occurs at t = 500 yr. It seems plausible that the two shocks prior to t = 500 yr interact through the cold, dense, thin wall, and we suspect the change in period, to effectively half of the previous one, signifies that the wall has become wide enough to allow the shocks to behave independently. 3.7 Single shock impacting a wall We have also verified that a wall model produces the same unstable structures as in the equivalent two-shock model, as shown in Fig. 14. Hence, our conclusions are independent of the chosen shock set up, given the limitations of one-dimensional shock waves. In two dimensions, the wall that forms between two shocks may well bend and fragment owing to induced transverse motions. 4 IMPLICA TIONS 4.1 Observational aspects We now calculate the effect of the instability on observable quantities. First, we chose molecular hydrogen emission lines to represent
C

the signatures from the infrared layer. The upper energy level associated with the H2 1 0 S(1) transition corresponds to an excitation temperature of nearly 7000 K. Emission arising from a cascade following molecule reformation is not accounted for here (we assume that the energy released on reformation is efficiently redistributed into the thermal component by collisions). Fig. 15 displays in grey-scale showing the location of the emission as a function of time. The integrated emission is shown in the panels on the right. We also show the 2 1 S(1) emission, which arises from more excited gas. The 2 1 S(1) spatial distribution, however, is very similar to 1 0 S(1) and so is not illustrated. The top panels demonstrate the characteristics of a steady shock, produced by taking low abundances. Note that the central wall slowly builds up in time with a constant velocity, producing the expected triangular shape. The H2 emission is spread over a large region, since the shock is dominated spatially by the slowly cooling gas between 4000 and 8000 K. The integrated emission of the H2 lines is constant, with the diagnostic ratio R = F [1 0 S (1)]/ F [2 1 S (1)] 16. This is consistent with the range of values from R = 2.5 for the standard steady case to R = 25 for the reduced abundances example [( O) = 1.0(-4) and (C) = 2.3(-4)] tabulated by Hollenbach & McKee (1989). The bottom row of Fig. 15 shows an example of large variations in the integrated flux of the H2 line in a simulation where the instability is active. The variations in 1 0 S(1) are over an order of magnitude. This is caused by the rapid appearances and disappearances of the whole 2000­4000 K section. This causes the ratio R to vary in the range 1.5­30. The time-scale for this variation is less than 2 yr, with some significant variations in a time-scale of months, in this high-density example. The CO flux from a low-lying rotational level (dashed line) is overwhelmingly generated in the cold downstream gas at the wall, and thus is simply related to the accumulated cold mass (Fig. 15). Variations in the H atomic line flux (not displayed) are present but relatively small since the hot dissociated part of the cooling layer is preserved throughout. The shock front speed does vary, however,

2003 RAS, MNRAS 339, 133­147


142

M. D. Smith and A. Rosen

Figure 12. The dependence of the instability on abundances of carbon and oxygen. The position of maximum temperature as a function of time for four combinations of nominal and depleted C and O abundances is plotted. In the top panel, simulation 6E40 is presented, which has the default abundances of (C) = 2.0(-4) and (O) = 5.0(-4). In the second panel (top to bottom), simulation 6D40-O is displayed, which has abundances of (C) = 2.0(-4) and (O) = 2.1(-4) and thus almost all of the oxygen will end up in CO and not H2 O. In the third panel is simulation 6D40-C with abundances of (C) = 1.0(-5) and (O) = 5.0(-4). In the final panel is simulation 6E40-CO with abundances of (C) = 1.0(-5) and (O) = 2.0(-5). Note that the simulations in the lower three panels were run with similar zone sizes, but that the decreasing cooling rate leads to larger cooling length-scales and a change in the resolution in the final case. The top two panels display data every 0.1 yr, and in the bottom two panels every 2 yr.

Figure 13. Change of period and amplitude of the instability in the carbon-depleted shock simulation. We display the position of one shock in simulation 6D40-C to late times, with 2-yr intervals between data. Note the dramatic change in the character of the instability at t = 500 yr. Just prior to this there is a series of paired oscillations, with one shallow dip followed by a deeper one. During this time period (350 < t < 500 yr), the other shock is completely out of phase with that shown, with its own shallow and deep collapse.

and this generates moderate changes in flux. A more recognizable signature of the instability is the pattern of variation of the mean radial velocity of the H line when viewed along the shock normal. The radial velocity goes through sporadic variations of over 10 km s-1 , accompanied by large dips and rises in flux. As with the integrated H2 line emission above, large variations occur in just a few months at a density of 106 cm-3 but take 30­50 yr at a density of 104 cm-3 .

4.2 Physical interpretation The complexity and time dependence of the total molecular cooling function generate complex unstable patterns. Cooling (atomic and H2 dissociative) following the shock front is extremely rapid, and can be treated as part of the shock front itself. These processes produce a narrow spike in the temperature profile (see Fig. 5). We have thus carried out a linear analysis by modelling the shock front
C

2003 RAS, MNRAS 339, 133­147


Shocks in molecular clouds

143

Figure 14. Variation of shock position for a single-shock simulation. Here we display the position of maximum temperature in simulation 6F40-1, which has only one shock and a reflecting boundary at the wall. We plot data once a year for t < 130 yr, and every 0.1 yr thereafter. Note the instability is still present and the variation of the shock position has a period and an amplitude similar to the two-shock simulation 6F40.

Figure 15. Spatial and temporal variations of molecular emission lines. We display the results for two simulations, with data every 2 yr for simulation 6E40-CO and every 0.1 yr for 6F40 in the bottom row. In each row, the first panel displays the spatial variation of the emission from H2 1 0 S(1), and in the second the variation of the CO 0 0 R(1) emission. In these two plots, darker shading indicates more emission, so most of the emission in H2 is from the part of the cooling layer closest to the central wall. Similarly, the bulk of the CO emission, contributions to which include a maximum speed cut-off at 15 km s-1 , is from the central wall. In the rightmost panel of each row, we show the temporal variation of the integrated emission from the two lines mentioned above(H2 1 0, solid; and CO, long dashed) and for the H2 2 1 S(1) line (short dashed).

and the rapidly cooling layer by jump conditions with a range of effective specific heat ratios for the shock front of 10­40. Hence, the modelled front treats the shock heating, atomic cooling and dissociative cooling as one discontinuity. For the subsequent cooling layer we set = 5 . A linear instability analysis solves the equations 3 as presented by Chevalier & Imamura (1982) but with the revised shock front boundary conditions. We find that the rapid cooling of the hot gas has a stabilizing influence on the shock (Smith & Rosen 2002), as expected from a previous linear analysis that included a second steep component (Saxton et al. 1998). We indeed find that negative values of the power-law index are required to produce linear instability (Smith & Rosen 2002). For example, with a shock compression factor of S = 20, corresponding to = 21/19, the first overtone grows for -0.4 and the fundamental mode only for -4. We thus expect to find significant regimes in which higher harmonics dominate.
C

The many different manifestations of the instability have been explored and classified by Walder & Folini (1996). They propose five different oscillation types. In particular, they find that the first overtone exists in a strong form that results in a high-amplitude quasi-periodic collapse and re-expansion of the shock with the appearance of weak subshocks (M-type). They also uncover a smooth sinusoidal mode (S-type). Some similarity of our long-term simulations with their simulations suggests that the first overtone may be responsible for much of the behaviour recorded here. Our simulations, however, generally demonstrate complex behaviour with several superimposed non-stationary wave forms usually being present. The simple oscillations encountered in the oxygen-depleted shock simulation have been examined for the harmonic signatures. Here, a simple wave form is present that reveals the harmonic mode. When the shock is at its maximum width, a high-amplitude velocity wave grows strongly near the wall. The wave represents infalling

2003 RAS, MNRAS 339, 133­147


144

M. D. Smith and A. Rosen
The application of the results are limited by the following assumptions. First, the magnetic field must be either parallel to the shock direction or must have a transverse component with an Alfven speed ´ 1 km s-1 for the cooling function to remain unchanged above 2000 K. Shock flows with transverse magnetic fields characterized by pre-shock Alfven speeds of 2.4 km s-1 (shock speed of 30 km s-1 ) ´ down to 1.0 km s-1 (shock speed of 70 km s-1 ) maintain the gradient in the cooling that initiated the instability in the oxygen-depleted simulation. We thus suspect that these flows will also exhibit the instability. This needs to be verified with magnetohydrodynamic (MHD) simulations. Secondly, the one-dimensional analysis does not allow for other instabilities that require transverse motions to be present. These instabilities depend on the specific problem and could well make the signatures more complicated to interpret. In this respect, the specific signatures from the 1D results should allow us to determine whether observed variations are caused by the instability. Thirdly, time-dependent abundances of H2 formation and dissociation are included but equilibrium CO and H2 O is assumed in the simulations. We have shown here that the cooling function derived with a non-equilibrium chemistry is smoother than, but has a shape approximately similar to, that assuming equilibrium for pre-shock densities n in excess of 3 â 104 cm-3 . Low-resolution numerical simulations should be treated with caution. When the sharp gradients in abundances are smeared out, apparently stable shock fronts may result. This is in addition to the coarse mesh problems noted by Walder & Folini (1996), which can eliminate particular modes. The instability occurs in high-density and low-speed dissociative shocks. For shocks with speed v 0 > 70 km s-1 , the UV radiation will inhibit the reformation of CO and H2 O molecules (Hollenbach & McKee 1989). Hence, we restrict our results to fast shocks below this speed, although we note that unstable conditions may also arise through OH cooling for v 0 80 km s-1 (see fig. 12 of Neufeld & Dalgarno 1989a). Large variations in shock-generated emission lines have been well documented. Recently, Chrysostomou et al. (2000) found variability of the 1 0 S(1) H2 line over 4 years in several locations within Herbig­Haro (HH) objects. Herbig & Jones (1981) documented 21 years of optical structural variations in HH 1. In the future, large flux variations will be increasingly documented as higher spatial resolution will permit observations of coherent shock sections. What would make it possible to associate these variations with a radiative cooling instability? A series of simultaneous spectroscopic observations of emission lines in the optical and infrared would provide a comprehensive data set. In the absence of such information, we have shown here that the interpretation of shocks as steady flows may lead to misleading results. A CKNO WLEDGMENTS The numerical calculations were run on computer, acquired through the PPARC participation. AR is funded by PPARC. many suggestions that helped us improve REFERENCES
Biham O., Furman I., Pirronello V., Vidali G., 2001, ApJ, 553, 595 Chanmugam G., Langer S.H., Shaviv G., 1985, ApJ, 299, L87 Chevalier R.A., Imamura J.N., 1982, ApJ, 261, 543 Chrysostomou A., Hobson J., Davis C.J., Smith M., Berndsen A., 2000, MNRAS, 314, 229
C

material and is accompanied by a strong mass influx on to the wall. This strong rarefaction wave propagates outwards, until it meets the collapsing shock front, at which point the collapse is halted. The shock front re-expands, accompanied by a much smoother velocity profile. This sequence repeats throughout the simulation. This behaviour is a characteristic of the first overtone, caused by the velocity (Chevalier & Imamura 1982) and mass inflow (Walder & Folini 1996) being out of phase with the shock front position. The oscillating modes have been described as standing waves, with zero velocity at a reflecting wall (the dense cold layer) and shaken at the shock front (Chevalier & Imamura 1982). The wave pattern is analogous to that in a half-open organ pipe (Saxton et al. 1998). This yields a series of resonances with the angular frequencies n = 2(2n + 1)/ t0 , (16)

where t 0 is the period of the fundamental mode. Calculated periods are 20.6­22.6 and 6.7­7.1 and 4.2, respectively, in units of L /v 0 , derived from a linear analysis (Chevalier & Imamura 1982) for in the range -1to + 1 and shocks with a compression factor of 2 S = 4. We derive the relationship t0 2( S + 1)/1.6 from a linear analysis (Smith & Rosen 2002); i.e. the (linear) instability period is approximately increased by the compression factor. Now we wish to determine whether the numerically derived results are consistent with the theoretical expectations. A very detailed comparison is not attempted here because the total cooling function is complex. The total cooling at the beginning of the unstable shock section at 6000 K is estimated as L ir 10-24 n 2 (v/40 km s-1 )2 0 erg cm-3 s-1 (see Fig. 2). The density at this location is n ir n 0 m p v 2 /(kT ) = 32.2 n 0 (v 0 /40 km s-1 ). Hence the cooling time-scale 0 at 6000 K is 1.5 n ir kT / L ir 120 (104 cm-3 /n 0 ) yr. This agrees in order of magnitude with the variation time-scales in the simulations. The oscillation dynamical time-scale is defined as t d = L c /v ir , where v ir is the speed of the gas at 6000 K and L c 7 â 1014 (104 cm-3 /n 0 ) cm is measured from the simulations. This yields t d = 200 (104 cm-3 /n 0 ) yr. Hence, the relevant cooling, dynamical and measured variation time-scales are all comparable. The theoretical linear oscillation period described above is, for the first harmonic, 1.3 SLc /v 0 = 190 (104 cm-3 /n 0 ) yr on substitution for S and L c . Hence the linear analysis is consistent with the numerical results.

5 SUMMAR Y AND CONCLUSIONS We have shown that the cooling layers of hydrodynamic molecular shocks with speeds between 30 and 70 km s-1 are subject to cyclic, potentially large-scale collapse and reformation. This is owing to a form of catastrophic cooling that produces an overstability caused by the rapid formation of CO and H2 O molecules. When C and O are both depleted, the instability is absent. Quasi-periodic oscillations dominate some sections of the shock simulations. The time-scale of the variations may also change. Periodic oscillations are present when either CO or H2 O dominate the cooling, as controlled by the C and O abundances. Large variations are predicted from those lines produced in the gas with temperatures in the range 1000­8000 K. We have shown this explicitly for near-infrared H2 lines but it also applies to high-J CO lines. Quite large radial velocity variations are expected from emission lines produced in gas moving with the shock front, such as those associated with recombination and common tracers of weak hydrodynamic shocks in the optical range, such as S[II].

the local FORGE superJREI initiative with SGI We thank the referee for the paper.

2003 RAS, MNRAS 339, 133­147


Shocks in molecular clouds
Dgani R., Soker N., 1994, ApJ, 434, 262 Dgani R., van Buren D., Noriega-Crespo A., 1996, ApJ, 461, 927 Draine B., Roberge W., Dalgarno A., 1983, ApJ, 264, 485 Field G.B., 1965, ApJ, 142, 531 Gaetz T.J., Edgar R.J., Chevalier R.A., 1988, ApJ, 329, 927 Herbig G.H., Jones B.F., 1981, AJ, 86, 1232 Hollenbach D., McKee C.F., 1979, ApJS, 41, 555 Hollenbach D., McKee C.F., 1980, ApJ, 241, L47 Hollenbach D., McKee C.F., 1989, ApJ, 342, 306 Imamura J.N., Wolff M.T., Durisen R.H., 1984, ApJ, 276, 667 Innes D.E., Giddings J.R., Falle S.A.E.G., 1987, MNRAS, 226, 67 Jacquet R., Staemmler V., Smith M.D., Flower D.R., 1992, J. Phys. B, 25, 285 Kwan J., 1977, ApJ, 216, 713 Langer S.H., Chanmugam G., Shaviv G., 1981, ApJ, 245, L23 Le Bourlot J., Pineau des Forets G., Flower D.R., 1999, MNRAS, 305, 802 ^ Lepp S., Shull J.M., 1983, ApJ, 269, 560 Mac Low M.-M., Norman M.L., 1993, ApJ, 407, 207 McCray R., Kafatos M., Stein R.F., 1975, ApJ, 196, 565 McKee C.F., Storey J.W.V., Watson D.M., Green S., 1982, ApJ, 259, 647 Neufeld D.A., Dalgarno A., 1989a, ApJ, 340, 869 Neufeld D.A., Dalgarno A., 1989b, ApJ, 344, 251 Neufeld D.A., Kaufmann A., 1993, ApJ, 418, 263 Pineau des Forets G., Flower D.R., Dalgarno A., 1988, MNRAS, 235, 621 ^ Saxton C.J., Wu K., Pongracic H., Shaviv G., 1998, MNRAS, 299, 862 Shapiro P.R., Kang H., 1987, ApJ, 318, 32 Silk J., 1983, MNRAS, 205, 705 Smith M.D., 1989, MNRAS, 238, 235 Smith M.D., 1993, ApJ, 406, 520 Smith M.D., 1994, MNRAS, 266, 238 Smith M.D., Rosen A., 2002, Ap&SS, submitted Stone J.M., Norman M.L., 1992, ApJS, 80, 753 Strickland R., Blondin J.M., 1995, ApJ, 449, 727 Sutherland R.S., Dopita M.A., 1993, ApJS, 88, 253 Suttner G., Smith M.D., Yorke H.W., Zinnecker H., 1997, A&A, 318, 595 Toth G., Draine B.T., 1993, ApJ, 413, 176 ´ Vishniac E., 1983, ApJ, 274, 152 Vishniac E., 1994, ApJ, 428, 186 Walder R., Folini D., 1996, A&A, 315, 265 Whitworth A., Clarke C., 1997, MNRAS, 291, 578 Wolff M.T., Gardner J., Wood K.S., 1989, ApJ, 346, 833 Wolfire M.G., Konigl A., 1991, ApJ, 383, 20 ¨

145

monoxide induced by collisions with molecular hydrogen and 11 is cooling through vibrational modes of carbon monoxide induced by collisions with atomic hydrogen. Therefore, we have excluded OH cooling ( 13 ) and [O I] fine structure cooling ( 12 ) from the present simulations. The form for dust (or gas­grain) cooling is taken from equation (2.15) in Hollenbach & McKee (1989):
1

= n 2 1 ,
-33

(A2)

where, on assuming standard dust properties, 1 = 3.8 â 10 T
1/2

(T - Tdust )
-1

â [1.0 - 0.8 exp(-75/ T )] erg s

cm3 .

(A3)

In the present calculations, we fix T dust = 20 K, and hence assume that the dust cools very rapidly after the shock front, as justified by Whitworth & Clarke (1997). H2 vibrational and rotational cooling is based on equations (7)­ (12) in Lepp & Shull (1983). We have shock tested these formula against the detailed subroutines presented by Le Bourlot, Pineau des Forets & Flower (1999) and found no significant differences to the ^ detailed shock properties. We thus take
2

=n

H2

L vH L rH , + 1 + ( L vH / L vL ) 1 + ( L rH / L rL )

(A4)

where the vibrational cooling coefficients at high and low density are (the first formula corrects a print error): L L
vH vL

= 1.10 â 10

-18 -13

exp(-6744/ T ) erg s-1 , exp(-6840/ T )

(A5)

= 8.18 â 10

â [n H kH (0, 1) + n H2 kH2 (0, 1)] erg s-1

(A6)

and the terms k H (0, 1) and kH2 (0, 1) are the v = 0 1 collisional excitation rates [and the exp (-6840/T ) term converts these to deexcitation rates], as follows: kH (0, 1) = 1.4 â 10 1.0 â 10
-13 -12

exp[(T /125) - (T /577)2 ] T
1/2

T < Tv T > Tv , (A7)

exp(-1000/ T )

APPENDIX A: THE COOLING FUNCTION Here we describe the cooling function used in the numerical code (equation 3) and in the present version of the steady model outlined in Section 2.3. The function is composed of 11 separate parts (one of which heats the gas), plus two further components were introduced in the steady shock models.
11

where T v = 1635 K, and kH2 (0, 1) = 1.45 â 10
-12

T

1/2

exp[-28728/(T + 1190)]. T < Tr T > Tr ,

(A8)

The rotational cooling rate coefficient at high density is L
rH

=

dex[-19.24 + 0.474x - 1.247x 2 ] if 3.90 â 10
-19

exp(-6118/ T )

if

(A9)

=

i i =1

.

(A1)

where T r = 1087 K, and x = log (T /10 000 K). For low density, the coefficient is
L rL Q (n )

The components are as follows: 1 is gas­grain (dust) cooling, 2 is collisional cooling associated with vibrational and rotational modes in molecular hydrogen, 3 is collisional cooling of atoms, 4 is cooling through rotational modes of water induced by collisions with both atomic and molecular hydrogen, 5 is cooling through vibrational modes of water induced by collisions with molecular hydrogen and 6 is for vibrational modes of water induced by collisions with atomic hydrogen, 7 is cooling from the dissociation of molecular hydrogen, 8 is heating resulting from the reformation of molecular hydrogen, 9 is cooling through rotational modes of carbon monoxide induced by collisions with both atomic and molecular hydrogen, 10 is cooling through vibrational modes of carbon
C

=

dex[-22.90 - 0.553x - 1.148x 2 ] if 1.38 â 10 n
0.77 H2 -22

T < Tl T > Tl ,

exp(-9243/ T )
0.77

if

(A10)

where T l = 4031 K, and Q (n ) = + 1.2 (n H ) . (A11)

Atomic cooling takes the expected form
3

= n 2 2 , H

(A12)

where we have used table 10 of Sutherland & Dopita (1993) (with Fe = -0.5) for the form of 2 and where an additional thermal bremsstrahlung term equal to 1.42 â 10-27 T 1/2 is added to 2 for T > 10 000 K.

2003 RAS, MNRAS 339, 133­147


146

M. D. Smith and A. Rosen
For CO vibrational cooling from collisions with H2 , we use (Neufeld & Kaufmann 1993)
10

Water rotational cooling takes the form:
4

= (n

H2

+ 1.39n H )n (H2 O)
-23

3

(A13) (A14)

3 = 1.32 â 10

(T /1000) ,



= 1.83 â 10

-26

n H2 n

CO

T exp(-3080/ T ) exp(-68/ T

1/3

). (A28)

where = 1.35­0.3 log(T /1000) fits the values tabulated by Neufeld & Kaufmann (1993). Cooling via water vibrational modes in collision with H2 and H is given by Hollenbach & McKee (1989):
5

For CO vibrational cooling through collisions with atomic hydrogen, we use
11

= 1.28 â 10 âT
1/2

-24

nHn

CO 3.43

= 1.03 â 10

-26

n H2 n (H2 O)
1/3

exp(-3080/ T ) exp[-(2000/ T )

].

(A29)

6

âT exp(-2325/ T ) exp(-47.5/ T = 7.40 â 10-27 n H 3 n (H2 O) âT exp(-2325/ T ) exp(-34.5/ T

), ).

(A15) (A16)

1/3

The cooling from the dissociation of molecular hydrogen follows the Shapiro & Kang (1987) modifications to the Lepp & Shull (1983) form:
7

Oxygen cooling through the 63-µm fine structure line has been included in the steady shock model and we intend to include it in further simulations in a form convenient for large-scale dynamical computations of molecular turbulence. We omit the calculation of the upper energy level that generates the much weaker 146-µm line. The cooling is then taken as
12

= 7.18 â 10

-12 -12

n

2 H2 k D,H2

= 2.82 â 10

-18

n

O

+ n H n H2 k

D,H

,

(A17)

1 , 1/ f H + A10 / r L

(A30)

where 7.18 â 10 k
D,H

erg is the 4.48 eV dissociation energy,

where A10 = 8.95 10-5 is the spontaneous transition rate, fH = 0.6 exp(-228/ T ) 1 + 0.6 exp(-228/ T ) + 0.2 exp(-326/ T ) (A31)

= 1.2 â 10

-9

exp(-52400/ T )
-1

â[0.0933 exp(-17950/ T )] cm3 s k
D,H2

(A18)

= 1.3 â 10

-9

â[0.0908 exp(-16200/ T )] cm3 s 1 1 - n2 n1 1 + n1

exp(-53300/ T )

is the fractional occupation of the 3 P1 level in local thermodynamic equilibrium and r L = rH + rH2 with rH = [4.37 â 10
-12

-1

T

0.66

(A19)

0.6 exp(-228/ T ) 0.2 exp(-326/ T )](n H + 0.48n H ) (A32)

with
-1

+ 1.06 â 10 and r
H2

-12

T

0.80

=

1.0 + n 2 f

(A20)

and finally, the critical densities for dissociation by collisions of molecular hydrogen with atomic hydrogen, n 1 , and with itself, n 2 , are fitted by the following approximations: n 1 = dex(4.0 - 0.416x - 0.327x ) cm
2 -3

= [2.88 â 10

-11

T

0.35

0.6 exp(-228/ T ) 0.2 exp(-326/ T )]n
H2

+ 6.68 â 10

-11

T

0.31

(A33)

(A21) (A22)

n 2 = dex(4.845 - 1.3x + 1.62x 2 ) cm-3 .

Heating from the reformation of hydrogen is including as a `cooling' term, in the following form:
8

are the collisional rates (values kindly provided by D. Flower), including rates calculated by Jacquet et al. (1992) for H2 collisions from which we have taken a single weighted index to combine ortho and para values for expediency, with an ortho to para ratio of 3, sufficiently accurate for the present purpose. OH cooling is taken as (see Hollenbach & McKee 1989)
13

= -7 ,
-12

(A23)

= 2.84 â 10

-28 2

nT

3/2

.

(A34)

where the rate 7 = k R is given in appendix B and = nn H (1 - )7.18 â 10 . (A24) Hence we employ to parametrize the fraction of released energy that is thermalized rather than radiated. Cooling via CO rotational modes through collisions with molecular or atomic hydrogen have been described in McKee et al. (1982). We base our cooling on their equations (5.2)­(5.5),
9

APPENDIX B: THE CHEMISTR Y We consider a limited network of chemical reactions, which we can test in these one-dimensional simulations and which is not too cumbersome to be included in three-dimensional MHD self-gravity simulations. These reactions are: (1) H + H H2 (2) H2 + H 3H (3) H2 + H2 2H + H2 (4a) O + H2 OH + H (4b) OH + H O + H2 (5a) OH + C CO + H (5b) CO + H OH + C (6a) OH + H2 H2 O + H (6b) H2 O + H OH + H2 . The first three reactions are included in a semi-implicit coolingchemistry step to calculate the temperature and H2 fraction.
C

=n

CO

n

kT vT 1 + n a / n cr + 1.5(n a / n cr )
.75

1/2

,

(A25) 8kT /(m H2 ) and (A26) (A27)

where the mean speed of the molecules vT = n cr = 3.3 â 106 T30 = 3.0 â 10
-16

cm-3 , cm-2 ,

T3

-1/4

where T 3 T /1000 K. We also have set the average density, n a = = 0.5(n H + 2n H2 ).

2003 RAS, MNRAS 339, 133­147


Shocks in molecular clouds
Reaction (1) takes place on grain surfaces with the rate taken from Hollenbach & McKee (1979): k R = (3 â 10 â
-18

147
(B4)

(4b) = 6.60 â 10

-13

1.53 T300 exp(-2970/ T ) cm3 s-1 ,

cm3 s-1 ) T
1/2

which yield an equilibrium of fa + 210-3 T + 810-6 T (B1) (O) n (H) -0 = 0.28T300 .4 exp(970/ T ) . (OH) n (H2 ) (B5)

1 + 0.04(T + Tg )1/2

2

with a factor f a = [1 + 10 000 exp(-600/ Tg )]-1 , (B2)

In comparison, the formula of Pineau des Forets et al. (1988) yields ^ (O) n (H) = 0.45 exp(1030/ T ) . (OH) n (H2 ) (B6)

which means that in our simulations, f a is quite close to unity. Since this formula is subject to ongoing debate (e.g. Biham et al. 2001), one can also alter f a to test the sensitivity. Dissociation rates (2) and (3) are given in Appendix A. We fix free oxygen and carbon abundances, 0 (O) and 0 (C). We then calculate the equilibrium O, OH and CO abundances, assuming that the CO reactions are much faster than the H2 O reactions. Then, the remaining free oxygen is distributed into O, OH and H2 O, according to equilibrium abundances. In this manner, we can calculate approximately the chemistry within fast shocks without overloading the hydrocode. At temperatures below 200 K, equilibrium abundances are unlikely to be reached within cooling layers. The main error will be in the H2 O abundance: in gas cooling below a few hundred Kelvin, the code will predict that the oxygen will be in O, whereas the H2 O that forms in the shock would remain as water considerably longer. CO and gas­grain cooling dominate, however, between 1000 and 8000 K so that the shape of the cooling function in this region is still approximately correct. Predicted line strengths of O and H2 O from the cold gas, however, will be in error. We find only small differences in published rates for the six reactions (4)­(6) (e.g. Hollenbach & McKee 1989; Pineau des Forets, ^ Flower & Dalgarno 1988). From Hollenbach & McKee (1989), the two rates (4) are (4a ) = 2.32 â 10
-12 1.93 T300

With the abundance of OH relative to O determined, the CO abundance is then given through (5a ) = 1.11 â 10 (5b) = 1.11 â 10
-10

T300 cm3 s-1 , T300 exp(-77700/ T ) cm3 s-1 ,
1/2

1/2

(B7) (B8)

-10

which combine to give an equilibrium CO abundance (CO) (OH) = , 0 (C) 1 + (OH) where = exp 77 700 T n . n (H) (B10) (B9)

With the CO abundance determined, the equilibrium H2 O, O and OH abundances are given by equation (B6) and (6a ) = 8.80 â 10 (6b) = 7.44 â 10
-13 1.95 T300 exp(-1429/ T ) cm3 s-1 , 1.57 T300 exp(-9140/ T ) cm3 s-1 ,

(B11) (B12)

-12

which combine to give an equilibrium: (OH) -0 = 8.45T300 (H2 O)
.38

exp(-3940/ T ) cm s ,

3

-1

exp

(B3)

-7711 T

n (H) . n (H2 )

(B13)

where T 300 = T /300 K, and, for simplicity, we assume the H2 is in the ground vibrational state, and

A This paper has been typeset from a TEX/L TEX file prepared by the author.

C

2003 RAS, MNRAS 339, 133­147