© ESO 2018

1 Introduction

Carbon monoxide, CO, is a critical molecule in the physics of molecular clouds and star formation. First, it is believed to be the next most abundant molecule in the Universe after molecular hydrogen, H₂. Second, as compared with H₂, it has a relatively high moment of inertia, and a permanent dipole moment, so it emits readily at the low temperatures found in molecular clouds, whereas H₂ does not; therefore CO and its isotopologues are often used to trace the structure and dynamics of molecular clouds (e.g. Roman-Duval et al. 2016). Third, because it emits readily, it plays an important role in the thermal balance of molecular clouds, helping to maintain their pervasive low temperatures. However, the physics underlying the formation/destruction of CO, and the net line emission from CO, is complicated.

The pathways to CO formation are very diverse and uncertain; they depend on the formation of several other species, on reaction rates that are only known approximately, and on the agency of cosmic rays to produce molecular ions. Likewise, the destruction of CO is difficult to model, largely because it involves line radiation; as a consequence, evaluating self-shielding is difficult, and is made more difficult because self-shielding has to compete with some lines also being blocked by H₂, with dust attenuation, and, at high densities, with CO freeze-out. A variety of schemes has been proposed to estimate the abundance of CO without considering all these details (e.g. Nelson & Langer 1997, 1999), and indeed these approximate schemes are the ones used in many interstellar chemistry codes (e.g. Glover & Clark 2012a). In the Solar vicinity, it seems that gas-phase CO is usually the dominant form of carbon at volume densities in the range $X\text{\_{\_CO}}$ and column-densities $X\text{\_{\_CO}}$. At lower volume- and column-densities, CO gives way to atomic and ionic carbon (C^o and C⁺), due to the slow rate of the two-body reactions leading to CO formation, and the strong ultraviolet irradiation which rapidly destroys CO. At higher volume-densities, CO appears to freeze out onto dust; parenthetically, CO’s role as a coolant also becomes less important at these higher densities, because the gas starts to couple thermally to the dust, and hence to cool by continuum emission from dust.

Line emission from CO is complicated by the fact that, in the range of volume- and column-density where CO tends to be abundant, transfer of population between the different levels on the rotational ladder of CO involves a balance between collisional and radiative excitation and de-excitation, and the details of this balance shift with changing volume-density, column-density, temperature and velocity dispersion. CO line emission is seldom in either of the asymptotic limits of very low volume- and column-density (hence low line-centre optical depth), or very high volume- and column-density (hence high line-centre optical depth), where its collective line emission can be described by simple algebraic equations (e.g. Goldreich & Kwan 1974).

Here we derive an approximate analytic formulation for the net cooling rate from CO and its isotopologues. The formulation is simple in the sense that it is a purely local function of state, depending only on four parameters: (i) the mass-density, ρ, or equivalently the volume-density of molecular hydrogen, $n{\text{\_{\_H}}_2\text{}}$; (ii) the isothermal sound speed, a, or equivalently the gas-kinetic temperature, T; (iii) the relative abundance of CO, $X\text{\_{\_CO}}\equiv$ $n\text{\_{\_CO}}/n{\text{\_{\_H}}_2\text{}}$; and (iv) the velocity divergence, ∇ ⋅υ. The formulation is obtained by first deriving the simple algebraic equations that describe the dependence of the CO cooling rate on density, temperature, CO abundance and velocity divergence, in the asymptotic limit of low volume-density and optical depth, and in the asymptotic limit of high volume-density and optical depth; hereafter we refer to these equations as the asymptotic forms. Their derivation is based on fundamental physical arguments, and forms the subject of Sect. 2. Then, in Sect. 3 the coefficients in front of the asymptotic forms, and the variation between them, are fit by reproducing the detailed results of GL78, using an ad hoc mathematical function. In Sect. 4 we illustrate the application of the approximate analytic formulation by using it to evaluate the thickness of the post-shock cooling layer behind a low-velocity steady J-shock in a non-magnetic molecular cloud. In Sect. 5 we summarise our conclusions.

For mathematical convenience, we use the mass-density, ρ and the isothermal sound speed, $a={(k\text{\_{\_B}}T/\overline{m})}^{1/2}$ (in place of the number-density of molecular hydrogen, $n{\text{\_{\_H}}_2\text{}}$, and the gas-kinetic temperature T) in much of the analysis. Here $k\text{\_{\_B}}$ is Boltzmann’s constant and $\overline{m}$ is the mean gas-particle mass. For the purpose of illustration, we use a reference temperature To = 10 K and a mean gas-particle mass $X\text{\_{\_CO}}$ (appropriate for molecular gas with elemental composition X = 0.70, Y = 0.28, Z = 0.02);¹ hence the reference isothermal sound speed is ao = 0.187 km s⁻¹. On the assumption that virtually all the hydrogen is in the form of H₂ , we also define the mean mass per hydrogen molecule, $X\text{\_{\_CO}}$, and hence $\rho =n{\text{\_{\_H}}_2\text{}}\overline{m}{\text{\_{\_H}}_2\text{}}$.

2 The CO cooling rate: asymptotic theory

In this section we derive, using basic physical arguments, the dependence of the CO line cooling rate on the density, on the isothermal sound speed (or equivalently the temperature), on the gas-phase CO abundance, and on the local velocity divergence. We treat the two asymptotic limits of (i) very low volume-density and optical depth (hereafter the “lo” limit), and (ii) very high volume-density and optical depth (hereafter the “hi” limit). We do not derive the absolute physical value of the cooling rate in these limits, but only the dependence on density, isothermal sound speed, abundance and velocity divergence. The coefficients converting these dependences into actual physical values are obtained in the following section (Sect. 3) by comparing the theoretical predictions with the detailed results of GL78.

2.1 CO cooling in the lo limit

For a linear molecule with moment of inertia $I\text{\_{\_MOL}}$, the rotational levels are characterised by quantum number J, and have energy $X\text{\_{\_CO}}$ therefore, if the gas-kinetic temperature is T, they are significantly excited up to $$X\text{\_{\_CO}}$$(1)

where $f\text{\_{\_EX}}$ is a factor greater than, but of order, unity. Radiative de-excitations from level J to level J − 1 release photons with energy $\text{UpDelta}E\text{\_{\_}J\text{}}=J{\hslash }^2/I\text{\_{\_MOL}}$. Therefore the mean energy of emitted photons is $$X\text{\_{\_CO}}$$(2)

Here, we have replaced a sum over levels with an integral. This is reasonable for higher temperatures, T≳20 K, but may introduce a small systematic error at lower temperatures where only a few levels are involved.

At sufficiently low volume-density, i.e. well below the critical density, the level populations are not thermalised. Most molecules sit in their ground state, and most collisional excitations are followed by radiative de-excitations. For example, at T ~ 10 K, the first rotationally excited level of CO, J = 1, has a critical density of order 10³ H_2 cm⁻³; for higher-J levels, higher temperatures and densities are required before they are thermalised. At sufficiently low optical depth (i.e. sufficiently low column-density and/or high velocity divergence), the lines are optically thin, so most of the emitted photons escape directly; consequently the thermal kinetic energy that caused the initial excitation is lost, and the gas is cooled. For example, at T ~ 10 K and with purely thermal velocity dispersion, i.e. no turbulence and no velocity gradient, the first rotational transition of CO (J = 1 → 0) becomes optically thick for column-densities $X\text{\_{\_CO}}$; this limit is increased if there is turbulence and/or a velocity gradient delivering extra velocity dispersion. The rate of collisional excitation per unit volume depends on the product of the number-densities of the molecule and the exciting collider (here presumed to be H₂ , see below) and their relative speed, so it is approximately proportional to $n\text{\_{\_MOL}}n{\text{\_{\_H}}_2\text{}}{T}^{1/2}$. The number of lines that are excited is approximately proportional to T^1∕2 (see Eq. (1)). And the mean energy of the emitted photons is approximately proportional to T^1∕2 (see Eq. (2)). Therefore, in the limit of very low volume-density and optical depth, the cooling rate per unit volume can be approximated by $$X\text{\_{\_CO}}$$(2)

In treating cooling by CO, we presume that H₂ is the dominant agent of collisional excitation, since, in regions where there is CO, H₂ is likely to be by far the most abundant species; GL78 make the same assumption.

On the basis of Eq. (3), we posit that, in the lo limit, the cooling rate due to CO is given by $$X\text{\_{\_CO}}$$(4)

The coefficient $\lambda \text{\_{\_CO}.\text{LO}}$ will be estimated in Sect. 3, but the above dependence of $\text{UpLambda\_{\_CO}.\text{LO}}$ on $X\text{\_{\_CO}}$, ρ and a is now fixed.

2.2 CO cooling in the hi limit

At sufficiently high volume-density, i.e. well above the critical density, the level populations are approximately thermalised, and radiative de-excitations account for only a small fraction of de-excitations. At sufficiently high optical depth (i.e. sufficiently high column-density and/or low velocity divergence), the intensity in the emission line cores is approximately equal to the Planck function, and the net frequency width of all the excited lines is $$X\text{\_{\_CO}}$$(5)

Here Δυ is the velocity spread along the line of sight and $f\text{\_{\_LW}}$ is the factor by which the effective width of an optically thick line with central frequency νo differs from² νo Δυ∕c. The effective width of the Planck spectrum is $\text{UpDelta}\nu \text{\_{\_BB}}=f\text{\_{\_BB}}k\text{\_{\_B}}T/h$ where $f\text{\_{\_BB}}$ = 14.4, and so the total cooling rate per unit volume in this limit can be approximated by $$X\text{\_{\_CO}}$$

In deriving Eq. (6), we have adopted the most generic geometry we could imagine, viz. a spherical cloud of diameter D expanding or collapsing homologously, so that Δυ∕D ≡|∇⋅υ|∕3; GL78 adopt the same configuration. However, in obtaining Eq. (7), this choice becomes immaterial, since we only retain the dependence on a and ∇ ⋅ υ. The effect of geometry is encapsulated entirely in the ∇⋅υ term, which measures how far, in different directions on the sky, line-photons have to travel before the Doppler-shift has separated their frequency from the frequencies that can be absorbed by the molecules they are passing.

In effect we are adopting the large velocity gradient (LVG) approximation to treat optically thick cooling radiation. The LVG approximation was introduced by Sobolev (1960), and developed further by Castor (1970) and Lucy (1971), in the context of stellar outflows, but it can also be applied to the cooling of molecular clouds (e.g. Goldreich & Kwan 1974). Basically it assumes that line radiation only experiences self-absorption in a local region whose linear size is of order σ |∇ ⋅υ|⁻¹, where σ is the local velocity dispersion (thermal plus turbulent) and υ is the local bulk velocity. Any material outside this region is moving at a sufficiently different velocity (bulk plus or minus the dispersion) from the emitting region that it cannot absorb its line radiation. Strictly, the LVG approximation requires that the large-scale spatial variation of the bulk velocity is monotonic along any line of sight.

On the basis of Eq. (7), we posit that, in the hi limit, the CO cooling rate is given by $$X\text{\_{\_CO}}$$(8)

the coefficient $\lambda \text{\_{\_CO}.\text{HI}}$ will be estimated in Sect. 3, but the above dependence of $\text{UpLambda\_{\_CO}.\text{HI}}$ on a and |∇ ⋅ υ| is now fixed — and $\text{UpLambda\_{\_CO}.\text{HI}}$ is independentof ρ and $X\text{\_{\_CO}}$.

We note parenthetically that the only parameters in Eq. (6) that depend on the specific molecule are $f\text{\_{\_EX}}$ and $f\text{\_{\_LW}}$, and that this dependence is in general rather weak, so Eq. (6) is approximately valid for many other linear molecules in this high volume-density, high optical depth limit. Thus, for example, the results presented by GL78 for cooling by O₂ are very similar to those for CO (at the same density and temperature).

2.3 In between the lo and hi limits

In the next section (Sect. 3) we show how Eqs. (4) and (8) can be combined to obtain an approximate analytic formulation for the total CO cooling rate in the intermediate regime (between the lo and hi limits) using $$X\text{\_{\_CO}}$$(9)

with β = β(ρ, a).

3 The CO cooling rate: calibration

In the preceding section (Sect. 2) we have obtained expressions for the CO cooling rate (in the lo and hi limits) in terms of the mass-density, ρ, the isothermal sound speed, a, the abundance of CO, $X\text{\_{\_CO}}$, and the velocity divergence, ∇⋅υ. This is because these are the most convenient variables to use when treating the equations of hydrodynamics. Here, we replace ρ with the density of molecular hydrogen $n{\text{\_{\_H}}_2\text{}}$, and a with the temperature, T, since these are the variables used by GL78. It is straightforward to switch between the two, using $n{\text{\_{\_H}}_2\text{}}=\rho /\overline{m}{\text{\_{\_H}}_2\text{}}$ and $T=\overline{m}{a}^2/k\text{\_{\_B}}$.

3.1 The detailed results of Goldsmith & Langer (1978)

GL78 consider a uniform-density spherical cloud with a linear radial velocity field, and determine the net cooling rate, $\text{UpLambda\_{\_CO}.\text{TOT}}$, from CO and its isotopologues. To do this, they solve in detail (i) the equations of statistical equilibrium, to determine CO level populations; then (ii) the equations of radiative emission (spontaneous and stimulated), to determine the net CO emissivity; finally (iii) they use an escape probability in lieu of solving the equation of radiative transfer. In their Fig. 2, which is the input we use to calibrate our approximate analytic formulation, they plot the quantity $\log \text{\_{\_}10\text{}}\{\text{UpLambda\_{\_CO}.\text{TOT}}/n\text{\_{\_CO}}\}$ against $\log \text{\_{\_}10\text{}}\{n{\text{\_{\_H}}_2\text{}}/{\text{cm}}^{-3}\}$, for different values of T and different values of³ $X\text{\_{\_CO}}$, specifically T = 10, 20, 40, and 60K and $X\text{\_{\_CO}}$. We note that $\text{UpLambda\_{\_CO}.\text{TOT}}/n\text{\_{\_CO}}$ is the cooling rate per CO molecule. It depends on $X\text{\_{\_CO}}/|\nabla \cdot \upsilon |$, because – all other things being equal – the optical depth increases with increasing $X\text{\_{\_CO}}$ and decreases with increasing |∇ ⋅ υ|. In the lo limit, where the emission is optically thin, $\text{UpLambda\_{\_CO}.\text{TOT}}/n\text{\_{\_CO}}$ is independent of $X\text{\_{\_CO}}/|\nabla \cdot \upsilon |$.

Given the complexity of the physics underlying the net CO cooling rate, and the consequent computational cost of calculating it properly, it would be useful to have an approximate analytic formulation that captured the dependence on local variables but avoided the associated computational cost. This formulation could then be used in situations where CO cooling played an important role, but was not the principal interest, as for example, contracting pre-stellar cores, or the accretion shock at the boundary of an assembling filament.

3.2 The fitting procedure

To fit the GL78 results analytically, we have read from their plots the values of $\log \text{\_{\_}10\text{}}\{\text{UpLambda\_{\_CO}.\text{TOT}}/n\text{\_{\_CO}}\}$ for each treated combination of T, i.e. T = 10, 20, 40, and 60K, and $X\text{\_{\_CO}}$, i.e. $X\text{\_{\_CO}}$ and − 4, at the densities $X\text{\_{\_CO}}$, thus a total of 140 discrete values. We estimate that the uncertainty which derives from our reading these values of GL78’s Fig. 2 by eye is <0.03 in $\log \text{\_{\_}10\text{}}\{\text{UpLambda\_{\_CO}.\text{LO}}/n\text{\_{\_CO}}\}$ i.e. at worst ± 7%, and less than the thickness of the lines in Fig. 1.

We then use a Metropolis-Hastings Markov Chain Monte Carlo algorithm to find the five fitting parameters $({\text{UpLambda}}^{\prime }\text{\_{\_CO}.\text{LO}},{\text{UpLambda}}^{\prime }\text{\_{\_CO}.\text{HI}},{\beta }^{\prime }\text{\_{\_O}},{{\beta }^{\prime }}_{n_{{\text{H}}_2}},{{\beta }^{\prime }}_T)$ in the formulation $$X\text{\_{\_CO}}$$

that give the best fit to these discrete values.

We note that the exponents from Eqs. (4) and (8) are not being allowed to vary, they retain the values derived in Sect. 2. Because we are here fitting $\text{UpLambda\_{\_CO}.\text{TOT}}/n\text{\_{\_CO}}$ (rather than $\text{UpLambda\_{\_CO}.\text{TOT}}$), the ρ² term in Eq. (4) has become n_H2¹, and the ρ⁰ term implicit in Eq. (8) has become $n{\text{\_{\_H}}_2\text{}^}-1$. The a³ term in Eq. (4) has become T^3∕2 and the a⁸ term in Eq. (4) has become T⁴.

Fig. 1 Continuous curves give the total CO cooling rate, as predicted by the approximate analytic formulation derived here (Eqs. (15) through (18)). The different panels correspond to (a) T = 10K ; (b) T = 20K; (c) T = 40K; and (d) T = 60K (these temperatures are given in the top lefthand corner of each panel). The different curves represent different values of $X\text{\_{\_CO}}$ (these values are given down the righthand margin of each panel). The rates obtained by the detailed computations of GL78 are given at discrete values of the density, $X\text{\_{\_CO}}$; plus signs, stars, open circles, crosses and open squares correspond respectively to $X\text{\_{\_CO}}$.

3.3 The best fit and its accuracy

The best fit to the results from GL78 is obtained with the following values for the fitting parameters: $$X\text{\_{\_CO}}$$(14)

In otherwords, the best fit is $$X\text{\_{\_CO}}$$

The continuous lines in Fig. 1 show the predictions of Eqs. (15) through (18) for T = 10, 20, 40, and 60K (respectively panels a, b, c, d), and $X\text{\_{\_CO}}$ (separate lines labelled on the righthand margin of each panel). Like GL78 we limit the plots to densities in the range $1\le \log \text{\_{\_}10\text{}}\{n{\text{\_{\_H}}_2\text{}}/{\text{cm}}^{-3}\}\le 7$. The continuous curves should be compared with the discrete symbols, which give the values that we have read from Fig. 2 of GL78 and that we used in the Metropolis-Hastings fitting procedure.

In general, the agreement is good, with the magnitude of the fractional offset being everywhere <0.18 in $\log \text{\_{\_}10\text{}}\{\text{UpLambda\_{\_CO}.\text{LO}}/n\text{\_{\_CO}}\}\left(\text{i}\text{.e}{\text{., at worst, }}_{-34\%}^{+52\%}\right)$, and on average ≲0.07 in $\log \text{\_{\_}10\text{}}\{\text{UpLambda\_{\_CO}.\text{LO}}/n\text{\_{\_CO}}\}\left(\text{mbox}i.e.,typically{,}_{-15\%}^{+18\%}\right)$. The fractional offset is worst in the intermediate regime. This is as expected, since the analytic forms that obtain in the asymptotic limits (lo and hi) are physically motivated, and therefore in some sense absolute. In contrast, the intermediate regime between these limits is being fit with an algebraic expression that has no physical motivation beyond the fact that it relaxes to the asymptotic forms in the corresponding limits; it is therefore unable to produce an exact fit. The error in reading the discrete values of GL78’s Fig. 1 by eye ( ± 3%) is much smaller than the fractional offsets cited above, and can therefore be ignored.

There are some additional uncertainties that should be noted. First, GL78 used the coefficients for collisional de-excitation of CO by H₂ computed by Green & Thaddeus (1976). More recent computations by Yang et al. (2010) have increased these coefficients somewhat, at the lowest temperatures, T≲20 K. Therefore the GL78 cooling rates, and our approximate analytic formulation based on them, may be a little low at these low temperatures. Second, the device of replacing with integrals, summations over the discrete contributions from individual levels, will be most inaccurate at low temperatures where only a few levels are involved.

We note that the rarer isotopologues of CO make little contribution to the cooling in the lo limit, but do contribute once the ¹²C¹⁶O lines start to become optically thick, because the rarer isotopologues remain optically thin to much higher column-densities. They also make an important contribution in the hi limit, by providing extra bandwidth at high column-densities; in fact, in the hi limit, each isotopologue makes essentially the same contribution to the net cooling rate.

3.4 Use and range of applicability of the approximate analytic formulation

Given this level of accuracy, the approximate analytic formulation (i.e. Eqs. (15) through (18)) provides a convenient and computationally very inexpensive way to estimate the net CO cooling rate: convenient because it is a local function of state, and computationally inexpensive because it entails only a small number of arithmetic operations. It can therefore be used to greatly speed up numerical and/or analytic integrations of the energy equation in hydrodynamic or hydrostatic simulations, provided only that the conditions are within the range explored by GL78.

The numerical results of GL78 cover the density range $X\text{\_{\_CO}}$, the temperature range 10 K ≤ T ≤ 60 K, and the abundance/velocity-divergence range $X\text{\_{\_CO}}$, and therefore this is the range over which our approximate analytic formulation is most secure. Given the regularity of the curves in GL78, and their apparent self-similarity, it may be safe to extrapolate further. It is unlikely that gas-phase CO is the dominant form of CO outside the density range defined above, so we only consider extrapolation of the temperature and abundance/velocity-divergence ranges. Specifically, we presume that the approximate analytic formulation can be extended to somewhat higher temperatures, T ~100 K, and somewhat lower abundances or higher velocity divergences, $X\text{\_{\_CO}}$.

4 Slow steady planar J-shocks in non-magnetic molecular clouds

As a simple example of the application of our approximate analytic formulation, we consider the role of CO in post-shock cooling layers in molecular clouds. There are many other possible applications, for example the role of CO cooling in collapsing and fragmenting pre-stellar cores, and the thermodynamics of the accretion shock bounding a growing filament. We plan to investigate these in future papers.

Molecular clouds are observed to be highly turbulent, and supersonic velocity dispersions are the norm on scales ≳0.1 pc (Larson 1981; Solomon et al. 1987; Goodman et al. 1998; Caselli et al. 2002; Heyer & Brunt 2004). Shocks are therefore endemic in molecular clouds, and a critical issue is then how quickly the post-shock gas cools back to something like its pre-shock value. In particular, if post-shock cooling is very quick, shocks can for many purposes be treated as isothermal discontinuities, and this may greatly simplify analysis of the dynamics on larger scales.

In the interests of simplicity, we focus on slow, steady, planar J-shocks (e.g. Brand 1989; Flower et al. 2003). We note (i) that, if the shock velocity is high, υo≳15 km s^-1, the associated kinetic energy, and resulting post-shock thermal energy, may be sufficient to drive significant chemical changes; (ii) that these chemical changes may influence the dynamical coupling between the gas, the dust and the magnetic field; and (iii) that the cooling radiation from the hot gas immediately behind the shock may have important effects on the chemistry and thermal state of the gas flowing into the shock, as will the ambient radiation field. However, these three considerations are probably not important for the low-velocity shocks, υo ≲1.5 km s^-1, considered here.

4.1 The hydrodynamics of a J-shock plus cooling layer

We consider a slow steady planar J-shock at x = 0, with gas flowing in from x < 0, having density ρ(x < 0) = ρo, velocity parallel to the x axis υ(x < 0) = υo, and isothermal sound speed a(x < 0) = ao. There is no magnetic field, so we are ignoring the possibility that the shock is cushioned by a lateral magnetic field; in that case, the shock would be a C-shock, and the immediate post-shock temperature would be somewhat lower, but, as shown by Pon et al. (2012), the cooling would still very likely be dominated by CO. We are also ignoring the possibility that the pre-shock gas is significantly heated by the radiation from the post-shock cooling layer; this is a reasonable assumption, since gas in the post-shock cooling layer is necessarily at a very different velocity to the gas flowing into the shock, and therefore the CO cooling radiation cannot be absorbed by the inflowing gas.

We assume that the gas flowing through the shock is and remains molecular, with constant adiabatic exponent γ = 5∕3. This is strictly only true in low-temperature regions, T≲100 K, where the rotational degrees of freedom of molecular hydrogen are not significantly excited (Whitworth & Clarke 1997; Boley et al. 2007). This requirement is naturally met, since the GL78 results (and hence our approximate analytic formulation) are limited to temperatures T≲60 K. If the pre-shock temperature is ~10 K, this means we must restrict our model to shocks with υo≲1.0 km s^-1. If we extrapolate the approximate analytic formulation to T≲100 K, we can treat shocks with υo≲1.5 km s^-1.

The gas emerging from the shock has density, ρs, velocity, υs, and isothermal sound speed, as, given by the Rankine–Hugoniot conditions: $$X\text{\_{\_CO}}$$

where the second expressions apply provided υo² ≫ ao², and we recall that here ao is the isothermal sound speed.

In the post-shock region (x > 0), conservation of mass requires constant ρ(x)υ(x), hence $$X\text{\_{\_CO}}$$(22)

conservation of momentum requires constant $\rho (x)\left({\upsilon }^2(x)+a^2(x)\right)$, hence $$X\text{\_{\_CO}}$$(23)

and conservation of energy requires $$\begin{array}{lll}\frac{\text{d}}{\text{d}x}\left\{\rho (x)\upsilon (x)\left(\frac{{\upsilon }^2(x)}{2}+\frac{\gamma a^2(x)}{(\gamma -1)}\right)\right\}=\hfill & \hfill & \hfill \\ \rho \text{subO}\upsilon \text{subO}\left\{\frac{5}{2}\left(\upsilon \text{subO}+\frac{a{\text{subO}}^2}{\upsilon \text{subO}}\right)-4\upsilon (x)\right\}\frac{\text{d}\upsilon }{\text{d}x}\hfill & =\hfill & -\text{UpLambda}\left(\rho ,a,\left|\frac{\text{d}\upsilon }{\text{d}x}\right|\right),\hfill \end{array}$$(24)

where Λ(ρ, a, |dυ∕dx|) is the cooling rate per unit volume, and in planar geometry ∇⋅υ →dυ∕dx.

Since, in the post-shock cooling layer, υ(x)≲υo∕4, Eqs. (23) and (24) approximate to $$X\text{\_{\_CO}}$$

where we have substituted |dυ∕dx| = −dυ∕ddx, since the gas here is decelerating.

4.2 The thermodynamics of a J-shock plus cooling layer

We assume that the abundance of CO is fixed at $X\text{\_{\_CO}}=3\times {10}^{-4}$. In other words, we assume that essentially all the carbon is in CO, that the abundance of C is 1.5 × 10⁻⁴ (cf. Sembach et al. 2000), and that CO is not destroyed in the shock, nor has it had time to freeze out significantly (e.g. Goldsmith 2001).

Substituting from Eq. (9) into Eq. (26), we obtain $$X\text{\_{\_CO}}$$(27)

the second expression is obtained by substituting for ρ and a from Eqs. (22) and (25).

In Eq. (27), the leading term (outside the braces) represents deceleration due to cooling in the low-density optically thin regime. The term in braces represents the correction due to the cooling trending towards the high-density optically thick regime, and this can in principle cause the post-shock cooling to stall. Since the gas must decelerate to υ = ao (see discussion following Eq. (29) below), it follows that post-shock cooling will not stall due to the CO cooling becoming optically thick, provided that $$X\text{\_{\_CO}}$$(28)

where $n{\text{\_{\_H}}_{\text{2}}.\text{O}}=\rho \text{subO}/\overline{m}{\text{\_{\_H}}_2\text{}}$ is the number-density of molecular hydrogen in the pre-shock gas. If this inequality is well satisfied, we can neglect the term in braces, and this is generally the case. For example, if the pre-shock density is high, say $X\text{\_{\_CO}}$, and the temperature is low, say To = 10 K, we require υo ≫ 0.012 km s⁻¹, which is a rather weak constraint on υo. In most cases of interest, the pre-shock density is lower than 10⁵ cm⁻³, and the temperature cannot be much below ~10 K, in which case the minimum velocity is even lower. Furthermore, if the pre-shock density is any higher than 10⁵ cm⁻³, the gas is probably so well thermally coupled to the dust that molecular line cooling is redundant (e.g. Glover & Clark 2012b).

4.3 The thickness of the post-shock cooling layer

The thickness of the post-shock cooling layer is $$X\text{\_{\_CO}}$$(29)

The lower limit on this integral is the immediate post-shock velocity, given by Eq. (20). If we were to follow the gas until it cooled right back down to $T\text{subO}=\overline{m}\text{subO}a{\text{subO}}^2/k\text{\_{\_B}}$, the upper limit on the integral would be υ = ao²∕υo. However, by this stage the velocity would be highly subsonic and the LVG approximation would no longer be valid (in the sense that it would be seriously underestimating the cooling rate). Therefore we set the upper limit to υ = ao; in other words, we follow the deceleration until the velocity becomes subsonic. We believe this is justified, since the integral in Eq. (29) is dominated by the lower limit (υ = υo∕4).

If we substitute from Eq. (27) in Eq. (29), neglecting the correction term for optical thickness, we obtain $$X\text{\_{\_CO}}$$(30)

We have replaced ≃ with ≲ in Eq. (30) – and also in Eq. (32) below – because (a) we have neglected the contribution from the upper limit in the integral of Eq. (29), and this contribution can be quite large for relatively slow shocks, thereby reducing $\text{UpDelta}x\text{\_{\_PSC}}$ further; (b) for fast shocks there will be a significant additional cooling contribution from dust (Glover & Clark 2012b), and possibly also other molecules (Neufeld et al. 1995), due to the high post-shock density, reducing $\text{UpDelta}x\text{\_{\_PSC}}$ still further.

The gas passes through $\text{UpDelta}x\text{\_{\_PSC}}$ in a time $$X\text{\_{\_CO}}$$

and the integrated flux in CO cooling lines from the shock is $$X\text{\_{\_CO}}$$

4.4 Caveats on shock model

As already stated, gas-phase CO appears to be the dominant form of carbon at volume densities in the range $X\text{\_{\_CO}}$, and therefore our shock model can only be applied to shocks involving pre-shock gas with density, $n{\text{\_{\_H}}_{\text{2}}.\text{O}}$, in this range. Many shocks arising in turbulent molecular clouds will involve gas with density in this range.

J-shocks in turbulent molecular clouds will only approximate to being steady and planar if the counter-flows that form them have coherence lengths significantly larger than the thickness of the post-shock cooling layer, $\text{UpDelta}x\text{\_{\_PSC}}$. If we adopt coherence lengths from Larson’s relations, i.e. L ~ 1600 $X\text{\_{\_CO}}$ (Larson 1981), this condition is very well satisfied. For example, if we set $X\text{\_{\_CO}}=3\times {10}^{-4}$, we have $X\text{\_{\_CO}}$. This suggests that the assumption of a steady planar shock is reasonable.

If we assume that the approximate analytic formulation can be extended to T ~ 100 K, and that the pre-shock gas has T ~ 10 K, then our shock model can only be applied to shocks with pre-shock velocity υo ≲1.5 km s^-1, otherwise the immediate post-shock gas-kinetic temperature (i.e. pre cooling) is too high. This is a significant, but not critical, restriction, since there are likely to be many shocks satisfying this constraint.

For example, column-density maps of low-mass cores (Könyves et al. 2015) and filaments (Palmeirim et al. 2013) have centrally condensed profiles. However, the column-density contrast between the background and the line of sight through the centre of the core or the spine of the filament is small, typically ≲10. This indicates that the ram pressure of the gas accreting onto the core or filament is low, otherwise the profile would be flat. It also suggests that the background volume-density is not hugely lower than that inside the core or filament, hence the Mach number of the approximately isothermal accretion shock at the boundary of the core or filament must be low, $M\text{la}5$, and the pre-shock velocity must also be low, υo≲1 km s^-1. In addition we note that gravitational acceleration of the material accreting onto a marginally Jeans-unstable core or Ostriker-unstable filament can only generate marginally trans-sonic velocities, so gravitational acceleration does not alter this conclusion.

To first order, we can approximate $$X\text{\_{\_CO}}$$(35)

using Eq. (30). If again we set $X\text{\_{\_CO}}=3\times {10}^{-4}$, then the limits on $X\text{\_{\_CO}^}-1|\nabla \cdot \upsilon |$, i.e. $X\text{\_{\_CO}}$ yield $$X\text{\_{\_CO}}$$(36)

Since we already have the constraint υo≲1.5 km s⁻¹, this reduces to $X\text{\_{\_CO}}$, which is only slightly more restrictive than the upper limit on $n{\text{\_{\_H}}_2\text{}}$ derived at the start of this section.

Our shock model assumes that the carbon in the pre-shock gas is mainly in CO, and that the CO is not destroyed in the shock, nor does it freeze out. If this is not the case, cooling will be provided by other species, for example C⁺ (see Glover & Clark 2012a), but the gas will not cool so fast, and will probably not be able to get down to ~ 10 K, unless the post-shock density becomes high enough for the gas to couple thermally to the dust.

4.5 Discussion of shock results

The thickness of the post-shock cooling layer, $\text{UpDelta}x\text{\_{\_PSC}}$, and the post-shock cooling time, $\text{UpDelta}t\text{\_{\_PSC}}$, are both likely to be very small, as compared with other length- and time-scales in the molecular cloud. With typical pre-shock densities in a molecular cloud, $X\text{\_{\_CO}}$, and setting $X\text{\_{\_CO}}\simeq 3\text{}\times \text{}{10}^{-4}$, we have $X\text{\_{\_CO}}$ and $X\text{\_{\_CO}}$. These values are similar to the values derived by Pon et al. (2014, 2016), although they treated magnetically cushioned C-shocks, in which the compression and heating are less abrupt and less extreme.

One reason why CO is an effective post-shock coolant is that it has a low dipole moment, and therefore its lines remain optically thin up to quite high column-densities. A second reason is that the steady deceleration of the post-shock gas ensures that there is a broad range of velocities, and hence wavelengths, over which each line can be emitted, and therefore the lines do not readily become optically thick.

When treating the dynamics of turbulent interstellar clouds – for example, their assembly, cloud/cloud collisions, the formation of sheets and filaments and cores – Eqs. (30) and (32) can be used to assess rather quickly whether it is acceptable to treat the associated shocks as isothermal discontinuities, or alternatively, what resolution is needed to resolve the post-shock cooling layers. That is, provided that the inherent assumptions in our analysis are valid, viz. steady slow non-magnetic planar J-shocks, in which CO dominates the cooling and is not destroyed or removed from the gas-phase (as discussed in Sect. 4.4).

Pon et al. (2012) and Pon et al. (2016) estimate that the volume fraction of a turbulent molecular cloud that is involved in shock dissipation (i.e. the fraction that is in post-shock cooling layers) is $f\text{\_{\_PSC}}~0.001$. Given the short cooling times estimated above (Eq. (32)), it follows that a representative fluid element in such a cloud should be shocked at least once every Myr. Indeed, Pan & Padoan (2009) identify such shocks as an important overall heating mechanism for molecular clouds, although our estimates indicate that this heating should be very localised and transient.

5 Conclusions

We have derived a simple approximate analytic formulation for the cooling rate due to rotational transitions of the CO molecule, $$X\text{\_{\_CO}}$$(37)

An equivalent formulation, but giving $\text{UpLambda\_{\_CO}.\text{TOT}}/n\text{\_{\_CO}}$ in terms of $n{\text{\_{\_H}}_2\text{}}$ and T (instead of ρ and a), is given in Eqs. (15) through (18). Our approximate analytic formulation extends from the low-density optically thin limit, to the high-density optically thick limit, and includes the contributions from the different isotopologues of CO. It is based on physical considerations, but is calibrated against the detailed numerical results of GL78, and reproduces those results well (see Fig. 1). It should be usable for $$X\text{\_{\_CO}}$$

Using this formulation, we have derived estimates of the thickness of post-shock cooling layers, $\text{UpDelta}x\text{\_{\_PSC}}$, and post-shock coolingtimes, $\text{UpDelta}t\text{\_{\_PSC}}$, under the assumptions (i) that the shocks are steady slow non-magnetic planar J-shocks, and (ii) that CO dominates the cooling, and is not destroyed or removed from the gas-phase in the shock or the post-shock cooling layer. With these caveats, and given the constraints on shock velocity, pre-shockdensity and pre-shock temperature detailed in Sect. 4.4, our estimates can be used to justify treating shocks as isothermal discontinuities when the primary concern is dynamics on larger scales; to evaluate the resolution required to model a shock; and to estimate the integrated CO flux from a post-shock cooling region. The values of $\text{UpDelta}t\text{\_{\_PSC}}$ suggest that a typical fluid element in a molecular cloud should be shocked at least once every Myr.

Acknowledgements

We thank the anonymous referee for a very thoughtful and constructive report, which greatly improved the original version of this paper. APW gratefully acknowledges the support of a consolidated grant (ST/K00926/1), and SEJ a PhD studentship, both from the UK Science and Technology Funding Council. This work was performed using the computational facilities of the Advanced Research Computing at Cardiff (ARCCA) Division, Cardiff University.

References

All Figures

Fig. 1 Continuous curves give the total CO cooling rate, as predicted by the approximate analytic formulation derived here (Eqs. (15) through (18)). The different panels correspond to (a) T = 10K ; (b) T = 20K; (c) T = 40K; and (d) T = 60K (these temperatures are given in the top lefthand corner of each panel). The different curves represent different values of $X\text{\_{\_CO}}$ (these values are given down the righthand margin of each panel). The rates obtained by the detailed computations of GL78 are given at discrete values of the density, $X\text{\_{\_CO}}$; plus signs, stars, open circles, crosses and open squares correspond respectively to $X\text{\_{\_CO}}$.

In the text