CO in the Draco nebula: the atomic–molecular transition

ABSTRACT

This paper presents maps of the J=2-1 transition of CO towards the Draco nebula intermediate-velocity cloud (IVC). The maps cover 8500 square arcmin with a velocity resolution of 0.33 km s|$^{-1}$| and angular resolution of 38 |${\rm arcsec}$|⁠, or 0.11 pc at the cloud distance of 600 pc. The mapped area includes all the emission detected by the Herschel satellite with 250 |$\mu$|m intensity >5 MJy sr⁻¹. Previously published observations of the far-infrared (far-IR) emission and the 21 cm line of HI are used to derive the column density distribution of H|$_2$| and the abundance ratio CO/H|$_2$|⁠, as well as the distribution of the molecular fraction of hydrogen, which approaches 90 per cent over much of the brighter parts of the nebula. The CO emission is highly clumpy and closely resembles the structures seen in far-IR images. The kinematics of the CO show supersonic motions between clumps but near-thermal to trans-sonic motions within clumps, consistent with model predictions that the scale length for dissipation of supersonic turbulence should be |$\sim 0.1$| pc, mediated by kinematic viscosity and/or ambipolar diffusion. Different parts of the nebula show evidence for a spread of molecular formation time-scales of a few 10|$^5$| yr, comparable with the dynamical time-scale of the infalling gas. The IVC will likely merge with the Galactic interstellar medium in |$\sim 10^7$| yr, and the densest clumps may form an unbound cluster of low-mass stars.

1. INTRODUCTION

This paper presents new observations combined with previously published results to address one of the key questions about the processes that occur in the interstellar medium, namely how neutral atomic gas is converted to the cold, dense molecular clouds which are the necessary preconditions for star formation. It has long been recognized that the gas between the stars is not confined to a quiescent layer in the disc of the Milky Way. Early observational studies of the 21 cm line of neutral atomic hydrogen in both emission and absorption showed that there are two dominant regimes of temperature which exist in approximate pressure equilibrium. These ‘phases’ are generally referred to as the warm neutral medium (WNM) with a typical kinetic temperature of |$\sim 8000$| K, and the cold neutral medium (CNM) with temperatures of |$\sim 80$| K. The physical explanation for these components was first given by Field, Goldsmith & Habing (1969), who showed that the shape of the cooling function for the interstellar gas allowed the simultaneous existence of two stable temperature regimes in approximate pressure equilibrium. Further studies of thermal bremsstrahlung emission at X-ray wavelengths, as well as far-ultraviolet spectral lines from highly ionized elements revealed a third component, the diffuse hot intercloud medium (HIM) with temperatures of |$\sim 10^6$| K, that appears to be pervasive within the Galactic disc. The hot gas likely results from the collective effects of supernovae explosions and the winds of massive stars. Since the cooling time-scale for gas at this temperature is very long, the HIM can co-exist with the other phases for at least millions of years, as shown by McKee & Ostriker (1977).

An additional component of ionized gas at a temperature of |$\sim 8000$| K, the warm ionized medium (WIM), has also been identified by various observational diagnostics such as hydrogen recombination lines and pulsar dispersion measures; however, for the object considered in this paper, there appears to be very little WIM along the line of sight. Therefore we will concentrate on the role of the WNM, CNM, and the formation of dense molecular clouds from the atomic gas.

The HIM may be part of the story at least indirectly through the action of supernovae in the distant past which created the conditions for a ‘Galactic fountain’. This concept was first proposed to explain the existence of the high-velocity clouds (HVCs) detected in early HI 21 cm line surveys, which revealed discrete clouds of neutral atomic hydrogen with large, mostly negative radial velocities. The number and inferred masses implied that these clouds were falling on to the Galactic disc with an improbably large mass flux (Shapiro & Field 1976). An alternative explanation was that the infalling gas had previously been part of the Galactic disc gas and was blown into the Galactic halo by supernova explosions (Bregman 1980). Eventually some of the gas cooled and condensed as a result of thermal instability and was drawn back towards the disc by the gravitational potential of the stars. Further studies of the HI emission at high Galactic latitudes showed that the clouds with negative radial velocities (i.e. moving in towards the Galactic disc) tended to fall into two categories: those with velocities more negative than about |$-100$| km s|$^{-1}$|⁠, referred to as HVCs; and those with radial velocities in the range of roughly |$-100$| to |$-20$| km s|$^{-1}$|⁠, referred to as intermediate-velocity clouds (IVCs). HI emission features with radial velocities near zero were attributed to low-velocity neutral atomic clouds (‘LVC’) within the Galactic disc and participating in differential Galactic rotation.

The object of this study, the Draco nebula (hereafter, ‘Draco’), is considered an archetype of the IVC category. In Galactic coordinates the centre of the nebula is near |$l = 91$||${}^{\circ }$|⁠, |$b = 38$||${}^{\circ }$|⁠. The overall extent of the nebula is about |$4^{\circ } \times 4^{\circ }$|⁠. The distance to the IVC has been estimated by various authors; in this paper, we follow Miville-Deschênes et al. (2017b) and adopt a distance of 600 pc, which places the nebula about 370 pc above the Galactic plane, with an overall size of |$\sim 40$| pc. The Draco IVC is therefore within the inner Galactic halo, somewhat above the bulk of the WNM in the disc (Dickey & Lockman 1990). For a more extensive description of the properties as well as a chronology of observational studies of Draco, see the paper by Miville-Deschênes et al. (2017b) which also motivated our study of CO emission.

Early surveys of CO J=1-0 emission (Mebold et al. 1985; Rohlfs et al. 1989) showed that Draco had a significant amount of molecular gas for an IVC HI cloud, compared with other IVCs. This anomaly suggested that some fraction of the neutral atomic gas had been cooled and condensed sufficiently to promote the formation of molecules, but the details of this process were unclear. In observations with the Herschel far-infrared (far-IR) space telescope; however, Miville-Deschênes et al. (2017b) found that the thermal dust continuum (which they mapped in the 250 |$\mu$|m band with the SPIRE instrument), showed numerous narrow structures extending towards the Galactic plane. These ‘fingers’ of emission exhibited features down to the resolution limit of the image (18|${\rm arcsec}$|⁠) and some were probably unresolved. Miville-Deschênes et al. (2017b) interpreted these structures as condensations of the HI gas in the IVC resulting from a Rayleigh–Taylor (RT)-type instability as the cold, denser HI cloud encountered the WNM of the Galactic disc while falling back into the Galaxy as part of the ‘Galactic fountain’ scenario. If this picture were correct, then the 250 |$\mu$|m structures ought to coincide with CO emission, assuming that the cold, dense features marked the regions where the neutral atomic gas was turning into a molecular state. The Herschel observations, however, showed only the thermal dust continuum, not the presence of molecular gas. An earlier study by Herbstmeier, Heithausen & Mebold (1993) revealed concentrations of CO J=1-0 emission that did coincide with the brighter peaks of the Herschel 250 |$\mu$|m image, but the sensitivity and areal coverage of the CO map were limited to the brightest parts of the far-IR emission. We therefore undertook a more extensive mapping programme, described in the next section, to look for quantitative relations between the thermal dust continuum (as a tracer of total hydrogen column density), the atomic hydrogen content as traced by the 21 cm emission line, and the molecular gas traced by CO emission.

2. OBSERVATIONS

2.1 CO J=2-1 mapping

All the new observations of CO J=2-1 emission presented in this paper were obtained with the Heinrich Hertz Submillimeter Telescope¹ on Mount Graham, Arizona (elevation 3200 m). The observations used the on-the-fly (OTF) method to obtain spatially fully sampled maps of 10 |${\rm arcmin}$| square ‘tiles’, similar to other molecular cloud maps we have published previously (see e.g. Bieging & Kong 2022 and references therein for details of the mapping and calibration). The final maps have been lightly smoothed to a resolution of 38 |${\rm arcsec}$| (full width at half-maximum; FWHM), or 0.11 pc at the adopted distance of 600 pc. The spectral resolution was 250 kHz, which corresponds to 0.33 km s|$^{-1}$| at the line frequency of 230 GHz. At this resolution, the typical rms noise level in one velocity channel is 0.1 K brightness temperature. The total extent of the analysed spectral coverage was from −43 to −8 km s|$^{-1}$|(local standard of rest; LSR) at this resolution. We also obtained spectral line data at 1 MHz resolution that covered |$\sim 330$| km s|$^{-1}$|⁠, which allowed us to search for CO emission in the velocity range of both the LVC and the HVC located in the general direction of Draco, as shown in the HI maps of Miville-Deschênes et al. (2017b). No CO was detected in either the LVC or the HVC velocity range towards the main Draco features.

Our observations were motivated by the Herschel map presented by Miville-Deschênes et al. (2017b), which covers an area of about 4 degrees on a side. A spatially complete map of CO over such a large area would have been prohibitively time-consuming. From initial observations of some of the brightest far-IR features, however, it was clear that CO is only detectable where the 250 |$\mu$|m surface brightness exceeds about 5 MJy sr⁻¹. We therefore placed the mapping tiles to cover only regions that met this criterion. The CO maps follow the far-IR structures in an irregular pattern. In total, we mapped some 85 tiles or about 8500 square arcmin. Fig. 1a shows the Herschel map from Miville-Deschenes et al. (2017a), convolved to a resolution of 38 |${\rm arcsec}$| to match the final resolution of our CO maps. The labelled white rectangles identify subregions that will be considered in detail below; we use these designations for convenience to refer to specific areas of the nebula. Fig. 1b shows the integrated CO J=2-1 brightness temperature in K km s|$^{-1}$|⁠, for all of the mapped tiles. The zero level is offset into the dark blue of the colour palette, so where CO is not detected, the dark blue background illustrates both the extent of the mapped tiles, and the generally uniform noise levels achieved with these observations. The white contour shows the 5 MJy sr⁻¹ intensity level of the 250 |$\mu$|m emission from Fig. 1a. It is evident that the CO emission is entirely confined to the sightlines where the 250 |$\mu$|m surface brightness exceeds 5 MJy sr⁻¹. The CO emission at the top of the map extends slightly beyond the white contour, because the Herschel map is cut off at that position—see Fig. 1a—even though there is detectable CO beyond the far-IR map edge.

As with our previous studies of Galactic giant molecular clouds (GMCs), we also mapped the |$^{13}$|CO J=2-1 transition simultaneously. In contrast to typical GMCs, however, the |$^{13}$|CO line in Draco is mostly undetected except at the very brightest CO emission peaks. We will therefore not include the mapping data for the |$^{13}$|CO line in this paper. Where |$^{13}$|CO was detected, the ratio of line intensities (CO/|$^{13}$|CO) will be used to estimate the optical depth of the CO J=2-1 line.

2.2 Published observations used in the analysis

Critical to this study are previously published observations which we use to derive: the distribution of total column density of hydrogen in all forms; the distribution of neutral atomic hydrogen in the radial velocity range of the IVC; and from those, to infer the column density distribution of molecular hydrogen in the IVC.

To derive the distribution of total hydrogen in all forms, we use the Herschel 250 |$\mu$|m map downloaded from Miville-Deschenes et al. (2017a) and apply the simple proportionality between far-IR surface brightness and gas column density prescribed in Miville-Deschênes et al. (2017b). Their observations could not be used to derive the dust temperature distribution and from that, assuming a standard dust opacity law, the distribution of total gas column density. They argue, however, that the lower resolution data from the Planck all-sky maps show that the thermal dust emission over this region has a very uniform temperature. Then, assuming a standard dust/gas opacity law and the Planck dust temperature, they derive a linear proportionality which we use to convert their 250 |$\mu$|m map to a map of total hydrogen column density, at the resolution of the Herschel data, i.e. 18 |${\rm arcsec}$| (FWHM). The Herschel I(250) surface brightness map (in MJy sr⁻¹), convolved to match the resolution of our CO data (38 |${\rm arcsec}$|⁠) is shown in Fig. 1a.

We also considered two other, more recent analyses of the Herschel images to derive the dust temperature and opacity distributions (Schneider et al. 2022; Singh & Martin 2022). We concluded, however, that the uncertainties in critical parameters, namely the dust/gas mass ratio and the dust absorption coefficient, were so great that we opted to use the simple relation from Miville-Deschênes et al. (2017b) to derive the total hydrogen column density. The implications of the alternative dust models are discussed in Appendix  A.

The distribution of the neutral atomic hydrogen can be derived more directly than the dust column density from observations of the HI 21 cm emission line. Miville-Deschênes et al. (2017b) also show maps of the HI 21 cm emission over three velocity ranges, comprising the HVC, IVC, and LVC, from the Green Bank survey of HI emission at intermediate Galactic latitudes (Martin et al. 2015). The angular resolution of this survey is |$\sim 9.6$||${\rm arcmin}$|⁠, i.e. much lower than the Herschel map, but the separation in velocity shows that there is very little HI in the velocity range of the LVC towards the Draco emission at 250 |$\mu$|m. There are some moderately bright HI peaks in the HVC velocity range towards the centre of Draco, but the dominant IVC HI emission is clearly well-correlated with the Herschel 250 |$\mu$|m map. We examined our CO J=2-1 spectra at 1 MHz resolution over the velocity range of the HVC from some of our map tiles located near the HVC peak from Miville-Deschênes et al. (2017b; see their fig. 2, top panel). There was no CO detected at those positions in the velocity range of the HVC, to a brightness temperature limit of |$\sim 0.1$| K. (The non-detection of CO in the HVC velocity range may indicate that the HVC is extragalactic in origin and of low metallicity.)

A map of the HI 21 cm line over most of the Draco region was obtained with the DRAO hydrogen line interferometer array, which combined with the Green Bank survey, was used to produce the DHIGLS survey (Blagrave et al. 2017), at an angular resolution of |$\sim 1$||${\rm arcmin}$|⁠. This survey provides HI maps separated into the LVC and IVC velocity ranges. We have convolved the Herschel 250 |$\mu$|m image to match that of the DHIGLS survey (which covers most but not all of the Herschel field). Then, assuming that the Draco is dominated by the IVC gas, we subtract the DHIGLS HI column density map from the corresponding map of total hydrogen gas density to obtain a residual map which represents the distribution of molecular hydrogen, limited to the DHIGLS field and angular resolution.

3. RESULTS AND ANALYSIS

3.1 CO J=2-1 maps and velocity moments

A map of the intensity of the CO J=2-1 line integrated over the velocity range of the IVC is shown in Fig. 1b for the entire mapped region. As noted in Section 2.1, the |$^{13}$|CO line is mainly undetected over the entire field which implies that the CO J=2-1 towards Draco is optically thin over the whole nebula. Therefore the integrated intensity of CO should be a good indicator of the CO column density distribution. In Section 3.2, we use an analysis that assumes local thermodynamic equilibrium (LTE) excitation conditions to derive quantitative values for the CO total column density and compare those with the column densities and spatial distributions of the molecular and neutral atomic hydrogen.

In their study of Draco, Miville-Deschênes et al. (2017b) show that the structures in the far-IR thermal dust emission are likely the result of an R–T instability that develops as the infalling IVC atomic hydrogen cloud encounters the layer of WNM in the Galactic disc. Such a dynamical process would be expected to show velocity effects as the infalling gas is decelerated upon encountering the WNM, likely with variations over the |$\sim 40$| pc extent of the IVC. The Herschel map is limited to the thermal dust continuum, however, so no velocity information can be derived from those data. Our CO spectral line data cube has a resolution on the velocity axis of 0.33 km s|$^{-1}$|⁠, allowing us to look in some detail at the velocity structure of the R–T condensations. These appear to be the regions where the atomic gas is transforming to a predominantly molecular state, of which CO is a tracer. Fig. 1c shows the velocity centroid (i.e. the intensity-weighted mean velocity) over the entire mapped field. The velocities, referred to the LSR, range over some 12 km s|$^{-1}$|⁠, from |$-15$| to |$-27$| km s|$^{-1}$|⁠. The overall appearance of the velocity centroid map is not entirely chaotic; along the sinuous east–west filament between declinations +62|$^{\circ }$| 30 |${\rm arcmin}$| and +63|$^{\circ }$|⁠, there appear to be two sections with roughly constant centroid velocities, but with a sharp boundary between them across an apparent break in the CO emission, amounting to a difference of 3 km s|$^{-1}$| (⁠|$-27$| versus |$-24$| km s|$^{-1}$|⁠). The most typical velocity centroid value is −24 km s|$^{-1}$|⁠. The region extending along the south-east side of the nebula, which is where Miville-Deschênes et al. (2017b) find the best examples of R–T condensations, has generally less negative centroid velocities, possibly consistent with a deceleration of the infalling gas as it encounters the Galactic layer of WNM.

In Fig. 1d, we show the distribution of the second velocity moment or dispersion, |$\sigma _{\rm V}$|⁠, of the CO line. (If the line shape is a Gaussian, then the FWHM intensity equals |$2.355 \sigma _{\rm V}$|⁠.) The most common velocity dispersion is |$\sigma _{\rm V} \approx 0.35$| km s|$^{-1}$| as shown in a histogram of the distribution of |$\sigma _{\rm V}$|⁠. The sound speed for the gas at T = 20 K is about 0.3 km s|$^{-1}$|⁠, so the internal motions of the molecular gas are generally not much above thermal. There is a tail of values in the distribution extending up to |$\sim 3$| km s|$^{-1}$|⁠. Some of the largest line widths appear near the brightest CO emission peaks, particularly in the south-east and south. These areas are likely where the R–T instability is most dominant according to Miville-Deschênes et al. (2017b). The instability may lead to increased line broadening via turbulence and multiple density components along the line of sight. The narrowest lines are generally in regions located away from the R–T features identified by Miville-Deschênes et al. (2017b), i.e. the dark blue pixels in Fig. 1d. In the following paragraphs, we examine spectral channel maps (Figs 2–6) and line profiles (Fig. 7) in some of these R–T features.

To illustrate the complexity of the velocity field in the CO emission, we selected five regions labelled A–E, as shown in Figs 1(a) and (d), to display individual channel maps of the CO emission with spectral resolution of 0.33 km s|$^{-1}$| (which is about the same as the sound speed in the gas at T = 20 K). These maps are in Figs 2–6; every second or third channel is shown, as indicated in the figure captions. Regions C, D, and E, and the upper left corner of region B, correspond to the R–T instability structures shown in fig. 4 of Miville-Deschênes et al. (2017b). The CO channel maps have a spatial resolution of 38 |${\rm arcsec}$| while the original Herschel map in Miville-Deschênes et al. (2017b) has 18 |${\rm arcsec}$| resolution. In the channel map figures we show as white contours the Herschel 250 |$\mu$|m map smoothed to match the CO map resolution. Also shown are labelled crosses that mark the pixels where CO spectra are displayed in Fig. 7. (The crosses are only marked in the map of the channel nearest the emission peak of each spectrum.) The CO channel maps clearly show that there is a complex velocity structure within the R–T instability ‘fingers’, even when the far-IR map is smoothed to 38 |${\rm arcsec}$| resolution. Here, we comment briefly on the velocity structure in each of the five selected fields:

1. Field A (Fig. 2): This field is the southern-most feature in Draco. The Herschel 250 |$\mu$|m map has several clumpy structures as shown in the white contours, and the CO channel maps make clear that these clumps have a spread of velocities, with some of the far-IR peaks showing CO emission in only one or two of the 0.33 km s|$^{-1}$| channels. Thus, the CO line width in each of the clumps is narrower than the spread of velocities over the entire field.

2. Field B (Fig. 3): This field contains a complex set of far-IR structures, and together with Field A, dominates the Herschel map of Draco (e.g. Fig. 1a). The top-left corner of the field is included in fig. 4 of Miville-Deschênes et al. (2017b), which illustrates the features they identify as resulting from the R–T instability. The CO channel maps show that this gas has radial velocities between about |$-18$| and |$-22$| km s|$^{-1}$| (LSR), while the more prominent far-IR emission peaks in the east-central and west-central parts of the field have velocities between about |$-22$| and |$-26$| km s|$^{-1}$|⁠. As in Field A, the many clumps generally appear in only a few CO channels, implying that the individual components have relatively narrow CO line widths, on an overall larger spread in systemic velocity.

3. Field C (Fig. 4): In the high-resolution fig. 4 in Miville-Deschênes et al. (2017b), this feature shows numerous compact clumps of dust emission. As with Fields A and B, the CO maps for Field C show that the dust clumps have radial velocities from about |$-15$| to |$-24$| km s|$^{-1}$|⁠, but the individual components appear in only a few CO velocity channels. There is no obvious large-scale organization to the clump velocities, suggesting that this field has a fairly chaotic velocity field rather than a more organized overall dynamical structure.

4. Field D (Fig. 5): This feature is one of those resulting from the R–T instability according to Miville-Deschênes et al. (2017b). Their Herschel map shows a clumpy knot with a narrow extension to the south. Our CO channel maps show that this field, like the others, has many bright CO peaks which appear in only a few velocity channels, so each peak has a CO line width of order 1 km s|$^{-1}$|⁠, while the overall velocity range of the emission extends from |$-17$| to |$-26$| km s|$^{-1}$|⁠. The narrow ‘tail’ extending to the south of the main emission appears mainly in only one channel, so the molecular gas in the narrow ‘tail’ must have an internal velocity spread <1 km s|$^{-1}$|⁠.

5. Field E (Fig. 6): This feature is prominent in the Herschel map, with a clumpy, almost cometary appearance. The far-IR emission shows a bright area stretched along the south-west side and bright head consisting of a clumpy structure at the south-east end. The CO channel maps show a more organized velocity structure compared with the other fields, in that the upper right (north-west) corner has an LSR velocity around |$-23$| km s|$^{-1}$|⁠, with a gradient to the south-east. The brightest part of the far-IR image has a radial velocity in the range |$-25$| to |$-26$| km s|$^{-1}$|⁠. At the intermediate velocity around |$-24$| km s|$^{-1}$|⁠, there appear to be two extended features stretching from the south-east to the north-west, suggestive of gas being ablated and decelerated from the main clump as it plunges into the WNM layer. Such a scenario is admittedly somewhat speculative but is suggested by the interpretation of an R–T instability forming the observed structures. It would be of interest to see if numerical hydrodynamic simulations could reproduce both the density and velocity structures that are observed.

Spectra at the individual pixels marked in Figs 2–6 are displayed in Fig. 7. Many of the spectra have FWHM line widths of |$\sim 1$| km s|$^{-1}$| (e.g. A1, A2, A4, E2, and E4), while the broader lines appear to consist of two or more individual components seen superimposed along the line of sight (e.g. A3, E1, and E3). A detailed analysis of the entire Draco CO data cube for turbulence-related properties is beyond the scope of this paper but will be a subject for a future study.

3.2 LTE analysis

The physical conditions of the molecular gas traced by CO emission are not well-constrained, so we must make assumptions to derive quantities such as total CO column density and abundance relative to molecular hydrogen. In particular, a non-LTE approach using a Large Velocity Gradient treatment of radiative transfer would have required many poorly constrained assumptions. Instead we opted for a simple LTE analysis, following the method described in Kong et al. (2015), to estimate the total column density of CO along each line of sight. We assumed that the gas temperature was 20 K, i.e. the same as that used by Miville-Deschênes et al. (2017b) for the nebular dust temperature. We caution, however, that in the less dense parts of the nebula, the gas and dust may not be in thermal equilibrium so the gas kinetic temperature is not 20 K. Fortunately, if the CO rotational level populations are close to LTE conditions, the total column density we derive from the J=2-1 transition is relatively insensitive to the assumed temperature. For the observed transition, the line optical depth and the source function have roughly opposite temperature dependencies over the range of plausible values, i.e. between 10 and 60 K the derived total CO column density varies by only |$\pm 20~{{\ \rm per\ cent}}$|⁠. Thus, we believe that T|$_{{\rm gas}} = 20$| K is a reasonable guess for the LTE analysis. See e.g. appendix A in Ginsburg, Bally & Williams (2011) for a comparison of the temperature dependence of the emissivity for the three lowest rotational lines of CO.

Assuming LTE applies and using the method in Kong et al. (2015), we find that the maximum optical depth of the CO J=2-1 line is 0.65, and the average over all detected pixels is 0.14 assuming a kinetic temperature of 20 K. Therefore, the line is well described as optically thin, which should allow a robust calculation of total CO column density. This low optical depth is corroborated by the lack of |$^{13}$|CO J=2-1 emission over virtually all of Draco. If we assumed a lower kinetic temperature, e.g. as low as 10 K, the optical depths would increase by only about a factor of 2, hence still optically thin. An assumed temperature higher than 20 K would result in even lower optical depths.

The LTE map of total CO column density, N(CO), is shown in Fig. 8 (left panel). There is considerable variation over the many discrete features detected in CO. This variation is not simply a result of variation in the total gas density. Miville-Deschênes et al. (2017b) argue that the far-IR surface brightness from their Herschel 250 |$\mu$|m map is directly proportional to the total column density of gas, i.e. hydrogen in all forms. We use their relation to derive a map of total hydrogen column density, N(H|$_{{\rm total}}$|⁠) and then divide our map of N(CO) by N(H|$_{{\rm total}}$|⁠) to find the distribution of CO abundance relative to total hydrogen. The result is shown in Fig. 8 (right panel). The relative abundance of CO varies by at least an order of magnitude over the Draco components where CO is detected. One notable feature is the relatively high abundance peak seen in selected area C (see Fig. 1), where the ratio N(CO)/N(H|$_{{\rm total}}$|⁠) |$\approx 5\times 10^{-6}$| (CO molecules/H atom). More typically, this ratio is in the range |$5\times 10^{-7}$| to |$1\times 10^{-6}$| over most of Draco.

3.3 Derivation of molecular hydrogen column density and related quantities

We want to compare the calculated CO column density with that of molecular hydrogen, of which CO is the most commonly used tracer. The combination of the Herschel 250 |$\mu$|m dust continuum map and the atomic hydrogen maps from the DHIGLS survey can be used to infer the molecular hydrogen column density for the velocity range of the IVC, if we make certain reasonable assumptions. Miville-Deschênes et al. (2017b) argue that their far-IR continuum map is a good tracer for the total column density of hydrogen atoms in all forms over the mapped region and derive a simple linear relation between the 250 |$\mu$|m surface brightness, |$I(250)$|⁠, and the column density of hydrogen in all forms, N(H|$_{{\rm total}}$|⁠). We have opted to use their relation to derive a ‘baseline’ estimate of N(H|$_{{\rm total}}$|⁠). [Two alternative analyses of the Herschel dust continuum maps by Singh & Martin (2022) and Schneider et al. (2022) have reached somewhat disparate conclusions for the total dust and gas column density distributions, based on spectral energy distribution (SED) fitting of the far-IR data. In Appendix  A, we discuss the possible implications of those studies compared with our baseline assumption.) For this discussion, we use the simple relation

from Miville-Deschênes et al. (2017b).

The DHIGLS map cube of atomic hydrogen from Blagrave et al. (2017) can be divided into velocity ranges that cover the IVC and LVC emission, separately. Based on the lower resolution map of the GHIGLS survey, as shown in Miville-Deschênes et al. (2017b), we assume that the HI emission in the LVC can be neglected over the area of the dust structures in Draco. Likewise, the HVC feature in the GHIGLS map largely avoids the Draco dust features, and since that range of velocities was not included in the DHIGLS map, we neglect the HVC as well. Then after convolving the Herschel map to match the resolution of the DHIGLS image, which is 61|${\rm arcsec}$||$\times$| 54|${\rm arcsec}$| (Blagrave et al. 2017), we subtract the HI column density for the IVC from the N(H|$_{{\rm total}}$|⁠) image. The residual is assumed to represent the column density of molecular hydrogen, N(H|$_2$|⁠). This map is necessarily limited to the area covered in the DHIGLS image, which omits the top-most part of the area we mapped in CO. The resulting map of molecular hydrogen column density in Draco is shown in Fig. 9a.

To compare the CO and H|$_2$| column density maps, we smoothed and regridded the LTE-derived map of N(CO) (Fig. 8, left panel) to match the resolution and area of the map of N(H|$_2$|⁠) over the DHIGLS field. The ratio of these maps, N(CO)/N(H|$_2$|⁠), is shown in Fig. 9b. The ratio varies by at least one order of magnitude over the map, from |$\sim 1 \times 10^{-6}$| to |$1 \times 10^{-5}$|⁠, with the most typical values around |$3 \times 10^{-6}$|⁠. The highest values are in selected field D (see Fig. 1), which stands out compared with the bulk of the molecular gas in having a relatively high CO abundance relative to H|$_2$|⁠. Even so, this peak abundance is about one order of magnitude lower than the typical CO/H|$_2$| abundances inferred for most Galactic GMCs. In Section 4.2, we discuss the implications of these CO/H|$_2$| abundances and their spatial variations in terms of time-dependent models for molecule formation.

It is unlikely that the spatial variations in N(CO)/N(H|$_2$|⁠) could result from small-scale spatial variations in the dust properties within and between clumps. The maps of dust optical depth and dust temperature (Singh & Martin 2022; Martin, private communication) show relatively small variations at the resolution of 36 |${\rm arcsec}$|⁠.

Another quantity of interest is the so-called X-factor, which is sometimes used to estimate the total molecular column density, N(H|$_2$|⁠), directly from the CO brightness temperature integrated over the line profile. Most commonly this factor has been derived for the CO J=1-0 transition, simply because most observations of CO have previously been made of that line. Here we calculate the corresponding quantity X(CO 2–1) as the ratio of the molecular hydrogen column density, N(H|$_2$|⁠), divided by the integrated CO J=2-1 line intensity, I(CO 2–1). We smoothed and regridded the integrated CO intensity shown in Fig. 1b to match the image of N(H|$_2$|⁠) shown in Fig. 9a. The ratio of the two images is defined as the X(CO 2–1) factor, i.e. X(CO 2–1) |$\equiv$| N(H|$_2$|⁠)/I(CO 2–1). The resulting map is shown in Fig. 9c. It is evident that the X-factor varies by nearly an order of magnitude over Draco, from low values of order |$1 \times 10^{20}$| to |$1 \times 10^{21}$| (H|$_2$| molecules per cm²) (K km s|$^{-1}$|⁠)|$^{-1}$|⁠. The mean and median over the whole image are |$5.2 \times 10^{20}$| and |$4.5 \times 10^{20}$| (H|$_2$| molecules per cm²) (K km s|$^{-1}$|⁠)|$^{-1}$|⁠, respectively, with a long tail to higher values. Lewis et al. (2021) found a comparably large variation in X(CO 2–1) for the California molecular cloud (CMC), where they derived total gas column densities from Herschel observations with SED fits to the far-IR emission. They find values of X(CO 2–1) in a range from |$(0.5-52) \times 10^{20}$| (H|$_2$| molecules per cm²) (K km s|$^{-1}$|⁠)|$^{-1}$|⁠, with an average over the whole CMC of |$(4.1 \pm 0.8) \times 10^{20}$| (H|$_2$| molecules per cm²) (K km s|$^{-1}$|⁠)|$^{-1}$|⁠. Their cloud-averaged mean is only about 20 per cent lower than the mean we derive for Draco. For the CMC, however, they find a strong dependence of X(CO 2–1) on the dust temperature. For sightlines with T|$_{{\rm dust}} \ge 25$| K, they calculate a fairly small dispersion of values with an average X(CO 2–1) of |$(1.3 \pm 0.4) \times 10^{20}$| (H|$_2$| molecules per cm²) (K km s|$^{-1}$|⁠)|$^{-1}$|⁠. We point out that our selected region D has X(CO 2–1) close to the average of Lewis et al. (2021) for the higher T|$_{{\rm dust}}$| pixels in the CMC. Region D also shows the largest CO/H|$_2$| abundance ratio in Draco.

The average abundance ratios N(CO)/N(H2) reflect the evolutionary state of the cloud, in terms of molecule formation. The Draco molecular features must be very young compared with the CMC or other typical GMCs. We find a lower average CO abundance N(CO)/N(H|$_2$|⁠) for Draco compared with the CMC, even though the averaged X(CO)-factors are similar. Time-dependent numerical simulations, e.g. Gong, Ostriker & Kim (2018) show, however, that the X(CO) factor depends on a variety of cloud properties as well as cloud age, including CO optical depth, gas temperature, incident radiation, cloud structure, etc. The similarity of the average X-factors between Draco and the CMC is probably related to these other properties, and simply reinforces the conclusion that the CO X-factor is a poor measure of total gas mass on sub-parsec scales.

An earlier study of CO in parts of the Draco by Herbstmeier et al. (1993) used observations of the CO J=1-0 transition, together with dust continuum maps by IRAS and HI 21 cm maps with various radio telescopes, to determine total gas column densities and also the CO X-factor. They derived X-factors that were an order of magnitude smaller than those we find here using the CO J=2-1 line. This discrepancy is unlikely to be entirely a result of the different CO lines used. Herbstmeier et al. (1993) used the IRAS maps with angular resolution of 3 |${\rm arcmin}$||$\times$| 5 |${\rm arcmin}$|⁠, but the Herschel 250 |$\mu$|m images of Miville-Deschênes et al. (2017b) show that the IVC features with detected CO emission are |$\sim$|1 |${\rm arcmin}$| or less in size. Consequently their measures of total column density for the IVC features are severely beam-diluted. The telescopes they employed for CO mapping had resolutions of |$\sim 4$||${\rm arcmin}$| with spectra sampled on a 4 |${\rm arcmin}$| grid, so the effective resolution was |$\sim 8$||${\rm arcmin}$|⁠. We suspect that these differences in map resolution are the principal reason for the discrepancy in X(CO) from Herbstmeier et al. (1993) compared with our values for X(CO 2–1).

In the context of molecule formation, we want to know the fraction of total hydrogen that has been converted to molecular form. We can only calculate a ratio of the column density-averaged fraction, 2N(H|$_2$|⁠)/N(H|$_{{\rm total}}$|⁠) which is the line-of-sight averaged fraction of H atoms in molecular form. The resulting map is shown in Fig. 9d. The image shows that a significant fraction of the Draco cloud has been converted to molecular hydrogen especially in the ‘leading edge’ of the filamentary structures seen in the Herschel 250 |$\mu$|m map, and in the long concentration of emission on the southern part of the nebula. In these regions the molecular fraction exceeds 90 per cent. Other parts of the filaments have molecular fractions on the order of 50 per cent to 70 per cent. Fig. 9 d shows that in fact a significant portion of Draco is molecular, contrary to what is sometimes assumed in recent studies (e.g. Schneider et al. 2022).

4. DISCUSSION

4.1 Uncertainties in derivation of column densities

Our analysis as summarized in Fig. 9 is based on three components of the gas content, all constrained to be the line-of-sight integrated column densities: (1) the CO molecule, derived from our observations of the J=2-1 transition, where the total CO column density N(CO) is obtained from an assumed LTE level population at a kinetic temperature of 20 K; (2) atomic hydrogen, where we rely on the DGHILS maps of Blagrave et al. (2017) of the 21 cm HI line as a robust measure of N(HI); and (3) the total gas column density of hydrogen, N(H|$_{{\rm total}}$|⁠), which we obtain from the Herschel 250 micron intensity map of Miville-Deschenes et al. (2017a) and apply the simple proportionality of Miville-Deschênes et al. (2017b).

If our assumption of LTE level populations for CO is correct, then the derived total CO column density is relatively insensitive to the assumed kinetic temperature. Moreover, the lack of detected |$^{13}$|CO and where detected, low brightness, confirms that the CO J=2-1 line is generally optically thin across the nebula. We assume that T|$_{{\rm kin}} = 20$| K, but for values in the range |$14{-}25$| K the LTE total column density varies by only a few per cent. Since the J=2 level lies at |$\sim 16$| K above ground, the levels J=1 and 2 should be nearly thermalized provided the gas density is high enough. If observations of the J=1-0 and 3-2 lines with sensitivity and resolution comparable with our maps were available, a comparison could confirm whether our assumption of LTE is correct. Absent such data, we regard our approach as the least subject to other uncertain quantities that would affect, e.g. a Large Velocity Gradient analysis.

The column density of atomic hydrogen is very well determined by the integrated intensity of the 21 cm HI line, since the emission is nearly independent of the gas temperature. We take the DHIGLS map of the IVC gas from Blagrave et al. (2017) as a firm lower limit on the total gas column density. The derived N(HI) values are somewhat compromised, however, by the modest resolution (⁠|$\sim 1\ {\rm arcmin}$|⁠) and more seriously, the limited sensitivity of their data. Even so, the HI map of the IVC emission is a critical constraint on the total gas column density for the area of Draco covered by the DHIGLS field.

We therefore argue that the maps of both N(CO) and N(HI) are reasonably secure. A more serious concern is the determination of the total gas column density as inferred from various analyses of the far-IR thermal dust emission measured by the Herschel PACS and SPIRE cameras. In this paper, we will adopt as a ‘baseline’ the simple proportionality from Miville-Deschênes et al. (2017b), namely that N(H|$_{{\rm total}}) = 2.49 \times 10^{20}~ I(250)$| H atoms per cm² where |$I(250)$| is the zero-corrected intensity in MJy sr⁻¹ at wavelength 250 |$\mu$|m from the SPIRE instrument. Since the publication of Miville-Deschênes et al. (2017b) there have been two other analyses of the Herschel maps of Draco which attempt a more refined treatment of the thermal dust emission and from that, to infer the distribution of total gas column density. In Appendix  A, we consider the implications of these studies with respect to our ‘baseline’ map described in Section 3.2. In the following discussion, however, we will apply the simple proportionality of Miville-Deschênes et al. (2017b) for the calculation of total hydrogen column density.

4.2 Time-scales for molecule formation

A significant objective of our Draco observations is to provide constraints on models of molecule formation in interstellar clouds. Constructing reasonably plausible simulations of this process is an exceedingly complex challenge for numerical computations, so most published results to date are restricted to simplified versions of either cloud structure and dynamics and/or of relevant chemical processes (typically including ion–molecule, neutral radical, and grain surface reactions, as well as photo-processes). Well-resolved time-dependent models in particular are rare. Moreover, the Draco cloud is substantially different from typical molecular clouds in the Galactic disc, as we noted above. Despite these caveats, the properties derived from our CO observations and related published work offer some useful data applicable to numerical simulations of molecule formation.

The first constraint we found is that CO only becomes detectable for a Herschel I(250 |$\mu$|m) |$\ge 5$| MJy sr⁻¹, which corresponds to a minimum total column density of H in all forms of |$\sim 1.3\times 10^{21}$| H atoms per cm|$^2$| for the onset of CO formation, if we assume the conversion factor in Miville-Deschênes et al. (2017b). Above this intensity threshold, there is only a fair correlation between I(250 |$\mu$|m) and the integrated CO J=2-1 line intensity when plotted over all the detected pixels in the entire field mapped in CO. As shown in Figs 2–6, however, the remarkably complex velocity structure of the molecular gas is not revealed in the far-IR continuum. The inability of using the thermal dust continuum to determine total gas column density as a function of radial velocity is perhaps the most important limitation on deriving relative molecular abundance distributions, e.g. CO/H|$_2$|⁠, for comparison with numerical models. The substantial uncertainties in the assumed gas/dust mass ratio and the dust mass opacity coefficient appropriate for Draco add additional potential sources of error in quantitative comparisons.

Despite these limitations, we found that correlation scatter diagrams of our derived column densities of CO versus molecular hydrogen show clear trends, when restricted to the smaller regions A–E shown in Fig. 1. These comparisons are illustrated in Fig. 10. In each region, there are one or more distinct ‘streams’ of points which are nearly linear with different slopes but intersecting the x-axis at about the same value of N(H|$_2$|⁠). The coloured lines show four examples with slopes as labelled. In each panel, all four lines have the same x-intercept, which differ from panel to panel but have a mean of N(H|$_2$|⁠) |$= 4 \pm 1 \times 10^{20}$| H|$_2$| molecules per cm². (These examples are by-eye estimates, not linear regressions, given the large scatter.) The x-intercepts imply that the CO J=2-1 line is detected only where the column density of molecular hydrogen corresponds to a Herschel I(250) |$\ge$| 3 MJy sr⁻¹, assuming all hydrogen is molecular and using the scale factor from Miville-Deschênes et al. (2017b). This lower limit is consistent with the value noted above of 5 MJy sr⁻¹ if, as is likely, some hydrogen is in atomic form at the outer edges of the molecular gas.

The approximately linear appearance of the ‘streams’ in Fig. 10 implies also that the abundance ratio of column densities, N(CO)/N(H|$_2$|⁠) rises from zero at the x-intercept value to an asymptotic limit given by the slope of each line. It is striking that three of the five selected fields, namely C, D, and E, appear to have one dominant feature, but with clearly different slopes between the three fields, while two fields, A and B, show two or three distinct ‘streams’ each with different slopes. Fields A and B are spatially larger and have more structure in velocity than the other three fields, suggesting that the multiple streams of points are associated with distinct spatial/kinematic components. The differences in slopes imply differences in the current molecular abundance ratio of CO/H|$_2$| between and within the five selected fields. The high-resolution time-dependent chemical modelling of Glover et al. (2010) suggests that the explanation for the different slopes lies in the time since molecule formation began in each feature. Their chemical evolution models for formation of CO (see their Fig. 5) show a steep initial rise in the mass-weighted CO/H|$_2$| abundance, increasing from |$10^{-6}$| to |$10^{-5}$| in only about 300 000 yr.

For comparison, we can estimate the dynamical time-scales for the individual features in Draco from our measured CO radial velocities as shown in Fig. 1. A typical CO line has a radial velocity of |$-20$| km s|$^{-1}$|⁠. If we assume that the spatial motion is directly towards the Galactic plane, at the latitude of 38|$^{\circ }$| the total spatial velocity is |$-33$| km s|$^{-1}$| and the transverse velocity is 26 km s|$^{-1}$|⁠. At our adopted distance of 600 pc, the transverse motion of a feature would be about 15 |${\rm arcmin}$| in |$10^5$| yr. This angular size is comparable with the dimensions of the selected fields illustrated earlier, e.g. Field D (Fig. 5). The multiple streams of points seen in Fields A and B (Fig. 10) may indicate that in these spatially larger fields, individual clumps of gas have undergone formation of CO with different time-scales which correspond to the dynamical time-scale for the infalling gas. In this picture, the streams with the steepest slopes (e.g. Field D) have been forming CO for longer times than those with shallower slopes (e.g. Field E). Since we do not know the detailed history of the gas motions in Draco, we cannot make specific conclusions about molecular (i.e. CO) formation in the individual features in our maps. Our observations do appear to be consistent with the time-dependent models of Glover et al. (2010), who find that molecule formation is rapid once gas densities reach high enough values, and that the CO/H|$_2$| abundance varies by large factors within a clumpy turbulent cloud.

4.3 Dynamics

An important motivation for our spectroscopic mapping is the capability to distinguish kinematic features in the radial velocity dimension, given the very high spectral resolution factor |$R = \nu _o/\delta \nu = (230 \rm {~GHz}/250 \rm {~kHz}) \approx 10^6$| of the data. This capability is shown in Figs 2–6, which demonstrate that the structures found in the Herschel far-IR continuum maps typically contain individual molecular components with relatively narrow radial velocity widths, but which span a much broader range of velocities over the entire structure. These dynamical properties are also evident in the corresponding spectra shown in Fig. 7. The first moment or velocity centroid maps of the five selected fields A–E span between 4 and 8 km s|$^{-1}$|⁠, while the corresponding line-widths (the second moment or velocity dispersion) at each pixel typically show a single peak at 0.5 km s|$^{-1}$| and extend to less than 2 km s|$^{-1}$|⁠. This kinematic behaviour is seen over the entire region mapped in CO, as shown in the lower two panels of Fig. 1. The many individual clumps (⁠|$\sim 5000$|⁠) described by Miville-Deschênes et al. (2017b) from their 250 |$\mu$|m image evidently have internal velocity motions close to the thermal width for gas temperatures of order 20 K, but these clumps have significantly larger (i.e. supersonic) systemic velocities relative to each other. This important kinematic information cannot be extracted from the thermal dust continuum emission, but is clearly a critical constraint on numerical simulations of the dynamics of the Draco cloud.

The narrow CO line-widths of individual gas clumps do appear to be consistent with the estimates by Miville-Deschênes et al. (2017b) for the viscosity and corresponding turbulent energy dissipation scales in Draco. They found that kinematic viscosity must be dominant compared with molecular viscosity, and from that they estimate a typical length of 0.1 pc for dissipation of supersonic turbulence. (They also estimate a comparable scale length of 0.1 pc for energy dissipation by ion-neutral friction, i.e. ambipolar diffusion, assuming reasonable values for the gas density and magnetic field strength.) The infalling IVC cloud has a Mach number of |$\sim 6$|⁠, according to Miville-Deschênes et al. (2017b). This value is comparable with the spread of velocities in the clumpy structures of the individual R–T features highlighted in selected fields A–E (Figs 2–6, also Fig. 7). We observe that the individual molecular clumps have CO line-widths close to thermal or slightly trans-sonic, where the spatial resolution of our CO map is 0.11 pc. If the turbulent dissipation scale were much less than 0.1 pc, we would expect to find CO line-widths well above the thermal value, since each resolution element in our maps should encompass a range of turbulent velocities along the line of sight. In this picture, the velocity structure of the R–T features was determined by supersonic turbulence as the IVC encountered the WNM layer, while the internal velocity spread of the individual molecular clumps has been reduced to near-thermal widths by energy dissipation mediated via kinematic viscosity and possibly ambipolar diffusion. If the scale length for dissipation of turbulence is in fact about 0.1 pc, we would predict that CO J=2-1 observations with higher angular resolution than ours, e.g. |$\sim$|10 |${\rm arcsec}$|⁠, would reveal some clumps with even narrower line-widths. The time-scale for turbulence dissipation within a molecular clump is evidently comparable with the chemical formation time-scale, i.e. on the order of a few times |$10^5$| yr, but with some variation between clumps as implied by the variation in slopes seen in Fig. 10.

We note that a numerical simulation by Saury et al. (2014) showed that compressive turbulence in the WNM can lead to rapid cooling and formation of dense structures within the WNM, with time-scales |$\lt 10^6$| yr. This simulation was not precisely comparable with the scenario proposed to explain the morphology of Draco, but does suggest that supersonic turbulence created as the infalling IVC encounters the WNM layer in the Galactic plane can indeed lead to the conditions needed for rapid formation of CO molecules as found by Glover et al. (2010). Our observations of CO emission could provide constraints on numerical simulations which combine the dynamics of the formation of R–T structures by infalling gas, together with the chemistry of H|$_2$| and CO.

4.4 Future evolution of Draco

We can speculate about the future evolution of the Draco IVC based on the available observations and analysis. If the nebula is now 370 pc above the Galactic plane and falling in with a velocity of |$\sim 40$| km s|$^{-1}$|⁠, and continues to move at that speed, it will take about 10 million yr to reach the mid-plane. If the gas is slowed as it encounters more of the WNM, as seems likely, the time to reach the minimum in the stellar gravitational potential would be longer. The denser clumps we see in the R–T structures may be less affected by such slowing than the more diffuse parts of the IVC, which would lead to some vertical stretching or separation of the IVC material. The denser parts might even overshoot somewhat before being drawn back by the gravity of the stellar disc, in a damped oscillatory motion, eventually settling into the Galactic mid-plane in a few tens of millions of years, i.e. short compared with the Galactic orbital time-scale for the Solar neighbourhood.

Miville-Deschênes et al. (2017b) used the clumpfind algorithm to identify about 5000 distinct condensations in their Herschel 250 |$\mu$|m map. They calculated a total mass for these condensations of |$\sim 5000$| M|$_{\odot }$|⁠, assuming the standard gas/dust mass ratio and dust opacity as discussed in Section 4.1 above. They also derived masses for the individual clumps in the range 0.1–20 M|$_{\odot }$|⁠, with a median mass of 0.53 M|$_{\odot }$|⁠. The great majority of these condensations were not gravitationally bound, however, with only 15 per cent of the total mass in bound structures. Provided these clumps remain together as an identifiable assembly, the future appearance of Draco may be something like the dark clouds we see in the Solar neighbourhood, at the lower end of the mass spectrum of Galactic molecular clouds, e.g. table 32.2 in Draine (2011). The bulk of the infalling gas will likely disperse into the Galactic interstellar medium on a time-scale of a few times 10|$^7$| yr. The surviving more massive molecular condensations may eventually form low-mass stars in a sparse unbound stellar cluster.

5. SUMMARY

We present new observations made with the Steward Observatory Heinrich Hertz Submillimeter Telescope, of the J=2-1 emission from CO towards the Draco IVC. The map covers about 8500 square arcmin, with an angular resolution of 38 |${\rm arcsec}$| and spectral resolution of 0.33 km s|$^{-1}$| over the velocity range of the IVC. The area mapped was selected to cover all of the thermal dust emission seen in the 250 |$\mu$|m Herschel image published by Miville-Deschênes et al. (2017b) with intensity I(250) |$\ge$| 5 MJy sr⁻¹, which marks the lower limit at which CO is detected. From the relationship in Miville-Deschênes et al. (2017b) between I(250) and total hydrogen column density, N(H|$_{{\rm total}}$|⁠) as a ‘baseline’, we derive a map of molecular hydrogen column density, N(H|$_2$|⁠) by subtracting from N(H|$_{{\rm total}}$|⁠) the atomic hydrogen column density in the IVC mapped in the DHIGLS survey (Blagrave et al. 2017). The column density of CO molecules is obtained from an LTE analysis, under the assumption that the CO J=2-1 line is optically thin over the nebula. This assumption is justified by the observed lack of |$^{13}$|CO emission except in a few positions. We emphasize also that the derived total CO column density, N(CO), is insensitive to the assumed gas kinetic temperature for the plausible range of values. From these column density maps, we present maps of the molecular abundance ratio N(CO)/N(H|$_2$|⁠).

In the Appendix, we also discuss two recent alternative analyses (Schneider et al. 2022; Singh & Martin 2022) of the thermal dust continuum maps by the Herschel and Planck satellites, in terms of the derived total gas column density. To reconcile these analyses, we find that our ‘baseline’ N(H|$_{{\rm total}}$|⁠) map would have to be scaled down by a factor of 0.3, assuming that their adopted gas-to-dust ratios and/or dust mass opacity coefficients are not correct for Draco compared with matter in the Galactic plane.

The Draco IVC shows a remarkably clumpy structure, with ‘fingers’ of dust emission extending towards the Galactic plane. These features have been interpreted by Miville-Deschênes et al. (2017b) as arising from an R–T instability that develops when the in-falling atomic gas in the IVC encounters the WNM at supersonic velocities. The high resolution in velocity of our CO map spectra allows a detailed examination of the kinematics of this process, and of the formation of CO in the in-falling gas. Our main conclusions are:

1. CO only becomes detectable for a Herschel I(250 |$\mu$|m) |$\ge 5$| MJy sr⁻¹, which corresponds to a minimum total column density of H in all forms of |$\sim 1.3\times 10^{21}$| H atoms per cm|$^2$| for the onset of CO formation, if we assume the conversion factor in Miville-Deschênes et al. (2017b). The derived CO column density does not depend on the total hydrogen column density, so if we scale down N(H|$_{{\rm total}}$|⁠) by a factor of 0.3 as suggested in the Appendix, the onset of CO formation would be where N(H|$_{{\rm total}}$|⁠) exceeds about |$4\times 10^{20}$| H atoms per cm|$^2$|⁠.

2. The ratio of CO to molecular hydrogen column densities, N(CO)/N(H|$_2$|⁠), varies by at least an order of magnitude across the nebula, ranging from |$10^{-6}$| to |$10^{-5}$| as shown in Fig. 9b. Even the highest ratio is about a factor of 10 less than is typically found in the Galactic GMCs. The fraction of hydrogen in molecular form (Fig. 9d) shows a gradient along the R–T ‘fingers’ with values of |$\sim 90~{{\ \rm per\ cent}}$| at the highest towards the SE and falling to |$\sim 30~{{\ \rm per\ cent}}$| away from the peaks. The large southern concentration and the sinuous filament in the northern part of Draco, however, have a high molecular fraction over most of their extent.

3. Correlation diagrams of individual pixel values for N(CO) versus N(H|$_2$|⁠), when plotted for isolated map fields as in Fig. 10, show well-defined linear features (‘streams’) with a range of slopes but similar x-axis intercepts. We interpret these streams as indicators of CO formation that is occurring in the dense gas clumps but with a spread in time-scales since the onset of chemical evolution. The inferred abundance ratios, N(CO)/N(H|$_2$|⁠), are consistent with the evolutionary models of Glover et al. (2010) with ages spanning a few times 10|$^5$| yr.

4. The very high spectral resolution of our data (⁠|$R \approx 10^6$|⁠) allows a detailed kinematic analysis of the CO emitting gas. Expressed as velocity moments (centroid and dispersion) we find that the spread in centroid velocities is large compared with the dispersion of spectra at individual map pixels. The spread in centroids over fields of |$\sim 15$||${\rm arcmin}$| in size is comparable with the predicted sonic Mach number of the R–T structures (Miville-Deschênes et al. 2017b) with |$M_{\rm sonic}\approx$| 6. The velocity dispersions show a peak at |$\sim$|0.5 km s|$^{-1}$|⁠, which is near the thermal line width for a kinetic temperature of 20 K, or at most is slightly trans-sonic. We argue that these values for the distribution of velocity centroids and dispersions are consistent with turbulent motions driven as the IVC gas falls in and encounters the WNM. The scale length for dissipation of turbulence is |$\sim$|0.1 pc, as predicted by Miville-Deschênes et al. (2017b) in their model for formation of R–T instabilities. The energy dissipation must then be mediated by kinematic viscosity and/or ambipolar diffusion, rather than molecular viscosity.

5. The Draco should settle on to the Galactic mid-plane in |$\sim 10^7$| yr; it will then resemble a dark cloud such as those seen in the Solar neighbourhood. The densest clumps may eventually form low mass stars in an unbound cluster, while the bulk of the gas will merge with the general interstellar medium.

ACKNOWLEDGEMENTS

The Heinrich Hertz Submillimeter Telescope is operated by the Arizona Radio Observatory, which is a unit of Steward Observatory at the University of Arizona. We thank P.G. Martin, M.-A. Miville-Deschenes, V. Ossenkopf-Okada, and N. Schneider for useful comments and for providing maps in digital form.

DATA AVAILABILITY

The Herschel 250 |$\mu$|m map used in this study is available at Miville-Deschenes et al. (2017a). The 21 cm HI data from the DHIGLS survey (Blagrave et al. 2017) are available at the website https://www.cita.utoronto.ca/DHIGLS/. The CO J=2-1 map cube is available in the Harvard Dataverse (https://dataverse.harvard.edu) with DOI https://doi.org/10.7910/DVN/7DO6EP.

Footnotes

REFERENCES

APPENDIX A: COMMENTS ON ALTERNATIVE ANALYSES OF THERMAL DUST CONTINUUM

Since the publication of Miville-Deschênes et al. (2017b) there have been two other analyses of the Herschel maps of Draco which attempt a more refined treatment of the thermal dust emission and from that, to infer the distribution of total gas column density. Here, we consider the implications of these studies as regards our ‘baseline’ map described in Section 2.2.

Singh & Martin (2022) made an extensive analysis of Herschel PACS/SPIRE maps of a number of molecular clouds to derive the distributions of dust temperature and optical depth in their Herschel Optimized Tau and Temperature (HOTT) survey, with data available at www.cita.utoronto.ca/HOTT. They made modified blackbody fits using the standard parametrization:

where the dust optical depth at a fiducial frequency, |$\tau _{\nu _o}$|⁠, is given by

Here r is the dust-to-gas ratio by mass, and |$\mathrm{ N(H)}$| is the column density of hydrogen atoms in all forms, |$\mu$| is the mean a.m.u. mass per gas particle, and |$m_\mathrm{ H}$| is the mass of the hydrogen atom. The dust mass opacity coefficient, |$\kappa (\nu)$| is assumed to have a power-law frequency dependence. Singh & Martin (2022) find the power-law exponent |$\beta$| by interpolating the lower resolution maps of |$\beta$| from the Planck XI dust model (Planck Collaboration XI 2014). Free parameters that were fitted to the intensity |$I_\nu$| in the 160 |$\mu$|m PACS and 250, 350, and 500 |$\mu$|m SPIRE images at each pixel were the P lanck dust temperature, |$T_{{\rm dust}}$| and the dust optical depth. Though their published results concentrated on the clouds in the Gould’s Belt Survey (André et al. 2010) and did not explicitly include the Draco field, that region was also analysed and kindly provided to us by P. G. Martin (private communication), in the form of maps of |$\tau _{1~\rm {THz}}$| and |$T_{{\rm dust}}$|⁠, as well as a map of total hydrogen column density, N(H|$_{{\rm total}}$|⁠), which assumes nominal values of the dust/gas mass ratio |$r = 1/100$|⁠, and the dust mass-opacity coefficient, |$\kappa _{1~\rm {THz}} = 10$| cm|$^2$| g|$^{-1}$|⁠. We convolved and regridded their map of N(H|$_{{\rm total}}$|⁠) to match the map of HI in the IVC from the DHIGLS survey, and examined the ratio of the two images. If the DHIGLS map is a firm lower limit to the total gas column density, as we assume, then the ratio N(H|$_{{\rm total}}$|⁠)/N(HI) must be greater than 1 everywhere. Applying the nominal values for gas/dust mass ratio and |$\kappa _{1~\rm {THz}}$| we find that many pixels have a ratio <1, implying that either one or both of the nominal values are incorrect for the Draco dust and gas. Singh & Martin (2022) found that the dust temperatures in their modified blackbody fits were typically in the range of 19–20 K and ranged down to 17 K for the highest opacity pixels. Thus, even at the higher angular resolution of Herschel compared with Planck, they found that the dust temperature is relatively uniform over the Draco field.

The DHIGLS N(HI) map does have a significant amount of scatter as a result of noise. We found that, after regridding the Singh & Martin (2022) map of N(H|$_{{\rm total}}$|⁠) to match the DHIGLS N(HI) map, a plot of individual pixel values of the ratio N(H|$_{{\rm total}}$|⁠)/N(HI) on the y-axis versus the integrated CO J=2-1 line intensity (a directly observed quantity) showed a clear trend such that the ratio N(H|$_{{\rm total}}$|⁠)/N(HI) increases with CO line intensity consistent with the assumption that CO traces H|$_2$|⁠. Given the condition that N(H|$_{{\rm total}}$|⁠)/N(HI) |$\ge 1$|⁠, the trend line must intersect the y-axis at a value of at least 1, where all the gas is atomic. Values greater than 1 imply a fraction of the gas has become molecular, as expected by the detection of CO. For the selected regions, we found that the trend lines intersect the (ratio) y-axis at values in the range of 0.3–0.7, i.e. well below the expected minimum value of 1. This discrepancy suggests that the nominal values assumed for |$\kappa _{\nu _o}$| or the dust/gas mass ratio or both, are not correct for the dust and gas in Draco. The nominal values used here are likely more appropriate for molecular clouds in the Galactic plane. If the matter in Draco was ejected from the plane by some combination of stellar winds and supernovae and is now falling back, i.e. the ‘Galactic fountain’ scenario, it is likely that the dust has been altered in the process to have properties rather different from its original state. If icy mantles were lost, e.g. the dust mass absorption coefficient |$\kappa _{\nu _o}$| may be changed from the assumed 10 cm|$^2$| g|$^{-1}$|⁠. If some fraction of the dust mass were destroyed then the assumed dust/gas mass ratio of 1/100 may well be too high (Draine & Salpeter 1979; Seab & Shull 1983; Jones, Tielens & Hollenbach 1996; Slavin, Jones & Tielens 2004). (In a recent study of 94 high latitude quasars, Shull & Panopoulou (2024) find that the average dust/gas ratio is a factor of |$\sim 1.5$| lower than the standard value for gas in the Galactic plane, with some sightlines showing values substantially lower than the average.) A combination of such effects could easily combine to make an increase by the factor of |$\sim 3$| necessary to bring the HOTT N(H|$_{{\rm total}}$|⁠) into agreement with the minimum value set by the DHIGLS N(HI) map for the IVC.

A study by Schneider et al. (2022) also analysed the same Herschel observations of the Draco field and reached somewhat different conclusions. They made fits to the SEDs using the 160 |$\mu$|m PACS and 250, 350, and 500 |$\mu$|m SPIRE images, all convolved to the lowest resolution (36.7 |${\rm arcsec}$|⁠). In contrast to Singh & Martin (2022), Schneider et al. (2022) fixed the frequency-dependent power-law index |$\beta$| at 2.0 and set the dust mass opacity coefficient to |$\kappa _{\nu } = (0.13/r) ({\nu }/{1.2 {\rm ~THz}})^\beta$| cm|$^2$| g|$^{-1}$|⁠, following Juvela et al. (2015) for |$\lambda = 250$||$\mu$|m (⁠|$\nu = 1.2$| THz). Schneider et al. (2022) then fitted as free parameters the dust temperature and column density. If we assume a dust/gas ratio of |$r = 1/100$|⁠, as in the HOTT column density map from Martin (private communication), the value of the dust mass opacity coefficient scaled to 1 THz, |$\kappa$|(1 THz) is 9.0 cm|$^2$| g|$^{-1}$| for Schneider et al. (2022), slightly less than 10 cm|$^2$| g|$^{-1}$| from Singh & Martin (2022), and the assumed frequency dependence is steeper than found by Singh & Martin (2022). Schneider et al. (2022) also find significantly lower dust temperatures compared with Singh & Martin (2022) with values distributed mainly in the range of 11–13 K. (Digital versions of the column density and dust temperature maps from Schneider et al. 2022 were kindly provided to us by N. Schneider and V. Ossenkopf-Okada).

The main difference between these two analyses is evidently the difference in adopted values for |$\beta$|⁠. Schneider et al. (2022) assumed a fixed value of |$\beta = 2$| while Singh & Martin (2022) use values which range between 1.4 and 1.7 with a mean of about 1.55, by interpolating from the sensitive but lower resolution maps of the Planck Collaboration. Since the assumed dust mass opacity coefficients are not very different (⁠|$\sim 10$| per cent) at the fiducial frequency of 1 THz, it is the difference in the adopted frequency dependencies that leads to different fits to the SED for dust temperature and optical depth. Lower opacity coefficients coupled with lower fitted dust temperatures (⁠|$\sim 12$| K versus |$\sim 19$| K), hence lower values of |$B_\nu (T_{{\rm dust}})$|⁠, require greater dust column densities to produce the observed far-IR intensities. As a consequence, Schneider et al. (2022) find higher total column densities than Singh & Martin (2022) for the same assumed gas/dust mass ratios. (See also appendix C in Shull & Panopoulou 2024 for further discussion of the uncertainties in dust models applied to far-IR images).

In conclusion, we have opted to use the map of total hydrogen column density derived from the analysis of Miville-Deschênes et al. (2017b) in our comparisons with the CO and HI distributions. The alternative approaches of Singh & Martin (2022) and of Schneider et al. (2022) have sufficiently large uncertainties in the critical dust-related parameters that we prefer the analysis rooted in the well-calibrated images from the Planck Collaboration as applied by Miville-Deschênes et al. (2017b) to their 250 |$\mu$|m image.