Molecular Tracers of Turbulent Shocks

Molecular Tracers of Turbulent Shocks in Giant Molecular Clouds

Abstract

Giant molecular clouds contain supersonic turbulence and simulations of magnetohydrodynamic turbulence show that these supersonic motions decay in roughly a crossing time, which is less than the estimated lifetimes of molecular clouds. Such a situation requires a significant release of energy. We run models of C-type shocks propagating into gas with densities around 10 cm at velocities of a few km s, appropriate for the ambient conditions inside of a molecular cloud, to determine which species and transitions dominate the cooling and radiative energy release associated with shock cooling of turbulent molecular clouds. We find that these shocks dissipate their energy primarily through CO rotational transitions and by compressing pre-existing magnetic fields. We present model spectra for these shocks and by combining these models with estimates for the rate of turbulent energy dissipation, we show that shock emission should dominate over emission from unshocked gas for mid to high rotational transitions (J ) of CO. We also find that the turbulent energy dissipation rate is roughly equivalent to the cosmic-ray heating rate and that the ambipolar diffusion heating rate may be significant, especially in shocked gas.

Subject headings:
ISM: clouds - ISM: molecules - shock waves - stars: Formation - turbulence

1. Introduction

Molecular line observations of giant molecular clouds (GMCs) yield line widths significantly larger than what would be expected from thermal motions alone (e.g., Larson 1981; Solomon et al. 1987). These large, nonthermal line widths are generally interpreted as being due to supersonic turbulence, with Mach numbers on the order of 10 (e.g., Zuckerman & Evans 1974; McKee & Ostriker 2007). Zeeman splitting measurements of magnetic field strengths in molecular clouds show that these supersonic motions are on the order of the Alfvén speed, which suggests that magnetohydrodynamic (MHD) waves may play a significant role in molecular clouds (Crutcher 1999; Crutcher et al. 2010).

Supersonic, hydrodynamic turbulence decays on the order of a free-fall time (e.g., Goldreich & Kwan 1974; Field 1978; Elmegreen 1985) and thus, maintaining the turbulent support of GMCs for their entire lifetimes, estimated to be between 2 and 30 times longer than the free fall timescale (e.g., Mouschovias & Spitzer 1976; Shu 1977; Blitz & Shu 1980; Shu et al. 1987; Williams & McKee 1997; Elmegreen 2000; Hartmann et al. 2001; Mac Low & Klessen 2004), is a significant problem. Based on theoretical calculations (Arons & Max 1975), it was believed that MHD turbulence would decay an order of magnitude slower than hydrodynamic turbulence, thereby preventing the dissipation of turbulent energy in GMCs; however, simulations of MHD turbulence show that MHD turbulence also decays on the order of a free-fall time at the driving scale (Gammie & Ostriker 1996; Mac Low et al. 1998; Stone et al. 1998; Mac Low 1999; Padoan & Nordlund 1999; Ostriker et al. 2001).

In MHD turbulence simulations, turbulent energy is dissipated via numerical viscosity and artificial viscosity in shock fronts. Under the assumption that the dissipated turbulent energy is lost as heat and rapidly radiated away, many MHD simulations are run with isothermal equations of state and thus, these simulations do not explicitly follow where the dissipated turbulent energy goes (e.g., Stone et al. 1998; Smith et al. 2000b). Basu & Murali (2001) made a first attempt to compare the CO luminosities of molecular clouds to what they predicted would be seen from molecular clouds based upon simple energetic arguments. Since then, however, little progress has been made in determining where this turbulent energy goes and whether there are any observational signatures of this dissipated energy.

Shocks increase the temperature and density of the shocked gas, which, in turn, can substantially alter the chemistry of the gas and the emission coming from the gas (e.g., Kaufman & Neufeld 1996b, a), thereby potentially providing a distinct signature and tracer of turbulent energy dissipation via shocks. For typical turbulent velocities and magnetic field strengths of molecular clouds, the magnetic field is capable of transmitting information about the presence of a shock to ions upstream of the shock front. This eliminates discontinuities in gas properties across the shock front and spreads out the thickness of the shock. In turn, this leads to lower temperatures in the shocked gas and prevents molecules from being dissociated. Such a shock is referred to as a continuous, or C-type, shock and is described in more detail in Mullan (1971), Draine (1980), and Draine & McKee (1993).

We run models of C-type, MHD shocks, based upon Kaufman & Neufeld (1996a), propagating into molecular gas with densities around 10 cm at velocities of a few km s, appropriate for the ambient conditions inside of a molecular cloud, to determine which species and transitions dominate the cooling and radiative energy release associated with shock cooling in turbulent molecular clouds. The shock velocities modeled are on the order of the typical turbulent velocity of molecular clouds, which are much lower than the velocities of protostellar outflows that have been the target of previous studies (e.g., Chernoff et al. 1982; Timmermann 1996; Kaufman & Neufeld 1996b, a). These shock models are combined with estimates for the rate of turbulent energy dissipation in molecular clouds to predict the integrated intensities of various shock excited lines coming from an entire molecular cloud and these integrated intensities are then compared to those from photodissociation region (PDR) models based upon Kaufman et al. (1999).

Typical scaling relations of GMCs are presented in Section 2.1 and the turbulent energy dissipation rate of molecular clouds is derived in Section 2.2. The shock and PDR models used in this paper are described in Sections 2.3-2.5 and the results of these models are presented in Section 3. In Section 4, the implications of these results are discussed, and in Section 5, the rate of turbulent dissipation is compared to other known heating mechanisms in molecular clouds. Finally, our findings are summarized in Section 6.

2. Setup

2.1. Scaling Relations of GMCs

Correlations between the size, density, and line-of-sight velocity dispersion of GMCs, as determined through CO observations, are well known and are collectively referred to as Larson’s laws (e.g., Larson 1981; Solomon et al. 1987; Heyer & Brunt 2004). The best fitting scaling relations found by Solomon et al. (1987) are:

(1)
(2)

where R is the effective radius (the radius of a spherical cloud with the same projected surface area as the observed cloud), is the one-dimensional velocity dispersion (which we assume is equal to the observed line-of-sight velocity dispersion), and is the average density. These relations suggest that a cloud with a mean density of 10 cm has a radius of approximately 2 pc, a mass of 2000 M, a total molecular hydrogen column density of cm (corresponding to a visual extinction of 12 through the entire cloud), and a one-dimensional velocity dispersion of about 1 km s.

The size-velocity relationship is fairly well established, although there is some evidence that the velocity dispersion of a molecular cloud may also depend upon the column density of that cloud (Heyer et al. 2009). The validity of the size-density relation, however, is much less certain, as the observed relationship may be only due to the limited dynamical range of current observations (Ballesteros-Paredes & Mac Low 2002). For the shock models used in this paper, Larson’s laws are only used to confirm that the simulated parameter range roughly corresponds to the properties of observed molecular clouds. Neither part of Larson’s laws is used to calculate the integrated intensity of the shock emission.

Solomon et al. (1987) also found a correlation between the velocity dispersion and CO J = 1 0 total luminosity, L, of molecular clouds. In units of K km s pc, Solomon et al. (1987) found that the CO 1 0 luminosity is:

(3)

The relationship between the magnetic field strength in a molecular cloud, B, and the number density of hydrogen nuclei, , is often expressed in the form

(4)

where b and k are the fitting parameters. A value of k = 0.5 corresponds to a constant magnetic energy density (McKee & Ostriker 2007) and is expected from ambipolar diffusion collapse models (Fiedler & Mouschovias 1993). A k = 0.5 relation is also expected if the turbulent velocity in a cloud is always roughly the Alfvén speed (e.g., Crutcher 1999). A value of k = 2/3, however, is predicted if magnetic fields are unimportant and a molecular cloud is able to maintain a roughly spherical shape during its collapse (Mestel 1966; Crutcher 1999; Crutcher et al. 2010).

Crutcher (1999) compiled Zeeman splitting observations and found b = 0.95 and k = 0.5 (McKee & Ostriker 2007). The MHD simulations of Padoan & Nordlund (1999) exhibit a k = 0.4 relation and the relation , is commonly adopted (e.g., Draine et al. 1983; Kaufman & Neufeld 1996b, a). Recently, Crutcher et al. (2010) have examined all of the available Zeeman splitting observations, including those used by Crutcher (1999), and found that the best fit for the maximum observed line-of-sight magnetic field strength comes from the relation:

(5)

For densities greater than 300 cm, the above relation corresponds to k = 2/3 and b = 0.25. Crutcher et al. (2010) also note that their data are consistent with having line-of-sight magnetic field strengths down to essentially zero.

The maximum magnetic field strengths that were fit by Crutcher et al. (2010) most likely correspond to cases where the magnetic field is highly aligned with the line-of-sight, such that the full magnetic field strength is measured. The average magnetic field strength along any random direction will thus be only half of the strength given by the above relation. For a cloud with an H density of 10 cm, the Crutcher et al. (2010) relation therefore predicts that the average magnetic field strength along any random direction is 17 G. Intrinsic scatter in the magnetic field strength between different clouds will also likely further reduce the average magnetic field strength.

2.2. Turbulent Energy Dissipation Rate

The turbulent energy density of a molecular cloud is approximately

(6)

Following the discussion in Basu & Murali (2001), the mean turbulent energy dissipation rate per volume can be written as , where is the dissipation timescale. We define the flow crossing time of the cloud as and introduce the ratio of the dissipation time to the flow crossing time as a new parameter: . The turbulent dissipation rate per volume is thus

(7)

As shown in Basu & Murali (2001), rather than writing the turbulent energy dissipation rate in terms of , the dissipation rate can be expressed in terms of the driving scale of the turbulence, :

(8)

where is a dimensionless parameter that is a function of the density, velocity dispersion, and driving wavelength. Comparing Equation (7) to Equation (8) gives the relation

(9)

Periodic box simulations of MHD turbulence have found that for a variety of initial conditions, has a value between 0.5 and 4 (Gammie & Ostriker 1996; Stone et al. 1998; Mac Low 1999; Ostriker et al. 2001). Unfortunately, there is no clear consensus on what scale turbulence is driven on. Protostellar outflows, which drive turbulence on small scales, appear to have enough energy to drive turbulence in active star forming regions (Quillen et al. 2005; Curtis et al. 2010; Arce et al. 2010). It is, however, unclear whether outflows are capable of driving turbulence across an entire molecular complex (Banerjee et al. 2007; Arce et al. 2010). Studies of density and velocity structure in molecular clouds find that the observed structures are only consistent with driving at size scales at, or above, the size of the cloud (e.g., Ossenkopf & Mac Low 2002; Brunt 2003; Heyer & Brunt 2004; Brunt et al. 2009; Padoan et al. 2009). Supersonic turbulence has also been observed in the Polaris Flare, which is devoid of any protostars (André et al. 2010).

For the remainder of this paper, a value of one will be adopted, for which the turbulent dissipation timescale is equal to the flow crossing time of the cloud. Via Equation (9) and the numerical factors for , this corresponds to a turbulent driving scale on the order of the size of a molecular cloud.

Equation (7) is a general result that can be applied to any cloud, given that the characteristic radius, density, and velocity dispersion are known. For this paper, the simplifying assumption that clouds are spherical will be made, such that the total turbulent energy dissipation rate is

(10)
(11)

If all of the dissipated turbulent energy is radiated away, the corresponding total integrated intensity, , is

(12)
(13)

This integrated intensity is independent of the size of the molecular cloud. For this paper, a mean mass per particle of g, or about 2.77 amu, is adopted, and with this value, the above equation becomes

(15)

2.3. Shock Code

To determine which species and transitions dominate the cooling and radiative energy release associated with shock cooling of turbulent molecular clouds, we run models of C-type shocks based upon the models of Kaufman & Neufeld (1996a) with initial conditions corresponding to that expected for roughly one parsec sized molecular clouds.

The code used first calculates the temperature, density, chemical abundance, and velocity profiles of a C-type shock. To do so, it calculates the cooling rates for rotational and vibrational transitions of HO, H, and CO (Neufeld & Kaufman 1993); collisions between the neutral gas and cooler dust grains (Hollenbach & McKee 1989); and H dissociative cooling (Lepp & Shull 1983; Hollenbach & McKee 1989). The freezing out of molecules on dust grains is, however, not calculated by the code as the freezeout timescale is expected to be much longer than the shock cooling timescale.

Once the shock structure is determined, the code then calculates the integrated intensities of each molecular transition of interest by solving the partial differential equations for the line emission at each point and then integrating the emission over the entire shock profile. This extra step of determining individual line strengths, rather than just determining an overall cooling rate for a particular molecule, is only used for CO. In this later step, the code only includes rotational line emission for gas down to 10 K, unlike in Kaufman & Neufeld (1996a), where emission is only included from gas above 50K. No such temperature limitation is used when calculating the overall cooling rates for each molecule in the first half of the code. For a more detailed description of how this code works, please see Kaufman & Neufeld (1996a).

By equating the total kinetic energy dissipated by shocks to the turbulent energy dissipation rate of a molecular cloud, given by Equation (11), we scale our shock models to predict the expected integrated intensities from each CO rotational line. That is, the integrated intensity of each line is set to the appropriate fraction of the integrated intensity given by Equation (2.2). It is assumed that the shock emission is coming from a region larger than the size of the beam and each line is scaled equally under the assumption that all of the lines are optically thin. The effect of the lower lying lines being optically thick is discussed further in Section 4.6.

2.4. Shock Code Parameters

For each shock model, the same, roughly solar, chemical composition as used in Kaufman & Neufeld (1996b, a) is used. In particular, the initial CO number density is set to be times that of the H nuclei number density and the initial HO abundance is set to 10. As shown in Section 2.1, a one parsec sized molecular cloud is expected to have a density of approximately 10 cm. Thus, the initial H density is set to be either 10, 10, or 10 cm in these models.

If the velocity distribution of gas particles in a molecular cloud is Gaussian in every direction, with a one-dimensional velocity dispersion of , then the distribution of relative velocities between two gas particles in the cloud will also be Gaussian in every direction with a one-dimensional velocity dispersion of . Since the energy dissipation rate of a shock scales with the third power of the shock speed, the mean speed at which energy is dissipated is the cube root of the mean cubed velocity difference between two gas particles, , which is roughly 2.4. The shock velocity at which the peak energy dissipation rate occurs is slightly higher, approximately 3.2. Thus, the characteristic shock velocity in a molecular cloud with a one-dimensional velocity dispersion of 1 km s, consistent with the size-velocity relation for a radius of 1 pc, is on the order of 2-3 km s. For the remainder of this paper, we assume that the one-dimensional velocity dispersion is a factor of 3.2 smaller than the shock velocity. The larger conversion factor of the two mentioned above is chosen so that the corresponding velocity dispersions, and thus the shock-integrated intensities calculated in Section 3.1, are smaller.

Models with shock velocities of 2 and 3 km s are computed. For a temperature of 10 K, these velocities correspond to Mach numbers of 12 and 17, respectively. While these velocities are appropriate for turbulent motions in a molecular cloud, they are much lower than the typical velocities of protostellar outflows and winds. Such higher velocity flows have been modeled extensively in the past and give rise to significantly higher post shock temperatures (e.g., Kaufman & Neufeld 1996b, a).

The strength of the magnetic field parallel to the shock front is initialized using the parameterizations k = 0.5 and b = 0.1 or 0.3, where b and k are as defined in Equation (4). The component of the magnetic field perpendicular to the shock front is always set to zero, as this component has no effect on the shock structure in our steady state, plane parallel models. Thus, the initial magnetic field strength ranges from 3 to 24 in the different models. For a weaker field parallel to the shock front, the shock thickness is smaller and the energy released in line radiation is relatively larger. This is why magnetic field strengths that are slightly lower than, although still generally consistent with, the average line-of-sight magnetic field strength given by the scaling relation of Crutcher et al. (2010) have been chosen.

The Alfvén speed is

(16)

For our shock models, the Alfvénic Mach number, given by , ranges from 4 to 16.

Twelve shock models are run, one for each combination of initial density, shock velocity, and magnetic field b parameter. Table 1 gives the shock velocity, magnetic field strength, Mach number, and Alfvénic Mach number for each model. A naming convention of nWXvYbZ is adopted, where W.X is the logarithm of the initial H number density in cm, Y is the shock velocity in km s, and Z is the magnetic b parameter.

Model log(n) v b B Mach
(cm) (km s) (G)
(1) (2) (3) (4) (5) (6) (7)
n25v2b1 2.5 2 0.1 3 12 11
n25v3b1 2.5 3 0.1 3 17 16
n25v2b3 2.5 2 0.3 8 12 4
n25v3b3 2.5 3 0.3 8 17 5
n30v2b1 3 2 0.1 4 12 11
n30v3b1 3 3 0.1 4 17 16
n30v2b3 3 2 0.3 13 12 4
n30v3b3 3 3 0.3 13 17 5
n35v2b1 3.5 2 0.1 8 12 11
n35v3b1 3.5 3 0.1 8 17 16
n35v2b3 3.5 2 0.3 24 12 4
n35v3b3 3.5 3 0.3 24 17 5

Note. – Column 1 shows the model name, while Columns 2 and 3 show the logarithm of the initial density and shock velocity of each model respectively. Column 4 represents the magnetic b parameter, as defined in Equation (4), and Column 5 shows the resulting initial magnetic field strength. Columns 6 and 7 show the Mach number and Alfvénic Mach number of the models, respectively.

Table 1Shock Model Properties

2.5. PDR Model

The shocked gas in a molecular cloud is not the only source of molecular line emission. The cool, well-shielded gas and the warmer gas in the PDR at the cloud’s surface, which is exposed to the interstellar radiation field (ISRF), will also contribute emission. To model this emission from unshocked gas, PDR models based on Kaufman et al. (1999) are used. As suggested by Kaufman et al. (1999), these plane parallel models are adjusted for spherical geometry by using the equation

(17)

where r is the radial distance from the center of the cloud and j(N) is the emissivity at a column N from the surface of the cloud. Kaufman et al. (1999) estimate that this procedure produces results that are within a factor of 1.5 from intrinsically spherical PDR models. Furthermore, for any optically thin line, the resulting integrated intensity is doubled to account for photons originally emitted radially inward.

The PDR models used have ISRFs of 3 Habing, where the average far ultraviolet ISRF in free space is 1.7 Habing or erg cm s (Tielens 2005), and microturbulent Doppler line widths of 1.5 km s, similar to the velocity dispersions of the shock models. It is assumed that the PDR emission fills the beam. These PDR models do not take into account the freezing out of CO onto dust grains.

A density of 10 H cm is used for all of the comparison PDR models, which is comparable to the median initial density in the shock models. We believe that this is an appropriate comparison density for the 10 cm shock models because a density gradient should be present within realistic molecular clouds, with the density decreasing toward the periphery of the cloud, in order for the clouds to remain in pressure equilibrium. Thus, it is expected that the warm outer layers of molecular clouds, from which most of the PDR emission comes from, are at lower densities than the bulk of the cool, CO-rich gas in the interior of the clouds from which most of the shock emission originates. In PDR models with densities below 10 H cm, CO only forms in the gas phase once the gas has cooled significantly, such that there is almost no emission in the mid to high J lines of CO. To be more conservative in our findings of when shock emission is stronger than PDR emission, we prefer using PDR models with densities of 10 H cm for comparison with the 10 and 10 cm shock models, even though these PDR models may over predict the CO emission from 10 cm gas. The results should be roughly consistent for the 10 cm shock models.

For the comparison PDR models, the size of the molecular cloud must also be known in order to know at what A to cut the model. For each shock model, the Solomon et al. (1987) size-velocity relation is used to determine an appropriate, typical size for a cloud and then the CO column density of that cloud is determined under the assumptions that the cloud is spherical and has a typical CO abundance of (Glover et al. 2010) throughout the entire cloud. This CO column density is then used to determine the appropriate depth of the comparison PDR model such that the PDR model has the same CO column density. The depths of the PDR models are chosen based upon a CO column density, rather than upon a hydrogen nuclei column density, because the initial chemical abundances used for the shock models are consistent with gas in which CO has already formed in the gas phase and because CO cooling dominates the energy budget of the shock models, as described later in Section 3.1.

Two shock models, n35v3b1 and n35v3b3, require CO column densities larger than the total CO column present at the maximum depth of the 10 cm PDR model. The contribution from the most deeply embedded layers of this PDR model to the total emergent flux does, however, drop to negligible values for all of the CO lines. This is because the lower lines are optically thick and the gas temperature at high A is too low for any significant emission in the higher lines. Thus, we use the full extent of the 10 cm PDR model as the comparison model for these two shock models and we do not believe that this failure to exactly match the CO column densities of these two models significantly affects our results.

2.6. Empirical CO 1 0 Luminosities

As described in Section 2.1, Solomon et al. (1987) found an empirical correlation between the velocity dispersion and CO J = 1 0 total luminosity of molecular clouds. For the velocity dispersion corresponding to each of the shock models, the 1 0 integrated intensity expected from this Solomon et al. (1987) relation is calculated using a cloud radius from the Solomon et al. (1987) size-velocity relation and the assumption that molecular clouds are spherical. While we have used the size-velocity relationship in determining the comparison PDR spectra and in calculating the empirically expected 1 0 integrated intensity, we re-emphasize that this relationship is not used at any time in the calculation of the expected shock spectra.

3. Results

3.1. Shock Line Emission

In the upper panels of Figure 1, the neutral velocity, ion velocity, and density profiles of models n30v2b1 and n30v3b1 are shown. As typical in MHD, C-type shocks, the density and velocity profiles show no sharp discontinuities and the ion velocity decreases before the neutral velocity does. Due to mass conservation, the density and neutral velocity are inversely related to each other, and thus, the maximum density reached in model n30v3b1 is larger than that reached in model n30v2b1. The magnetic field strength and ion velocity are similarly inversely correlated.

The lower panels of Figure 1 show the temperature profiles of models n30v2b1 and n30v3b1 as well as the cooling profiles due to CO, H, and HO lines and gas-grain coupling. During the initial stages of a shock, where the gas temperature is still increasing, the rate of CO cooling increases in close tandem to the increase in temperature. After the gas has reached its peak temperature, the CO cooling rate remains higher than it was at the same temperature earlier in the shock. This is due to the higher gas densities in these later regions of the shock that allow for more efficient population of higher J CO states and thus, more effective CO cooling. The CO cooling rate is the least temperature sensitive of any of the plotted cooling terms, whereas the H cooling rate shows a strong dependence on temperature. The H cooling rate is strongly peaked around the temperature peak and shows the most significant change between the two shock models. While the gas-grain cooling rate is temperature dependent, it is also clearly larger at higher densities as the gas-grain cooling curves are skewed toward the higher density sides of the shocks. High frequency noise, which is likely numerical in nature, appears toward the end of both models in the CO cooling rate and thus, a boxcar smoothing algorithm has been applied to the CO cooling rates for distances larger than 0.06 pc in model n30v2b1 and for distances larger than 0.04 pc in model n30v3b1. This noise does not significantly affect our results as it only occurs over very limited spatial scales and occurs only when the cooling rates have decreased significantly, such that the noise does not significantly affect the total cooling rate. Furthermore, this noise only occurs when the temperature has dropped below 10 K, at which point the lack of cosmic-ray heating in our models becomes important (see Section 4.1). The other shock models are qualitatively similar to the ones shown in Figure 1 and thus, are not shown. The temperatures, densities, and timescales of the shocks are too low for any significant chemical changes to occur within the gas in any of the models and thus, the small changes in chemical abundance across the shock models are also not shown.

Figure 1.— Various profiles of models n30v2b1 and n30v3b1. The top row shows density, neutral velocity, and ion velocity profiles as the solid (black), dotted (blue), and dashed (red) lines, respectively. The velocity axis is given on the left-hand border while the density axis is given on the right-hand border. The bottom row shows temperature profiles as the solid (black) lines and cooling profiles due to CO, H, gas-grain interactions, and HO as the dotted (blue), dashed (red), dash-dotted (green), and dash-triple-dotted (yellow) lines, respectively. The cooling rate axis is given on the left border and the temperature axis is given on the right border. The CO cooling profiles have been boxcar smoothed beyond a distance of 0.06 pc in model n30v2b1 and past 0.04 pc in model n30v3b1 due to the presence of high frequency noise. This noise is likely numerical in nature and should not significantly affect our results (see the text in Section 3.1). The left-hand column shows profiles of the n30v2b1 model and the right-hand column shows profiles of the n30v3b1 model. The x-axes of all four boxes are the same and the y-axis scaling is the same for both models. Please see the online version for a color version of this figure.

The dominant molecular coolant in all of these slow shock models is CO, with 40%-80% of the dissipated energy going into CO rotational lines. A significant fraction, 15%-60%, of the dissipated energy is not radiated away, but rather, goes toward compressing the magnetic field. In the models with the weakest shocks, those with b = 0.3 and a shock velocity of 2 km s, the conversion of kinetic energy into magnetic energy is the most significant mechanism for dissipating kinetic energy. Molecular hydrogen rotational lines are the second most effective molecular coolant, but dissipate less than 1% of the shock energy in all but the models with b = 0.1 and a shock velocity of 3 km s, which are the models with the strongest shocks. In these stronger shock models, H lines account for between 7% and 21% of the energy dissipated, with H cooling being more important at lower densities. All other cooling mechanisms are very minor in these shock models. A summary of where the energy goes in each model is given in Table 2. It should be noted that the sums of the cooling functions do not exactly equal the total kinetic energy dissipation rates. The total cooling rate, however, is never more than 5% discrepant from the kinetic energy dissipation rate. We believe that this discrepancy arises from difficulties in extending our shock cooling functions to low temperature but do not believe that this small discrepancy significantly affects our results.

Model
(1) (2) (3) (4) (5) (6)
n25v2b1 76 23 0.8 2 0.02
n25v3b1 61 17 21 2 0.02
n25v2b3 42 56 0.1 0.1 0.01
n25v3b3 54 44 0.9 1 0.01
n30v2b1 76 26 0.2 2 0.02
n30v3b1 68 18 14 3 0.03
n30v2b3 43 54 0.1 0.2 0.01
n30v3b3 55 41 0.3 2 0.02
n35v2b1 74 25 0.1 3 0.04
n35v3b1 72 19 7 4 0.04
n35v2b3 41 53 0.1 1 0.03
n35v3b3 53 42 0.1 3 0.03

Note. – Column 1 shows the model names. Columns 2-6, respectively, list the percentage of the kinetic energy of the shock that is dissipated via CO rotational lines, increasing the magnetic field strength, H lines, gas-grain collisions, and HO lines.

Table 2Sources of Energy Dissipation in the Shock Models

The 12 modeled shocks are relatively weak shocks and produce density enhancements of at most a factor of 20, of the order of the Mach number. Such compressions are still smaller than the density contrast between the ambient material in molecular clouds and dense cores, which have densities of 10 cm or above (e.g., Di Francesco et al. 2007). The maximum temperature in the shocked gas varies significantly, with the stronger shock models achieving maximum temperatures of approximately 150 K and the weaker shock models not even warming up to 20 K. The maximum density and temperature reached in each model is given in Table 3.

The cooling length of each shock is taken to be the full width at quarter maximum of the total cooling function profile. The cooling lengths range from 0.01 pc to 0.35 pc, with the high magnetic field strength and low-density models having the largest cooling lengths. The corresponding cooling timescales range from years to years, with the longer cooling timescales corresponding to larger cooling lengths.

The volume filling factor of shocked gas in a molecular cloud, , can be calculated from

(18)

where is the turbulent energy dissipation rate per volume (given by Equation (7)), is the cooling length of the shock, and is the kinetic energy dissipated per shock front area. Both and are calculated by the shock code but the cloud radius is required to calculate . For this filling factor calculation, we use the relatively well-established size-velocity relation of Larson’s laws, Equation (1), to determine the appropriate cloud radius for each shock model. We reiterate, however, that the integrated intensities presented in Figures 2-4 and in Table 4 are derived independently from any part of Larson’s laws. The volume filling factor of shocked gas is always between 0.02% and 0.5% of the cloud volume, except for the three weakest shock models, where the volume filling factor becomes as large as 2%. The cooling times, cooling lengths, and filling factors for all 12 models are given in Table 3.

Model log(n) ff
(cm) (K) ( years) (pc) (%)
(1) (2) (3) (4) (5) (6) (7) (8)
n25v2b1 3.7 56 6.4 0.06 0.38 14 54
n25v3b1 3.8 145 3.0 0.04 0.11 22 104
n25v2b3 3.1 11 28.5 0.35 2.15 4 10
n25v3b3 3.3 60 9.8 0.17 0.46 7 19
n30v2b1 4.1 54 2.6 0.03 0.16 14 57
n30v3b1 4.3 154 1.3 0.02 0.05 22 106
n30v2b3 3.6 13 11.8 0.14 0.87 4 9
n30v3b3 3.8 60 3.9 0.07 0.18 7 18
n35v2b1 4.6 53 1.1 0.01 0.07 14 55
n35v3b1 4.7 157 0.5 0.01 0.02 21 108
n35v2b3 4.1 17 4.9 0.06 0.36 4 9
n35v3b3 4.3 61 1.7 0.03 0.08 7 18

Note. – Column 1 shows the model names while Columns 2 and 3 show the logarithm of the maximum density reached and the maximum temperature reached, respectively. Column 4 shows the cooling time of the shocked gas and Column 5 shows the cooling length. The volume filling factor of shocked gas for a cloud compatible with the size-velocity relationship of molecular clouds is given in Column 6. Columns 7 and 8 show the factors by which the magnetic field strength and the ambipolar diffusion volume heating rate increase between the initial and shocked gas.

Table 3Shock Model Properties

Figure 2 shows the integrated intensities of CO rotational transitions as calculated from the shock models with densities of 10 cm. The CO spectra from the corresponding comparison PDR models, as well as the estimates for the J = 1 0 integrated intensities from the Solomon et al. (1987) scaling relation, are also shown in Figure 2. Figures 3 and 4 show the shock spectra from the models with densities of 10 cm and 10 cm respectively, as well as the corresponding PDR spectra and Solomon et al. (1987) J = 1 0 integrated intensities. All integrated intensities are given in units of erg s cm arcsec. The y-axis ranges in Figures 2, 3, and 4 are identical to facilitate comparisons between the different shock models. As explained in further detail in Section 4.1, CO spectra have not been calculated for any shock models with b = 0.3 and a shock velocity of 2 km s.

Figure 2.— Integrated intensities of various CO rotational transitions for shock models with densities of 10 cm in units of erg s cm arcsec. The shock velocity and magnetic field b parameter used for each shock model are given on the top and left of the grid, respectively, while the model name is given in the top left hand corner of each box. The green (lightest) lines show the Solomon et al. (1987) CO line strengths, the blue (darkest) lines show the CO shock spectra, and the red (medium) lines show the CO PDR spectra. The ten lowest rotational transitions of CO are labeled in the lower right grid panel. Note how the shock spectra dominate over the PDR spectra for high J transitions.
Figure 3.— Integrated intensities of various CO rotational transitions for shock models with densities of 10 cm in units of erg s cm arcsec. The shock velocity and magnetic field b parameter used for each shock model are given on the top and left of the grid, respectively, while the model name is given in the top left hand corner of each box. The green (lightest) lines show the Solomon et al. (1987) CO line strengths, the blue (darkest) lines show the CO shock spectra, and the red (medium) lines show the CO PDR spectra. The ten lowest rotational transitions of CO are labeled in the lower right grid panel. Note how the shock spectra dominate over the PDR spectra for high J transitions.
Figure 4.— Integrated intensities of various CO rotational transitions for shock models with densities of 10 cm in units of erg s cm arcsec. The shock velocity and magnetic field b parameter used for each shock model are given on the top and left of the grid, respectively, while the model name is given in the top left hand corner of each box. The green (lightest) lines show the Solomon et al. (1987) CO line strengths, the blue (darkest) lines show the CO shock spectra, and the red (medium) lines show the CO PDR spectra. The ten lowest rotational transitions of CO are labeled in the lower right grid panel. Note how the shock spectra dominate over the PDR spectra for high J transitions.

In all of the cases considered, the PDR emission is significantly stronger, by approximately an order of magnitude, than the shock emission in the three lowest CO rotational transitions. The Solomon et al. (1987) CO J = 1 0 integrated intensity is also larger than the predicted shock J = 1 0 integrated intensity in all of the models. On the other hand, in all of the models, the shock emission is stronger for all of the high J transitions. In the strongest, high-density shock models, the shock-integrated intensity becomes larger than the PDR-integrated intensity in the J = 5 4 line. The PDR emission in the higher J lines drops much more rapidly than the shock emission and the shock emission is an order of magnitude stronger in the J = 7 6 line in most models.

While the shock code does not calculate the individual strengths of the different H rotational transitions, the lowest rotational levels should be in local thermodynamic equilibrium (LTE) at the modeled densities. Thus, it is assumed that all of the H line ratios are given by their LTE values. At 50 K, the ratio of S(1) to S(0) is only , but at 150 K, this ratio increases to 2.3. The S(2) and higher lines are negligible in comparison to the S(1) and S(0) lines at these temperatures. For the three models with significant H emission, the expected integrated intensities of the S(1) and S(0) lines are given in Table 4.

While the PDR models used do not include the expected line strengths for H rotational transitions, Kaufman et al. (2006) found that the S(1) and S(0) H lines have integrated intensities of erg s cm arcsecond and erg s cm arcsecond, respectively, in a PDR with a density of cm and an ISRF of 3 Habing. Therefore, the PDR H-integrated intensities are expected to be comparable to or slightly lower than the shock H-integrated intensities.

Model S(0) S(1)
(10 erg s (10 erg s
cm arcsec) cm arcsec)
(1) (2) (3)
n25v3b1 1.5 3.5
n30v3b1 3.0 6.9
n35v3b1 5.0 11

Note. – Column 1 shows the model names while Columns 2 and 3 show the integrated intensities of the S(0) and S(1) lines, respectively, for the three shock models with significant H cooling.

Table 4Predicted H-Integrated Intensities

4. Discussion

4.1. CO Lines

As mentioned in Section 3.1, no CO spectra are calculated for the models with b = 0.3 and a shock velocity of 2 km s. This is because the maximum post shock temperatures in these three models, 11, 13, and 17 K, are on the order of the temperature expected from cosmic-ray heating alone (e.g., Goldsmith 2001; Pan & Padoan 2009). The shock code used does not contain any additional heating sources, such as cosmic-ray or photoelectric effect heating, and thus, underestimates the temperatures that would be achieved in these weak shocks. The combination of shock heating and cosmic-ray heating must produce temperatures higher than that produced by cosmic-ray heating alone. Note that in Figure 1, the gas temperature is initially less than 0.1 K and falls below 10 K by the end of the model due to the lack of these extra heating terms. As such, we do not consider the models with b = 0.3 and a shock velocity of 2 km s to be valid. The lack of cosmic-ray heating does not significantly affect the other models because the peak temperatures are much larger than 10 K in these models, indicating that shock heating is significantly stronger than what cosmic-ray heating would be.

Figures 2-4 show that almost all of the emission from a molecular cloud in the lowest three rotational transitions of CO comes from unshocked gas. These figures also show, however, that most of the emission coming from the mid to high J CO lines (lines at or above J = 7 6) comes from shocked gas. Not only should the integrated intensities of these lines be higher than predicted from PDR models, but the excitation temperatures derived from the line ratios of these lines should also be higher than a PDR model would predict. Thus, these mid to high J CO lines should serve as observational diagnostics of turbulent energy dissipating via shocks. Since CO accounts for the majority of the cooling via shocks, these mid to high J transitions of CO should be the best tracers of where the majority of the energy goes in turbulent shocks in molecular clouds.

Shock emission dominates at higher J transitions because the CO in the shocked gas is warmer than the majority of the CO in the rest of the cloud. While the outer layers of a molecular cloud can be quite warm, due to the incident ISRF, there is little CO in these warm outer regions. The majority of the CO flux in the PDR models comes from gas that is below 20 K. Since our low-velocity shocks are C shocks, CO survives the shock and radiates from gas at temperatures in excess of 50 K.

The Solomon et al. (1987) CO J = 1 0 integrated intensity is larger than the predicted shock-integrated intensity in all of the models, which confirms that shock emission is not the major source of emission in the 1 0 transition. While the comparison PDR models for the high-density shock models over predict the 1 0 integrated intensity, compared to the Solomon et al. (1987)-integrated intensity, the comparison PDR models for the 10 cm shock models have J = 1 0 integrated intensities in remarkably good agreement with the Solomon et al. (1987)-integrated intensity. This slight discrepancy between some of the PDR models and the Solomon et al. (1987) relation is likely due to the characteristic line width, depth, and density of the clouds observed by Solomon et al. (1987) being slightly different than the values used in the PDR models.

4.2. Variation across Parameter Space

The maximum temperature reached in the shocked gas increases with increasing Alfvénic Mach number and Mach number of the shock. All of the models with a shock velocity of 3 km s, and thus a Mach number of 17, have higher maximum temperatures than all of the models with a shock velocity of 2 km s, corresponding to a Mach number of 12. For models with the same Mach number, models with higher Alfvénic Mach numbers, those with lower magnetic b parameters, have higher maximum temperatures. A larger Alfvénic Mach number alone, however, does not necessarily imply a larger maximum temperature, as the models with b = 1 and a shock velocity of 2 km s have roughly the same maximum temperature as the models with b = 3 and a shock velocity of 3 km s despite having almost twice as large Alfvénic Mach numbers. The maximum temperature reached significantly affects the CO shock profile, as higher temperatures excite higher rotational states, which leads to considerably more emission in higher J transitions. As described further in Section 4.4, the effectiveness of H cooling is also highly sensitive to the maximum temperature reached in the gas.

The fraction of energy going toward compressing the magnetic field is primarily dependent upon the Alfvénic Mach number of the shock, with more energy going into the magnetic field in the models with lower Alfvénic Mach numbers. Larger Alfvénic Mach numbers also produce smaller cooling lengths, cooling times, and filling factors.

The turbulent energy density of a molecular cloud is dependent upon both the density and turbulent velocity dispersion of the cloud. Thus, the models with higher densities and higher shock velocities have larger integrated intensity scaling factors (see Equation (2.2)) and, consequently, higher integrated intensities for all CO transitions.

The critical densities for all the CO lines above the J = 3 2 line are greater than 10 cm (Schöier et al. 2005). Therefore, larger initial densities make it easier for higher rotational states to be populated and, as such, more emission comes out in higher lying lines in the higher density models. This effect is clearly seen in the shock models with b = 0.1 and a shock velocity of 2 km s, as the line with the largest integrated intensity shifts from the 3 2 line to the 4 3 and then finally to the 5 4 line as the density increases from 10 cm to 10 cm and then to 10 cm. The peak of the PDR spectra remains at the 3 2 transition in all of the PDR comparison models because the same density is used for all of the PDR comparison models. This difference between the shock and PDR models causes the transition at which shock emission becomes larger than PDR emission to move to slightly lower transitions as the density of the shock model increases.

As described above, changes to the magnetic field strength, shock velocity, or initial gas temperature can significantly alter the shape and the scaling of the shock CO spectrum. None of these changes, however, alter the key result that shock emission dominates at mid to high J transitions, particularly from J = 7 6 and up.

The fraction of energy dissipated via CO cooling is not strongly correlated with any shock property. Only in the strongest shock models does the fraction of energy dissipated via CO weakly depend upon the initial density of the gas. For these strong shock models, the fraction of energy that is emitted via CO rotation lines increases by approximately 10% as the initial density increases from 10 cm to 10 cm.

An increase in the initial density of the gas weakly increases the effectiveness of gas-grain cooling, but this cooling term accounts for at most 4% of the shock cooling in any of the models. Much higher densities, densities closer to 10 cm, are required before gas-grain coupling becomes reasonably efficient. Going from a density of 10 cm to 10 cm also decreases the cooling lengths, cooling times, and filling factors of all of the models by approximately a factor of six.

4.3. Magnetic Field Compression

In our low velocity shock models, the energy that goes into compressing the magnetic field is of the same order of magnitude as the energy radiated away in CO rotational lines. It is, however, unclear where this injected magnetic energy would go. A local increase in magnetic field strength could further drive MHD waves, which would subsequently shock. In this case, more shocks would be required in order for all of the cloud’s turbulent energy to be dissipated by CO cooling and our predicted line integrated intensities would have to be increased by a factor of two.

Alternatively, this magnetic energy may slowly leak out of the cloud via magnetic coupling with the external medium, as described by Elmegreen (1985) and seen in simulations by Eng (2002). This magnetic energy may also be dissipated on small scales via a process such as ambipolar diffusion.

4.4. H Lines

Molecular hydrogen does not have a permanent dipole moment and thus, radiates through weak quadrupole transitions (J = 2). Furthermore, since hydrogen is so light, the rotational energy levels of molecular hydrogen are relatively widely spaced. These two effects make H rotational emission highly temperature sensitive and temperatures in excess of 100 K are required for significant emission. This temperature sensitivity can be seen in the shock models, as H emission is essentially negligible in all but the strongest shock models, which are the only models where the maximum temperature exceeds 100 K. In these models, the S(1) and S(0) H lines have comparable integrated intensities to the CO J = 7 6 line, although the H lines are relatively stronger in the lower density models.

The lack of H rotational emission from cool gas means that, aside from shocked gas, the only significant source of H emission in a molecular cloud is the thin outer edge of the cloud’s PDR where the temperature is high and H is not rapidly photodissociated. Thus, H rotational emission could be a very useful tracer of shocked gas in a molecular cloud. In particular, while the predicted S(0) shock-integrated intensities are on the order of the S(0) PDR-integrated intensity, the S(1) shock-integrated intensities range from being five times to over fifteen times larger than the S(1) PDR-integrated intensity. The ratio of S(1)/S(0) at the 150 K maximum temperature of the strongest shocks, approximately 2.3, is also significantly different from the ratio of these lines in the Kaufman et al. (2006) comparison PDR models, roughly 0.3.

H cooling may be even more significant in gas that is lacking in gas phase CO, as this shocked gas is likely to reach higher temperatures with the effectiveness of CO cooling reduced. This increase in H shock emission in CO sparse gas, such as the “dark gas” in the periphery of a molecular cloud (e.g., Wolfire et al. 2010), may naturally produce a limb brightening effect for H rotational emission in molecular clouds, as observed in Taurus by Goldsmith et al. (2010).

In all of our shock models, cooling from vibrational transitions of H is negligible because temperatures on the order of a few thousand Kelvin are required to excite higher energy vibrational states (Kaufman & Neufeld 1996a).

4.5. Other Shock Tracers

Water is known to be an effective coolant in high-velocity shocks (Kaufman & Neufeld 1996a) but water cooling is negligible in all of our shock models. This is because water is only formed in gas with a temperature of a few hundred Kelvin (Elitzur & de Jong 1978; Elitzur & Watson 1978) and is only efficiently liberated from dust grains, due to sputtering, in 15 km s or faster shocks (Draine et al. 1983). Low-velocity shocks are also very ineffective at heating dust grains, meaning that thermal sublimation of water off of dust grains is completely negligible in our low-velocity shocks (Draine et al. 1983).

At the low densities of our shock models, the interaction timescale of gas with dust grains is long compared to the cooling time such that only a few percent of the energy is ever liberated via gas-grain interactions in the models. Densities closer to 10 cm are needed before gas-grain coupling becomes effective.

While CO line radiation and the compression of magnetic fields are the dominant coolants in our low-velocity shocks, other molecular lines, which have not been included in our shock models, may still be valuable tracers of shocked gas. In particular, molecular lines that are sensitive to increased temperature could provide shock tracers in different wavelength regimes. Molecular transitions that are sensitive to density may also be useful shock tracers, but are less likely to be as useful as temperature sensitive lines since the maximum densities reached in the shocked gas are less than the gas density of prestellar cores.

4.6. Additional Caveats

In scaling the shock models to predict the total strength of the shock emission from an entire molecular cloud, it was assumed that all of the turbulent energy of the cloud is dissipated at one particular shock strength. In reality, energy will be dissipated through a variety of different strength shocks, either due to different shock velocities or different strengths of the magnetic field parallel to the shock front. Furthermore, Smith et al. (2000a) show that high-velocity shocks dominate energy dissipation in driven turbulence while Smith et al. (2000b) show that if the turbulence is decaying, low-velocity shocks dissipate most of the energy. If energy is dissipated through lower velocity shocks, then our calculated line-integrated intensities will overestimate the actual emission from molecular clouds in lines with higher excitation temperatures. If turbulence is driven at scales much smaller than the size of the cloud, however, then our line strengths should be increased by a factor of (see Section 2.2 for a discussion on how relates to the driving scale of turbulence).

Another factor of two uncertainty in the shock-integrated intensities comes from the assumption that the energy in magnetic field fluctuations in the cloud is negligible. While some MHD simulations have shown that the turbulent kinetic energy dominates over the energy in magnetic field fluctuations, particularly for smaller initial magnetic fields (Padoan & Nordlund 1999; Padoan et al. 2000; Heitsch et al. 2001), other simulations find that these magnetic field fluctuations have energy on the order of the kinetic energy of the turbulence (Stone et al. 1998; Ostriker et al. 2001). If the energy in magnetic field fluctuations is on a par with the kinetic energy of turbulence, the line-integrated intensities would have to be scaled up by a factor of two to account for the dissipation of this additional energy.

In scaling up the shock models to estimate the total integrated intensities from an entire cloud, the integrated intensity of every line in a particular model has been multiplied by the same factor and no optical depth effects have been taken into account. These optical depth effects, however, should not be significant for the higher rotational transitions of CO, where the CO transitions are effective at tracing shock emission, as the PDR models indicate that the CO lines are only optically thick up to, and including, the J = 5 4 transition. As for shock emission in the lower rotational transitions of CO, this emission is likely to be absorbed by the ambient gas and thereby will serve as a heating source for the ambient gas. Note, however, that the expected CO-integrated intensity at these low transitions is much less than the PDR-integrated intensity, which indicates that this extra shock heating will have only a very minor effect on the temperature, and thus the spectrum, of the ambient gas.

In deriving the total turbulent energy dissipation rate of a molecular cloud, it was assumed that turbulence decays in a crossing time. Recently, Basu et al. (2009) and Basu & Dapp (2010) discovered a long-lived magnetic-tension-driven mode in their thin disk simulations of flattened molecular clouds, arising from interactions between the disk and an external magnetic field, which was able to preserve a significant fraction of the turbulent energy of the cloud for much longer than the crossing time. The existence of such a long-lived MHD mode would significantly reduce the required energy dissipation rate of a molecular cloud and, therefore, our predicted line-integrated intensities as well. This long-lived mode, however, has not been noticed in any further simulations, including the three-dimensional simulations of collapsing cores done by Kudoh & Basu (2011), who did look for this particular MHD mode. Kudoh & Basu (2011) suggest that the absence of this mode in their simulations may be due to the very small density contrast between the disks and surrounding gas in their simulations.

The line ratios from these shock models are independent of the scaling of the lines and thus, are not affected by the above uncertainties. The line ratios of the shock models are also independent of the assumption of spherical geometry.

The turbulent energy of molecular clouds may also not be dissipated completely through shocks. Smith et al. (2000b) and Stone et al. (1998) find, in their simulations of turbulence, that only 50% of the turbulent energy is dissipated through artificial viscosity, due to the presence of shock fronts, while the other 50% is dissipated through numerical viscosity, representative of small-scale dissipation distributed relatively uniformly across the cloud. It is possible that some of this uniformly dissipated energy should have been dissipated in weak shocks or vortices that were not resolved in these simulations and thus, we consider that 50% is only a lower limit for the fraction of energy dissipated in intermittent structures (i.e., shocks).

An alternative model for turbulent dissipation, where energy is dissipated through magnetized vortices, has also been put forward (Godard et al. 2009). In the Godard et al. (2009) turbulent dissipation region (TDR) model, small vortices on the order of a few tens of AU heat gas to temperatures of nearly 1000 K via ion-neutral friction. The spectral signature of such a TDR model should be significantly different from our low-velocity shock models because the TDR model produces temperatures much higher than what can be obtained from slow shocks.

All of the shock models have been run with standard, roughly solar, metal abundances. In lower metallicity clouds, the abundance of CO will be reduced. The gas phase CO abundance will also be lower at higher densities, due to freeze out of CO onto dust grains (e.g., Goldsmith 2001), and in the outer layers of molecular clouds, where CO is readily photodissociated (e.g., Wolfire et al. 2010). The reduction of gas phase CO may lead to higher post shock temperatures, which would change the resulting CO spectrum and could affect which cooling mechanism is dominant. H line cooling and the deposition of energy into magnetic fields may also be more important for dissipating kinetic energy in CO poor gas. Further low-metallicity shock models, however, are needed to confirm this.

The CO spectrum predicted from the PDR model is dependent upon the chemistry put into the models. In particular, any chemical effects which would increase the temperature in the outer CO layers, such as an increase in polycyclic aromatic hydrocarbon heating, would shift the PDR spectrum toward higher rotational transitions. The PDR models used also have a relatively low ISRF of 3 G. In active, high-mass star-forming regions, the ambient ISRF is likely to be much higher, due to the presence of previously formed, massive, young stars and thus, the PDR emission from these regions is likely to be much stronger and peaked toward much higher J transitions. As such, our model comparisons should only be used for relatively quiescent, low-mass star-forming regions in which the ambient ISRF is relatively low, rather than in strong UV environments such as the Orion Bar.

4.7. Observational Potential

The total energy dissipated by shocks cannot be directly observed because some of the turbulent energy is not radiated away, but rather goes into increasing magnetic field strengths. Furthermore, many of the low lying CO lines are not readily observable, because these transitions are dominated by emission from unshocked gas. This ambient gas emission, however, does drop rapidly at higher J numbers and the CO J = 6 5 (691.47308 GHz) and 7 6 (806.651806 GHz) transitions are dominated by shock emission in most of our models. Thus, these higher rotational CO transitions act as shock tracers and by fitting shock models to the observed strengths of multiple high J CO transitions, the total shock luminosity of a cloud can be estimated.

The Herschel Space Observatory’s Heterodyne Instrument for the Far Infrared (HIFI) and Spectral and Photometric Imaging Receiver (SPIRE) both have the necessary wavelength coverage and sensitivity to be able to detect the CO J = 5 4, 6 5, and 7 6 lines, if they are as bright as we predict in our shock models. These two Herschel instruments also cover the J = 8 7 wavelength but our predicted line strengths for this transition are at the limits of what could be detected within a few hours of observing time. From the ground, the James Clerk Maxwell Telescope’s (JCMT) receiver W is capable of observing at the wavelength of the CO J = 6 5 transition while the Atacama Pathfinder Experiment (APEX) and the Caltech Submillimeter Observatory (CSO) have instruments that operate in the appropriate wavelength regimes to detect both the 6 5 and 7 6 lines. The CO J = 6 5 and 7 6 lines also lie within bands 9 and 10, respectively, of the Atacama Large Millimeter Array (ALMA). ALMA, with its superb resolution, may resolve individual shock fronts and thus, provide information regarding the properties of individual shocks.

Some care must be taken when choosing a location in a molecular cloud to observe, as there are other sources of high J CO line emission that have not included in the comparison PDR models. High-velocity protostellar outflows will generate large shocks that can easily produce temperatures of hundreds of Kelvin. Such strongly shocked gas will radiate significantly in high J lines (e.g., Kaufman & Neufeld 1996b, a). Embedded high-mass stars will also significantly heat nearby gas and lead to high J line emission. Thus, while turbulent dissipation should occur in active, high-mass star-forming regions, the spectral signatures of low-velocity, turbulence-induced shocks may only be readily detectable in quiescent regions of low-mass star-forming molecular complexes. Water emission can be used as a discriminant between low-velocity, turbulent shocks and the stronger shocks produced by outflows, as little water emission is predicted from our models while stronger shocks are expected to cool significantly via water lines (e.g., Kaufman & Neufeld 1996a).

The temperature sensitivity of the S(0) and S(1) lines of H, at 28.2 m and 17.0 m, respectively, makes these lines potential shock tracers. The atmospheric transmission at 28.2 m is, however, very poor, and thus, it is extremely difficult to observe the S(0) line with ground-based facilities. The expected weakness of these two lines makes observations of H emission even more problematic.

5. Global Heating

5.1. Cosmic Rays

In the cold, well-shielded, central regions of molecular clouds, cosmic-ray heating is believed to be the dominant heating term (e.g., Goldsmith 2001), but the general cosmic-ray ionization rate in the galaxy is not particularly well known. Estimates for the cosmic-ray ionization rate in dense gas vary from s (Caselli et al. 1998; Liszt 2003; Doty et al. 2004) to s (Caselli et al. 1998, 2002; Flower et al. 2007; Hezareh et al. 2008) but commonly lie between s and s (Williams et al. 1998; van der Tak & van Dishoeck 2000; Doty et al. 2002; Wakelam et al. 2005; Bergin et al. 2006; Maret & Bergin 2007; Goicoechea et al. 2009). While the spread in cosmic-ray ionization rates may simply be due to spatially varying rates (e.g., van der Tak et al. 2006), some of the scatter is likely due to uncertainties in the chemical models used to determine the ionization rates, as changes in known reaction rates (e.g., Dalgarno 2006) and the inclusion of non-equilibrium chemistry (Lintott & Rawlings 2006) have changed the results of ionization rate calculations.

For this paper, a cosmic-ray ionization rate range of s to s is adopted and, since each cosmic-ray ionization deposits approximately 20 eV into the gas (Goldsmith 2001), the volume heating rate by cosmic-rays is taken to be

(19)

This cosmic-ray heating rate is shown in Figure 5.

5.2. Turbulent Heating

If the turbulent energy of a cloud is not dissipated through localized shocks, but rather is dissipated relatively uniformly across the cloud, then this turbulent energy dissipation will act more as a general heating mechanism for the cloud and will yield the heating rate given by Equation (7). For the sole purpose of comparing the turbulent energy dissipation rate to the cosmic-ray heating rate, both the size-velocity and density-size relationships from Solomon et al. (1987) are used to rewrite the turbulent energy dissipation rate in terms of only and the gas density. Using these two scaling relations, the turbulent energy heating rate can be written as

(20)

where is the fraction of the dissipated turbulent energy that acts as a global heating mechanism. The simulations of Smith et al. (2000b) and Stone et al. (1998) both suggest that roughly half of the turbulent energy of a cloud may be dissipated uniformly across a cloud instead of in localized shocks and thus, suggest that (see Section 4.6 for further discussion on ). Furthermore, the magnetic energy injected by shocks may also decay and lead to a global heating of the cloud. Figure 5 shows the turbulent heating rates, as a function of density, for and . As before, a fixed value of 1 is used.

5.3. Ambipolar Diffusion

One potentially important energy dissipation mechanism that is usually not included in simulations (e.g., Stone et al. 1998; Mac Low 1999) is ambipolar diffusion. From their numerical simulations, Padoan et al. (2000) find that the ambipolar heating rate in well shielded molecular gas is given by

(21)

where is the volume-averaged magnetic field strength, is the Alfvénic Mach number of the turbulence, and is the volume-averaged number density. Applying Larson’s laws and a density-magnetic field strength scaling relation with k = 0.5 converts this equation to the form

(22)

Figure 5 shows this ambipolar diffusion heating rate, as a function of density, evaluated for the two b parameters previously used for our shock models, b = 0.1 and 0.3.

Figure 5.— Various heating rates for well-shielded molecular gas. The shaded (green) region shows the range of cosmic-ray heating rates, the dark (purple) solid line shows the total shock turbulent energy dissipation rate, the light (blue) solid line shows 50% of the shock turbulent dissipation rate, and the dark (red) and light (yellow) dotted lines show the ambipolar diffusion heating rates for b values of 0.3 and 0.1, respectively. Please see the online version for a color version of this figure.

5.4. Discussion

For gas with a density larger than 10 cm, the cosmic-ray heating rate is the dominant heating term. The turbulent energy dissipation rate is comparable to the cosmic-ray heating rate around a density of 10 cm and becomes larger than the cosmic-ray heating rate at lower densities.

Goldsmith (2001) examines the thermal balance of molecular clouds and finds that the ambient temperature of well-shielded gas, at a density of 10 cm, is 10 K. We took the Goldsmith (2001) prescriptions for cooling and heating and increased the cosmic-ray heating by a factor of two to reproduce the effect of having a turbulent energy dissipation heating rate equivalent to the cosmic-ray heating rate. With this increased heating rate, we find a slightly higher equilibrium temperature of 13 K. Pan & Padoan (2009) present a more detailed thermal balance model that also includes a turbulent heating term and, for a 1 pc sized cloud, find similar gas temperatures of 13-17 K, depending upon the cosmic-ray heating rate used.

This small change in temperature is well within the intrinsic scatter of observed gas temperatures in molecular clouds (e.g., Bergin & Tafalla 2007). Since the cosmic-ray heating rate is also uncertain by at least a factor of two, the finding of previous thermal balance studies that the ambient temperatures of molecular clouds are roughly consistent with heating by cosmic-ray ionization alone (e.g., Bergin et al. 2006) is not in conflict with the presence of heating from turbulent dissipation at the rate calculated above. The above agreement between observations and models does, however, constrain the turbulent heating rate to not be significantly greater than estimated above. This implies that cannot be much less than one, similar to the findings of Basu & Murali (2001). Padoan et al. (2000) also note that if turbulent heating is significant in a molecular cloud, then a positive temperature-velocity dispersion relation is expected, as observed by Jijina et al. (1999).

The significance of ambipolar diffusion to the thermal balance of a cloud is highly dependent upon the strength of the magnetic field. For a strong field, b = 0.3, ambipolar diffusion should be the dominant heating process in well-shielded gas with an average density of 100 cm, and should be comparable to other heating processes at a density of 10 cm. Gas at densities of 100 cm, however, may not necessarily be well shielded and may instead be dominated by the ISRF. If magnetic field strengths are slightly lower, corresponding to b = 0.1, the ambipolar diffusion rate is negligible for all densities greater than or equal to 100 cm. It should also be noted that the ambipolar diffusion heating rate becomes larger at lower densities, unlike the turbulent and cosmic-ray heating rates, because the typical collision speed between ions and neutrals increases with decreasing density.

Shocks will cause both the density and magnetic field strength to increase. The change in the density will be proportional to the change in the magnetic field strength, given the flux freezing approximation made in the shock models. The Alfvén speed in the gas will thus also increase by a factor proportional to the square root of the change in the magnetic field strength. If the magnetic field strength is locally increased by a factor of q by a shock and the velocity dispersion of a cloud remains relatively unchanged, then the ambipolar diffusion heating rate in the shocked gas will be

(23)

where is the ambipolar diffusion dissipation rate in the unshocked gas. Thus, the postshock gas will have its ambipolar diffusion heating rate increased.

Table 3 gives the factors by which the magnetic field strength and the ambipolar diffusion heating rate increase in each of the shock models. The magnetic field increases in strength by a factor of 4-22 and the ambipolar diffusion heating rate increases by one to two orders of magnitude. Even for the weakest magnetic field models, the ambipolar diffusion heating rate in the shocked gas is greater than the rate at which energy is injected into the magnetic field through the shock and thus, this enhanced ambipolar diffusion rate may provide an efficient mechanism for dissipating the injected magnetic energy. Furthermore, this enhanced ambipolar diffusion rate should add an extra heating source to the shocked gas, thereby increasing the CO flux that will emerge at higher J transitions.

The above discussion regarding heating rates is only relevant for the well-shielded centers of molecular clouds, as photoelectric heating due to the ISRF will dominate the heating in the outer PDR zones of molecular clouds (e.g., Kaufman et al. 1999). The heating rates shown in Figure 5 were also derived using Larson’s laws, which, as discussed in Section 2.1, have recently been called into question.

6. Conclusions

We have run models of MHD, C-type shocks, based on Kaufman & Neufeld (1996b), for shock velocities of 2 and 3 km s and initial densities between 10 and 10 cm. CO is found to be the dominant molecular coolant with 40%-80% of the shock energy being emitted in CO rotational lines. All other line cooling processes are negligible, except for H line cooling in the models with the very strongest shocks, in which H cooling accounts for 5%-20% of the total shock cooling. Between 20% and 60% of the shock energy also goes into compressing the magnetic field.

The expected CO spectrum from each of the shock models has been calculated and PDR models, based on Kaufman et al. (1999), were used to determine the expected contribution of CO emission from not only the cold, well-shielded interior of a molecular cloud, but also from the warm outer layers of the cloud. The PDR emission dominates for low J transitions of CO but the shock emission is larger at mid to high J transitions. In all models the shock emission is larger than the PDR emission in the J = transition and the shock emission can dominate as low as the J = 5 4 transition, depending upon the shock model. The J = 6 5 and 7 6 should serve as shock tracers and should be detectable with current observational facilities.

The turbulent energy dissipation rate is larger than the cosmic-ray heating rate for densities less than 10 cm. The presence of such an additional heating term is, however, still consistent with previous thermal balance studies, given the uncertainty in the cosmic-ray heating rate and the range of observed molecular cloud temperatures. The ambipolar diffusion heating rate is negligible at high densities and low-magnetic field strengths but can be dominant at lower densities and higher magnetic field strengths. Ambipolar diffusion is also enhanced in the shocked gas and may provide a mechanism for the dissipation of energy injected into the magnetic field by a shock.

We thank Shantanu Basu for helpful suggestions on the role of magnetic fields in molecular clouds as well as our anonymous referee for many useful changes to this paper. A.P. is partially supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) graduate scholarship program and D.J. is supported by a NSERC Discovery Grant. This research has made use of NASA’s Astrophysics Data System.

Footnotes

  1. affiliation: Department of Physics and Astronomy, University of Victoria, P.O. Box 3055, STN CSC, Victoria, BC V8W 3P6, Canada; arpon@uvic.ca
  2. affiliation: NRC-Herzberg Institute of Astrophysics, 5071 West Saanich Road, Victoria, BC V9E 2E7, Canada; Douglas.Johnstone@nrc-cnrc.gc.ca
  3. affiliation: NRC-Herzberg Institute of Astrophysics, 5071 West Saanich Road, Victoria, BC V9E 2E7, Canada; Douglas.Johnstone@nrc-cnrc.gc.ca
  4. affiliation: Department of Physics and Astronomy, University of Victoria, P.O. Box 3055, STN CSC, Victoria, BC V8W 3P6, Canada; arpon@uvic.ca
  5. affiliation: Department of Physics, San Jose State University, One Washington Square, San Jose, CA 95192-0106, USA; mkaufman@email.sjsu.edu
  6. affiliation: Space Science and Astrobiology Division, MS 245-3, NASA Ames Research Center, Moffett Field, CA 94035, USA

References

  1. André, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102
  2. Arce, H. G., Borkin, M. A., Goodman, A. A., Pineda, J. E., & Halle, M. W. 2010, ApJ, 715, 1170
  3. Arons, J., & Max, C. E. 1975, ApJ, 196, L77
  4. Ballesteros-Paredes, J., & Mac Low, M.-M. 2002, ApJ, 570, 734
  5. Banerjee, R., Klessen, R. S., & Fendt, C. 2007, ApJ, 668, 1028
  6. Basu, S., Ciolek, G. E., Dapp, W. B., & Wurster, J. 2009, New A, 14, 483
  7. Basu, S., & Dapp, W. B. 2010, ApJ, 716, 427
  8. Basu, S., & Murali, C. 2001, ApJ, 551, 743
  9. Bergin, E. A., Maret, S., van der Tak, F. F. S., et al. 2006, ApJ, 645, 369
  10. Bergin, E. A., & Tafalla, M. 2007, ARA&A, 45, 339
  11. Blitz, L., & Shu, F. H. 1980, ApJ, 238, 148
  12. Brunt, C. M. 2003, ApJ, 583, 280
  13. Brunt, C. M., Heyer, M. H., & Mac Low, M. 2009, A&A, 504, 883
  14. Caselli, P., Walmsley, C. M., Terzieva, R., & Herbst, E. 1998, ApJ, 499, 234
  15. Caselli, P., Walmsley, C. M., Zucconi, A., et al. 2002, ApJ, 565, 344
  16. Chernoff, D. F., McKee, C. F., & Hollenbach, D. J. 1982, ApJ, 259, L97
  17. Crutcher, R. M. 1999, ApJ, 520, 706
  18. Crutcher, R. M., Wandelt, B., Heiles, C., Falgarone, E., & Troland, T. H. 2010, ApJ, 725, 466
  19. Curtis, E. I., Richer, J. S., Swift, J. J., & Williams, J. P. 2010, MNRAS, 408, 1516
  20. Dalgarno, A. 2006, Proceedings of the National Academy of Science, 103, 12269
  21. Di Francesco, J., Evans, II, N. J., Caselli, P., et al. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil (Tucson, AZ: Univ. Arizona Press), 17
  22. Doty, S. D., Schöier, F. L., & van Dishoeck, E. F. 2004, A&A, 418, 1021
  23. Doty, S. D., van Dishoeck, E. F., van der Tak, F. F. S., & Boonman, A. M. S. 2002, A&A, 389, 446
  24. Draine, B. T. 1980, ApJ, 241, 1021
  25. Draine, B. T., & McKee, C. F. 1993, ARA&A, 31, 373
  26. Draine, B. T., Roberge, W. G., & Dalgarno, A. 1983, ApJ, 264, 485
  27. Elitzur, M., & de Jong, T. 1978, A&A, 67, 323
  28. Elitzur, M., & Watson, W. D. 1978, A&A, 70, 443
  29. Elmegreen, B. G. 1985, ApJ, 299, 196
  30. —. 2000, ApJ, 530, 277
  31. Eng, C. 2002, PhD thesis, University of Illinois at Urbana-Champaign
  32. Fiedler, R. A., & Mouschovias, T. C. 1993, ApJ, 415, 680
  33. Field, G. B. 1978, in IAU Colloq. 52: Protostars and Planets, ed. T. Gehrels, 243
  34. Flower, D. R., Pineau Des Forêts, G., & Walmsley, C. M. 2007, A&A, 474, 923
  35. Gammie, C. F., & Ostriker, E. C. 1996, ApJ, 466, 814
  36. Glover, S. C. O., Federrath, C., Mac Low, M., & Klessen, R. S. 2010, MNRAS, 404, 2
  37. Godard, B., Falgarone, E., & Pineau Des Forêts, G. 2009, A&A, 495, 847
  38. Goicoechea, J. R., Pety, J., Gerin, M., Hily-Blant, P., & Le Bourlot, J. 2009, A&A, 498, 771
  39. Goldreich, P., & Kwan, J. 1974, ApJ, 189, 441
  40. Goldsmith, P. F. 2001, ApJ, 557, 736
  41. Goldsmith, P. F., Velusamy, T., Li, D., & Langer, W. D. 2010, ApJ, 715, 1370
  42. Hartmann, L., Ballesteros-Paredes, J., & Bergin, E. A. 2001, ApJ, 562, 852
  43. Heitsch, F., Mac Low, M., & Klessen, R. S. 2001, ApJ, 547, 280
  44. Heyer, M., Krawczyk, C., Duval, J., & Jackson, J. M. 2009, ApJ, 699, 1092
  45. Heyer, M. H., & Brunt, C. M. 2004, ApJ, 615, L45
  46. Hezareh, T., Houde, M., McCoey, C., Vastel, C., & Peng, R. 2008, ApJ, 684, 1221
  47. Hollenbach, D., & McKee, C. F. 1989, ApJ, 342, 306
  48. Jijina, J., Myers, P. C., & Adams, F. C. 1999, ApJS, 125, 161
  49. Kaufman, M. J., & Neufeld, D. A. 1996a, ApJ, 456, 611
  50. —. 1996b, ApJ, 456, 250
  51. Kaufman, M. J., Wolfire, M. G., & Hollenbach, D. J. 2006, ApJ, 644, 283
  52. Kaufman, M. J., Wolfire, M. G., Hollenbach, D. J., & Luhman, M. L. 1999, ApJ, 527, 795
  53. Kudoh, T., & Basu, S. 2011, ApJ, 728, 123
  54. Larson, R. B. 1981, MNRAS, 194, 809
  55. Lepp, S., & Shull, J. M. 1983, ApJ, 270, 578
  56. Lintott, C. J., & Rawlings, J. M. C. 2006, A&A, 448, 425
  57. Liszt, H. 2003, A&A, 398, 621
  58. Mac Low, M. 1999, ApJ, 524, 169
  59. Mac Low, M., & Klessen, R. S. 2004, Reviews of Modern Physics, 76, 125
  60. Mac Low, M., Klessen, R. S., Burkert, A., & Smith, M. D. 1998, Physical Review Letters, 80, 2754
  61. Maret, S., & Bergin, E. A. 2007, ApJ, 664, 956
  62. McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565
  63. Mestel, L. 1966, MNRAS, 133, 265
  64. Mouschovias, T. C., & Spitzer, Jr., L. 1976, ApJ, 210, 326
  65. Mullan, D. J. 1971, MNRAS, 153, 145
  66. Neufeld, D. A., & Kaufman, M. J. 1993, ApJ, 418, 263
  67. Ossenkopf, V., & Mac Low, M. 2002, A&A, 390, 307
  68. Ostriker, E. C., Stone, J. M., & Gammie, C. F. 2001, ApJ, 546, 980
  69. Padoan, P., Juvela, M., Kritsuk, A., & Norman, M. L. 2009, ApJ, 707, L153
  70. Padoan, P., & Nordlund, Å. 1999, ApJ, 526, 279
  71. Padoan, P., Zweibel, E., & Nordlund, Å. 2000, ApJ, 540, 332
  72. Pan, L., & Padoan, P. 2009, ApJ, 692, 594
  73. Quillen, A. C., Thorndike, S. L., Cunningham, A., et al. 2005, ApJ, 632, 941
  74. Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369
  75. Shu, F. H. 1977, ApJ, 214, 488
  76. Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23
  77. Smith, M. D., Mac Low, M., & Heitsch, F. 2000a, A&A, 362, 333
  78. Smith, M. D., Mac Low, M., & Zuev, J. M. 2000b, A&A, 356, 287
  79. Solomon, P. M., Rivolo, A. R., Barrett, J., & Yahil, A. 1987, ApJ, 319, 730
  80. Stone, J. M., Ostriker, E. C., & Gammie, C. F. 1998, ApJ, 508, L99
  81. Tielens, A. G. G. M. 2005, The Physics and Chemistry of the Interstellar Medium (Cambridge, UK: Cambridge University Press)
  82. Timmermann, R. 1996, ApJ, 456, 631
  83. van der Tak, F. F. S., Belloche, A., Schilke, P., et al. 2006, A&A, 454, L99
  84. van der Tak, F. F. S., & van Dishoeck, E. F. 2000, A&A, 358, L79
  85. Wakelam, V., Selsis, F., Herbst, E., & Caselli, P. 2005, A&A, 444, 883
  86. Williams, J. P., Bergin, E. A., Caselli, P., Myers, P. C., & Plume, R. 1998, ApJ, 503, 689
  87. Williams, J. P., & McKee, C. F. 1997, ApJ, 476, 166
  88. Wolfire, M. G., Hollenbach, D., & McKee, C. F. 2010, ApJ, 716, 1191
  89. Zuckerman, B., & Evans, II, N. J. 1974, ApJ, 192, L149
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
174733
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description