THE THREE-MM ULTIMATE MOPRA MILKY WAY SURVEY. I. SURVEY OVERVIEW, INITIAL DATA RELEASES, AND FIRST RESULTS

Even today, most radio telescopes have single-pixel detectors in their focal plane. For single-dish antennas operating over the ∼1–100 GHz range, angular resolutions are typically limited to about an arcminute. Without heterodyne focal-plane detector arrays (which are still rare and have most recently tended to focus on narrow-field astrochemical studies) to speed up the mapping, even this low resolution takes an extraordinary time to completely map large areas of sky at reasonable sensitivity. The conventional wisdom for some time has been that the dedicated, small-telescope, CfA/Nanten model was the only practical way to produce useful wide-field molecular maps.

Motivated by this conundrum, we devised two new techniques specifically optimized for the 22 m diameter Mopra⁹ antenna to dramatically improve the telescope's speed for wide-field mapping. We also implemented a third technique (pioneered in the CfA survey) that improves the image fidelity for all mapping projects. A detailed exposition of these techniques has been given by Barnes & Dame (2009), but for completeness we highlight the essential features here (summarized in Table 1).

Table 1.  Mopra Mapping Techniques

Mapping	Sampling	Finest	Atmospheric Stability	Relative	Relative
Mode	Interval	Resolution	Accommodated	Mapping Speeda	rms Noisea
Regular	2.000 s	33''	cannot adjust	1	1 (var.)
FM	0.256 s	33''	cannot adjust	8	$\sqrt{8}$ (var.)
BSM	any	66''	cannot adjust	4	1 (var.)
AM	any	any	>1 hr	≤1	1
ThrUMMS	0.256 s	66''	>1 hr	≤32	$\sqrt{8}$

Note.

^aThese factors are multiplicative: for ThrUMMS, all of FM+BSM+AM are used.

Download table as:  ASCIITypeset image

2.1.1. Fast Mapping with Mopra

2.1.2. Beam Sampled Mapping with Mopra

2.1.3. Active Mapping with Mopra

The third technique, new to Mopra but used at some other facilities, is a sky-related modification of BSM, but which can also be applied to regular "slow" OTF mapping. For most modern mm-wave radio telescope systems, the system temperature T_sys is dominated by the contribution from the atmosphere, T_atm. Apart from weather-related factors such as humidity, this of course depends on the airmass at the elevation e of the observation above the horizon, T_atm ∝ cosec(e). It turns out that we can easily (and exactly) accommodate the predictable, elevation-dependent component of T_sys by scaling ${\nu }_{{\rm{N}}}$ to T_sys(e)/T_sys(90°), since the slightly slower mapping at a larger ${\nu }_{{\rm{N}}}$ will increase the integration time in exact proportion to that required to adjust for the higher T_sys(e) compared to T_sys(90°). Thus,

$$\begin{eqnarray}&&{\nu }_{{\rm{N}}}(e)={\nu }_{{\rm{N}}}(90^\circ )\left(\displaystyle \frac{{T}_{\mathrm{sys}}(e)}{{T}_{\mathrm{sys}}(90^\circ )}\right),\end{eqnarray} \tag{ 1 }$$

with the minimum ${\nu }_{{\rm{N}}}$(90°) chosen to be 55 GHz. We call this technique "sky-active mapping," or AM.

In practice, the minimum T_sys(90°) for the scaling of ${\nu }_{{\rm{N}}}$ was chosen to be a conservative approximate global minimum in the best conditions (400 K at 115 GHz in our case), rather than estimated separately for each day at e = 90°. In that case, slowly varying (∼hourlong) and day-to-day weather changes could also be accommodated by this scheme just as easily. Only in the most variable weather is the AM method susceptible to bad spectral baselines, producing "weather stripes" in the maps, but this is true for any OTF mapping project, and the weather in such cases is, on average, not good enough to observe at all. This sky-active technique at least evens out much of the weather variation, and usually the scheme works quite well: an illustrative example of the combination of two of these techniques (AM+BSM) using Mopra's OTF system is shown in Figures 2 and 3, and is compared with a similar mosaic made with regular OTF mapping.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Mopra gives us another important advantage over mapping surveys made elsewhere: the MOPS spectrometer. This further combination makes ThrUMMS rather unique compared to previous molecular surveys of the Milky Way.

For "slow" OTF mapping, MOPS simultaneously records up to 16 high spectral resolution (∼0.1 km s⁻¹) IF "zoom modes" of width 137 MHz each (∼400 km s⁻¹ at 3 mm), which together can semi-arbitrarily cover an 8 GHz wide portion of the 3 mm band. Especially in the line-rich 90 GHz range, Mopra can therefore map 16 or more spectral lines extremely efficiently, as in CHaMP (Barnes et al. 2011) and other surveys of dense molecular cores and clumps. During FM, however, the data transfer rate from the spectrometer hardware to mass storage is maximized, and imposes a limit of four IF zooms on the data-taking. This limitation, however, does not affect the science goals of ThrUMMS, since there are only four transitions bright enough (at the sensitivity of FM) with a practical astrophysical application in the 107–115 GHz range. These transitions are listed in Table 2. We show in Figure 4 some sample spectra to illustrate the spectral data quality during VFM+AM.

Zoom In Zoom Out Reset image size

Table 2.  ThrUMMS 3 mm Band Spectral Lines

Species	Transition	Frequency	Efficienciesa		Resolutionb	rmsc	Science Topics
		(GHz)	ηb	ηc	( km s−1)	(K)
12CO	$J=1\to 0$	115.271	0.42	0.55	0.33	1.3	largest scales, GMC envelopes, temperatures, outflows
CN	$J=1\to 0$	113.5 (5 hf)d	0.42	0.56	0.33	0.9	dense clumps, Zeeman candidates
13CO	$J=1\to 0$	110.201	0.43	0.57	0.34	0.7	PDRs, fractionation, optical depth, velocity dispersion
C18O	$J=1\to 0$	109.782	0.43	0.57	0.34	0.7	column density, depletion, dense gas finder

Notes.

^aFollowing Kutner & Ulich (1981), the beam efficiency η_b is appropriate for obtaining ${T}_{\mathrm{mb}}$ of compact (∼beam-sized) sources, while the coupling efficiency η_c is appropriate for computing ${T}_{\mathrm{mb}}$ of more extended (several arcmin) sources. ^bThis is for a four-channel binning: the intrinsic velocity resolution is 0.085 km s⁻¹. ^cThis is for a four-channel binning: the intrinsic rms noise per full-resolution channel is twice as large. ^dThis line has nine hyperfine components, of which five are recorded in the IF zoom band.

Download table as:  ASCIITypeset image

In 2012 we completed (except as noted below) all mapping for the equatorial degree-wide strip (i.e., $| b| \lt 0\buildrel{\circ}\over{.} 5$) for two-thirds (360° > l > 300°) of the Fourth Quadrant (hereafter 4Q). In 2013 we began expanding the coverage to $| b| \lt 1^\circ$ over the same longitude range, and anticipate that this will be completed in 2015. While optimizing the data taking for efficiency and low noise, the filling of our 60° × 2° target area during the observing is somewhat random, but large contiguous portions of the 4Q are now available. A small number of fields are still affected by poor weather or equipment problems, such that they will require re-observation to maintain uniform sensitivity; other fields were not finished in one session. Despite this patchwork progress, there are no calibration "edges" between fields observed at different times or seasons, since the calibration is robust (see Sections 2.4 and 2.5). Moreover, because of our variable-${\nu }_{{\rm{N}}}$ approach to noise management, the rms varies only modestly between fields, with very little (if any) effect of such blended data being apparent in the merged data cubes.

All ThrUMMS maps are made with the VFM OTF raster direction in longitude (i.e., raster lines are at constant latitude). The absolute pointing accuracy is monitored and corrected to within 6'' using various 86 GHz SiO maser sources (R Car,IRSV 1540, AH Sco, VX Sgr; Indermuehle et al. 2013) between each map. Pointing variations rarely exceed 20'' at the end of each OTF map. (With the 24'' pixel size in each data cube, such pointing corrections, although slightly larger than those for smaller OTF maps at Mopra, are inconsequential here.) At the beginning of each daily observing session, a calibration scan is taken of the bright molecular clump BYF 99 in the η Carinae Nebula complex (Figure 5), and the T_sys measurements made to start the ${\nu }_{{\rm{N}}}$ scaling for subsequent VFM. OFF positions for each VFM field are taken from the list of Dame et al. (2001). Under good conditions each 05 field takes about 2^h of clock time, or somewhat longer as the T_sys rises. Fields are prioritized for mapping so that, whenever possible, they start at an hour angle (HA) near −1^h, and finish near HA = +1^h, minimizing the airmass and maximizing the mapping efficiency. This procedure is controlled by a robotic telescope observing script, implemented in the control software since 2010, and requiring only occasional oversight functions from our observing team.

Zoom In Zoom Out Reset image size

We use an augmented version of the Livedata/Grid-zilla package (Barnes et al. 2001), the usual standard for ATNF single-dish data reduction. This package provides a convenient environment for basic spectral calibration, baselining, separation of different IFs and spectral line data, gridding and convolution, and the formation of standard FITS (l, b, V) cubes. We describe our procedures next, including custom improvements to the package.

In Livedata, raw data for each observed field are bandpass-corrected and scaled by the measured ${T}_{\mathrm{sys}}$ to a relative antenna temperature scale ${T}_{{\rm{A}}}^{*}$, using the "reference" method (with OFF positions as described in Section 2.4). A linear baseline from a robust fit to signal-free channels (not including 100 channels at each end of each 4096-channel spectrum) is also subtracted, the frequency axis is rescaled for all topocentric to ${V}_{\mathrm{LSR}}$ corrections, and the result for each MOPS IF-zoom spectrum is written to a separate SDFITS file.

Most of our custom augmentations to this pipeline are effected between the Livedata and Gridzilla stages. First, absolute flux calibration is performed on the SDFITS files produced by Livedata, as illustrated in Figure 5. Next, automated removal of known observing or processing artifacts is performed. For example, in DR2 and earlier data releases, some artifacts (failed bandpass corrections for single 0.256 s samples, bad calibrations due to observing malfunctions, or bad-weather "stripes" where baseline subtraction becomes less reliable) were not corrected, and were manifest in the data products. While the fraction of data so affected is small, the effects on the cubes and maps can be quite large, since the individual errors can be of a large magnitude. Left uncorrected, such data produce noticeable deleterious effects in the maps, such as "bright" or "dark" spots, strong "edges," warped baselines, and so on. However, these errors were easy to identify and then either flag or correct. Thus in DR3, we have completely removed the single-pixel artifacts, systematic calibration errors, and some of the weather stripes, by new, custom, semi-automated routines written as a combination of IDL and Unix c-shell scripts. These additional calibration and editing steps improved the overall data quality to the point where the resulting noise levels were always consistent with theoretical expectations, given the observing conditions.

From this point the data processing resumes with the normal operation of Gridzilla, but with mapping parameters that are guided by the scientific motivation of defining the Galactic CO and CN emission to arcminute and sub-km s⁻¹ resolutions, given the unique features of the BSM technique. During the observations, one integration cycle occurs approximately every half-arcminute interval, in both the scanning and scan-step direction (i.e., l then b). This is comparable to the native resolution of the Mopra telescope at 115 GHz (approximately 33''; Ladd et al. 2005), so our mapping strategy undersamples that scale. To form images with an effective resolution that is well-sampled, we define gridding parameters in Gridzilla that are ∼twice as large as those for normal OTF maps. Thus, the map pixels are of size 24'' and the Gaussian convolving beam is 68'', truncated within a 36'' gridding kernel, giving an effective resolution of 72'' in the maps, with a typical contribution of 40 × 0.256 s samples per pixel to the final spectra. Therefore, the effective integration time is about 10 s pixel⁻¹, and the output from Gridzilla preserves the intrinsic velocity resolution of MOPS, as listed in the footnote to Table 2, with the typical noise level also given there. The final spatial sampling compares well with the typical spacing of ∼20''–25'' (≈1 pixel) between raster lines in our VFM OTF maps. Our mapping and data reduction procedures therefore return a product that is better than Nyquist-sampled at 72'' resolution, and cover an area of the Galaxy not feasible otherwise.

Table 2 also summarizes some parameters of the data cubes for the four spectral lines we are mapping. The per channel rms noise values drop with frequency since these lines are at the edge of the atmospheric 3 mm band; thus, the ${T}_{\mathrm{sys}}$ for ¹³CO and C¹⁸O are almost a factor of 2 less than for the ¹²CO data. Each species' cube is produced on the ${T}_{{\rm{A}}}^{*}$ scale, and must be converted to ${T}_{{\rm{r}}}^{*}$ = T_mb for physical analyses. Table 2 gives both the main beam (η_b) and inner error beam (η_c) efficiencies we recommend for each spectral line, based on the measurements of Ladd et al. (2005). These efficiencies are respectively appropriate (Kutner & Ulich 1981) for compact sources of ∼one intrinsic Mopra beam or less, and for slightly more extended (a ∼ few arcmin) sources. For most emission seen in the ThrUMMS maps, η_c will be more often used, and the DR cubes are all scaled to ${T}_{{\rm{r}}}^{*}$ by η_c (although see the discussion in Section 2.6 about refinements to this scaling).

The next step in our augmented processing pipeline after the Gridzilla stage is to calculate moment maps for the line emission. We describe these methods in Section 3.3.

The final step in the data reduction pipeline is to check that we have constructed maps that are at least consistent with previous surveys, e.g., the CfA data of Dame et al. (2001). As can be seen in Figure 6, the ThrUMMS data do indeed compare well: averaging over a variety of regions, we derive an overall scale factor between the two surveys of ${T}_{{\rm{r}}}^{*}$(ThrUMMS) = (0.90 ± 0.05) $\times \;{T}_{{\rm{r}}}^{*}$(CfA), which is within the overall calibration uncertainties given in Figure 5. We consider this result (close to unity, with a small variance) a strong verification of the quality-control procedures described above.

Zoom In Zoom Out Reset image size

Even so, the conversion factor is slightly but systematically below 1. This may be due to two effects. One is that the effective coupling between the sky emission and the inner error beam (∼5') of the Mopra telescope is not as uniform as is implicitly assumed by the use of η_c = 0.55–0.57 in Table 2. The conversion suggests instead that an efficiency of 0.50 may be a better value, i.e., about 10% below η_c, but 15% above η_b. This would be true if the sky emission only partly fills (perhaps to ∼60% of) the inner error beam, giving an effective η between these nominal values. However, the appearance of the maps (especially the ¹²CO from which this result is obtained) is not completely consistent with this idea: over large areas the ¹²CO cloud structure is very extended, although it is compact in other areas, so we might expect a more variable conversion between the ThrUMMS and CfA brightness scales, rather than the fairly systematic trend evident in Figure 6.

An additional factor may be that our BSM technique is systematically missing a few percent of the line flux from all mapped areas, due to the undersampling of the emission structure (see Section 2.1.2). That is, between raster lines there may be small, bright, compact (i.e., beam-sized) structures whose total flux would not be fully measured by the beam-sampling, yielding an apparently lower brightness scale for the maps after interpolating between the raster lines. Another way to understand this issue is by analogy with interferometry: there are scales to which our maps are less sensitive than if they were properly Nyquist-sampled. Tests of this effect in ¹²CO maps of the CHaMP areas (which are Nyquist-sampled; P. Barnes et al., in preparation) suggest that it may be relevant, but we defer that discussion to that work. For the purposes of this paper, we use these explanations as working assumptions to account for the 0.90 conversion factor. Thus, while we quote all results herein on the ${T}_{{\rm{r}}}^{*}$ scale as derived with the tabulated η_c factors, the reader should be aware that these brightness levels may need to be increased by a further ∼10%.

From its inception, ThrUMMS was designed to be an open-access project, with frequent Data Releases as progress is made in the observing and data reduction. As of 2014 May, we have placed a large amount of data on the website http://www.astro.ufl.edu/thrumms, freely available to the community. For the ¹²CO and ¹³CO cubes, the data are at DR3 level and organized into 6°-wide "Sectors" (i.e., each including up to 48 contiguous half-degree fields; see Figures 7–11). For the C¹⁸O and CN data, we currently have available the DR2 cubes, organized into 2°-wide blocks; DR4 is planned to include the augmented artifact-correction for the latter two species, as well as wider area coverage from the 2014 observing season, and is anticipated to be available in mid-to-late 2015. At the present time, these data cubes cover ∼80% of the 60° × 2° target area for all four species. With the anticipated completion of observing in 2015, later Data Releases will continue to be provided as processing continues and further data products can be assembled.

Zoom In Zoom Out Reset image size

We use the Miriad package (Sault et al. 1995) to process the data cubes into moment maps and LV diagrams, and perform other basic analyses, with visualization provided by the Karma package (Gooch 1997). Figure 7 gives some sample comparisons between ThrUMMS maps and other surveys, while Figure 8 shows sample VFM data for our four molecular transitions. Note the high S/N of the ¹²CO and ¹³CO maps, and even though the C¹⁸O and CN maps have noticeably smaller S/N, the peaks in these maps (around the positions of the brightest ¹²CO or ¹³CO emission) are still detected with good reliability (S/N ≳ 5–10). In Figure 9 we show further examples of higher moment and LV maps in the ¹²CO line, demonstrating the very high image fidelity achieved across extended emission regions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The example moments maps of Figure 9 were processed with a smooth-and-mask (SAM) algorithm in order to improve the fidelity of these products. As part of this algorithm, the DR3 cubes are binned by four channels in velocity (see Table 2), which improves the S/N in the public data cubes and reduces the file sizes. (The full-resolution cubes are, however, still available on request.) Our SAM algorithm is similar to those described by Rots et al. (1990), Dame et al. (2001), or Dame (2011). We form a data cube that has been smoothed/convolved in all three dimensions by a factor of 2 (i.e., beyond the four-channel binning of the public data), a 5σ_rms threshold was calculated, and voxels (where a voxel is defined here as a single (l, b, V) element) below this level in the smoothed cube were used to define a blanking mask for voxels in the unsmoothed data cube. Moments are then calculated only with the unmasked voxels.

Construction of moment maps, however, depends strongly on the particular investigation one wishes to pursue with the data cubes. There is a large range of input parameter combinations possible for such maps, including angular extent, velocity range, projections, SAM parameters, uncertainty estimation, etc. For this reason, moment maps using the above techniques and based on the public data cubes are not currently part of the DRs, but are available upon request (such requests are processed on a best-efforts basis).

We do plan, however, to include sets of standardized moment maps in future DRs, once we have re-observed most of the weather-damaged areas and made the data cubes' noise statistics more uniform. Meanwhile, users can use many standard astronomical or image processing packages to perform their own analyses, including complex sorting and assessment algorithms such as Principal Component Analysis (Heyer & Schloerb 1997; Lo et al. 2009), ISM structural analyses like Velocity Component Analysis (Lazarian & Pogosyan 2000), or Spectral Correlation Functions (Rosolowsky et al. 1999).

The full survey data cubes provide an extremely rich data set for current and future investigations. This richness means it is difficult to render, whether on the printed page or a web browser, the detailed physical insights contained in the data. Nevertheless, the reader may gain some intuitive feeling for this richness by means of color-compositing the different species' integrated intensity maps. We show examples of this in Figure 10, combining the DR3 ¹²CO, ¹³CO, and C¹⁸O integrated intensities, and in Figure 11 for pre-DR4 versions of the ¹²CO, ¹³CO, and C¹⁸O data, each showing clear positional variations in the line ratios as variations in the rendered colors.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Similarly, compressed image versions of the spectral line integrated intensities across the full survey area, in all four species, are given in the Appendix. The reader should note, however, that the Appendix images, and the data used in the analysis presented below, have been manually edited from the DR3 versions available on the ThrUMMS web site to exclude most of the artifacts generated during poorer observing conditions ("weather stripes"). These excisions are manifest in the Appendix figures as half-degree-wide gaps in certain fields; the gaps themselves are not the weather stripes. These blanked areas cover only ∼2% of the mapped area, however.

With our comprehensive coverage of a large portion of the Milky Way in multiple species, we have an unprecedented opportunity to use the molecular lines to investigate global GMC conditions. We examine now the variation in the CO line ratios (and consequently, the variation in the physical conditions in GMCs) with position in the Galaxy, as shown visually in Figures 10 and 11. This rendering can be intuitively understood in terms of radiative transfer physics in the following way.

The clouds rendered in redder colors show where the ¹²CO is relatively much brighter on average than either ¹³CO or C¹⁸O, i.e., where the ¹²CO/¹³CO and ¹²CO/C¹⁸O brightness ratios are near their maximal values. This must therefore indicate areas of relatively low optical depth, since the intrinsic abundance ratios are also quite large (any line saturation = high optical depth would produce brightness ratios significantly less than the respective abundance ratios). At the same time, the ¹²CO lines in such areas are quite bright, typically with ${T}_{{\rm{r}}}^{*}$ ≳ 30 K. So the "red" clouds are typically those of low τ and high T_ex, or "warm and translucent" clouds.

In contrast, green or even cyan/blue areas are where either ¹³CO, C¹⁸O, or both are much brighter compared to ¹²CO than on average. This is only likely to occur where the optical depth is high enough to produce some level of saturation in either line, so that the ¹²CO/¹³CO and ¹²CO/C¹⁸O brightness ratios are closer to unity than average. At the same time, the ¹²CO emission in such areas is not particularly bright; coupled with the higher optical depth, this means that the excitation temperature in such clouds is likely to be somewhat low. So the "green" and "blue" clouds will typically represent areas of high τ and low T_ex, or "cold and opaque" clouds.

Therefore, we understand that the line ratios provide a simple yet powerful way to derive basic physical conditions in Galactic GMCs from an elementary radiative transfer analysis. The rendered colors make these relationships intuitively clear: they are set so "white" or "neutral" colors represent the average line ratios and not where the line ratios are close to 1.

Other workers have used a similar approach to derive cloud properties from CO line ratios; examples include Myers et al. (1983), Wong et al. (2008), and Hindson et al. (2013). Each of these, however, uses some simplifying assumptions about the excitation temperature, abundance ratios, or other radiative transfer parameters, and typically for relatively small samples of clouds covering at most tens of parsecs. Ours is the first such study of which we are aware that attempts this analysis uniformly across many kiloparsecs of the Galaxy; at the same time, we retain the pc-scale resolution necessary to characterize the clumps that give rise to individual star clusters.

In the analysis that follows, we use the full 3D iso-CO data cubes, but binned in velocity to 1 km s⁻¹ channels, to compute all physical quantities. We do this in order to fully track the variations of ${T}_{\mathrm{ex}}$. τ, and other parameters with position and velocity, but also to boost the S/N in the C¹⁸O data especially. In what follows, therefore, the notation $\int {dV}$ refers only to these 1 km s⁻¹ integrals, which are still voxels in the binned cubes, unless otherwise stated.

For abundance ratios R₁₃ = [¹²CO]/[¹³CO] and R₁₈ = [¹³CO]/[C¹⁸O], the optical depths in the different isotopologues will follow τ₁₂ = R₁₃τ₁₃ and τ₁₃ = R₁₈τ₁₈. Then the radiative transfer equation for each spectral line (in the Rayleigh–Jeans limit, where J_ν ∝ T),

$$\begin{eqnarray}&&{T}_{\mathrm{mb}}=({T}_{\mathrm{ex}}-{T}_{\mathrm{bg}})(1-{e}^{-\tau }),\end{eqnarray} \tag{ 2 }$$

can be combined to convert the line ratios to an optical depth, with the further assumption of a common excitation temperature T_ex. In general, this last condition may not hold everywhere: the lower ¹³CO optical depth means that this species will more often be subthermally excited compared to ¹²CO, which has a lower "effective" density due to strong radiative trapping (Evans 1999; Heyer & Dame 2015). Nevertheless, this assumption will only act to underestimate the resultant ¹²CO column density (see Section 4.2).

With these assumptions, the brightness ratio in the ThrUMMS maps should obey

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{13}}{{T}_{12}}=\displaystyle \frac{1-{e}^{-{\tau }_{13}}}{1-{e}^{-{R}_{13}{\tau }_{13}}},\end{eqnarray} \tag{ 3 }$$

and a similar equation applies for C¹⁸O and ¹³CO. For typical values of R and τ (especially for R₁₃), the denominator in Equation (3) will often be close to unity. Also, for any reasonable value of R₁₃, the ¹²CO will be virtually opaque almost everywhere, and T₁₂ = $({T}_{\mathrm{ex}}-{T}_{\mathrm{bg}})$.

The visual representations of line ratio variations, as seen in Figures 10 and 11, can be quantified and understood in Figure 12. Here we see that, in the left panel, there is a very wide range of ratios log(${I}_{{}^{13}\mathrm{CO}}$/${I}_{{}^{12}\mathrm{CO}}$) = log($\int {T}_{{}^{13}\mathrm{CO}}$dV/$\int {T}_{{}^{12}\mathrm{CO}}$dV), from ∼+0.2 in the upper left of this panel (seen as green areas in Figures 10 and 11) to ∼−1.1 in the lower right of Figure 12 (redder areas in Figures 10 and 11). Similarly in the right panel, we see a range of ratios for log(${I}_{{{\rm{C}}}^{18}{\rm{O}}}$/${I}_{{}^{13}\mathrm{CO}}$) = log($\int {T}_{{{\rm{C}}}^{18}{\rm{O}}}$dV/$\int {T}_{{}^{13}\mathrm{CO}}$dV), from ∼+0.5 in the upper left of this panel (seen as blue or cyan areas in Figure 10) to ∼−1.2 in the lower right of Figure 12 (green or yellow areas in Figure 10).

Zoom In Zoom Out Reset image size

Now we can apply Equation (3) to the data in Figure 12, where in each panel we show lines of constant τ for an assumed R. As expected, low optical depth for the line pairs occurs to the right of each panel, while high optical depth lies to the left. Interestingly, in Figure 12 we see several groups of high-brightness emission peaks (i.e., the brightest clumps in Figures 10, 11) as separate extensions of points, parallel to lines of constant τ, in the distribution to the upper right of each panel. These extensions apparently correspond to groups of bright clumps (complexes) with different regional τ and T_ex (see Section 4.3).

The plots in Figure 12 can be combined into a single "CO-color" plot (Figure 13), where a voxel's position in this plot corresponds to its color rendering in images like Figures 10 or 11. This can be compared with the analytical relations derived above (Equations (2) and (3)) which are shown as overlaid lines in Figure 13. From this comparison it is immediately obvious that we cannot model the most extreme "green" and "blue" voxels with the assumptions used above (i.e., a common T_ex and fixed ratios R > 1). For this small number of voxels, the model would require either (1) ratios R < 1 (e.g., due to strong fractionation or other selective effects), or (2) a lower T_ex in the more abundant and optically thick species (e.g., due to a colder foreground layer and warmer cloud center). Fractionation effects for both C¹⁸O and especially ¹³CO (Glassgold et al. 1985) do occur, as does self-absorpton from colder cloud envelopes, so undoubtedly these both contribute to deviations from the simple radiative transfer picture presented here. Note also that single clouds may be represented by a range of loci in Figure 13. So, while we can understand most of the line ratio variations in Figures 10 and 11, i.e., the lower and leftmost boundaries of the distribution in Figure 13, eventually a more sophisticated radiative transfer approach will be needed to fully understand all the line ratios, because the interplay between optical depth, excitation, and abundance is more complex than the above treatment allows.

Zoom In Zoom Out Reset image size

Nevertheless, we can use our treatment to estimate the molecular column density distribution in the Galactic Plane (Section 4.2). We use the line ratios to evaluate τ via Equation (3) at each 1 km s⁻¹-binned voxel, and then invert Equation (2) to also obtain T_ex at each voxel, producing all-4Q cubes of both parameters. Finally, we form moment maps from each of these cubes in the normal way, e.g., integrating over large velocity ranges to obtain the total τ, peak or average T_ex, or total N_CO (see also Section 4.2). We show examples of these results in Figure 14; although these are moment maps, the structures they exhibit are typical of the cubes over which they were integrated, although there are fewer visible complexities in any given velocity channel.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The τ panel in Figure 14 (top) shows a very high-contrast image: all the highest optical depth structures are quite compact, typically having size scales of ∼1'–5'. It may not be coincidental that this is also the size scale (∼a few pc at 3 kpc distances) of dense clumps, IRDCs, and clusters. Lower optical depth structures are not quite so spatially compact, and are somewhat widely distributed.

In contrast, the ${T}_{\mathrm{ex}}$ panel of Figure 14 (middle) shows a dramatic difference: nearly all structures are relatively large-scale, and somewhat amorphous and/or low-contrast. While the warmest structures are quite extended, at least 5', moderate-temperature clouds extend over even larger degree-wide (∼50 pc at 3 kpc) scales. Again, it is tempting to draw a connection between these scales and the similar sizes of GMCs and classic H ii regions, 10–100 pc.

We note an interesting consequence of this radiative transfer approach. In Figure 13, τ₁₈ maximizes near a value of 0.6; for any higher optical depth, all CO isotopologues would begin to saturate at T_mb = ${T}_{\mathrm{ex}}-{T}_{\mathrm{bg}}$. Such clouds would have τ indistinguishable from infinity, whatever the T_ex (all of the R₁₈ curves converge to where both line ratios are 1). Even so, τ₁₈ = 0.6 corresponds approximately to τ₁₂ ∼ 250 for any reasonable values of the abundance ratios R. This means that there are a large number of clouds in the ThrUMMS maps with very large column densities (e.g., near the top-right of the plots in Figure 13, with blue, cyan, or green ratios), some of which could masquerade as lower column densiy clouds if one were to convert from I₁₂ with a standard X-factor.

A subject of much study in the literature is the conversion of observed molecular emission into estimates of molecular column density or mass, which are important for a number of questions in astrophysics, e.g., in star formation when measuring cloud stability against gravitational collapse, or across disk galaxies for studies of the Kennicutt–Schmidt relation, or for estimating the contribution of the molecular ISM to a galaxy's dynamics, etc. Dame et al. (2001) summarized this topic and, when converting CO emission into mass, obtained a best average estimate for the X-factor of

$$\begin{eqnarray}{X}_{\mathrm{CO}} & = & {N}_{{H}_{2}}/{I}_{\mathrm{CO}}\\ & = & 1.8\times {10}^{24}\;{{\rm{H}}}_{2}\;\mathrm{molecules}\;{{\rm{m}}}^{-2}/({\rm{K}}\;\mathrm{km}\;{{\rm{s}}}^{-1})\end{eqnarray} \tag{ 4 }$$

in the solar neighborhood, based on high-latitude data in the CfA survey. But as Dame et al. (2001) point out, this value is expected to vary locally in the Galactic midplane (i.e., at low latitudes) due to various factors.

Having derived τ and T_ex at each voxel, we can use the ThrUMMS data to make a unique contribution to this exercise with a CO column density map, via

$$\begin{eqnarray}&&{N}_{\mathrm{CO}}=\displaystyle \frac{3\;h}{8{\pi }^{3}{\mu }^{2}}\;\displaystyle \frac{Q({T}_{\mathrm{ex}}){e}^{{E}_{{\rm{u}}}/{{kT}}_{\mathrm{ex}}}}{{J}_{{\rm{u}}}(1-{e}^{-h\nu /{{kT}}_{\mathrm{ex}}})}\;\int \tau {dV},\end{eqnarray} \tag{ 5 }$$

where μ is the CO dipole moment, Q is the rotational partition function, E_u and J_u are the energy and quantum number of the upper level of the transition at frequency ν, and the integral is over the velocity range of the emission line. Here, instead of using a constant X_CO factor to convert the observed I_CO into a molecular gas mass, as is sometimes done in the literature despite Dame et al.'s (2001) caveat against doing so, we can use the (τ, T_ex) estimates above (Section 4.1) to compute a more realistic CO column density. This can then be compared with other estimates of molecular gas mass (e.g., SED fits to the Hi-GAL data) to derive CO abundance maps.

The computation of Equation (5) across a spatially resolved map with ∼8 × 10⁸ independent (l, b, V) resolution elements, complete with a calculation of the partition function at each voxel and appropriately propagating uncertainties in each quantity, is a non-trivial exercise. Moreover, in many locations the ¹³CO and especially C¹⁸O can be quite weak, so there are S/N limitations in applying this method to our full-resolution data cubes. Nevertheless, as described in Section 4.1 we have performed this computation on cubes that have been binned in velocity to 1 km s⁻¹ channels. After so computing the τ and T_ex in each such voxel, we use Equation (5) to derive cubes of column density per channel.

We also describe here a simpler calculation, which gives a very good approximation to Equation (5) in cases where insufficient computational power is available. The temperature-dependent terms in Equation (5) can be collected and written as a closed function of T_ex, which we call the modified temperature function ${{\mathcal{T}}}_{\mathrm{mod}}$:

$$\begin{eqnarray}&&{{\mathcal{T}}}_{\mathrm{mod}}=\displaystyle \frac{Q({T}_{\mathrm{ex}}){e}^{{E}_{{\rm{u}}}/{{kT}}_{\mathrm{ex}}}}{{J}_{{\rm{u}}}(1-{e}^{-h\nu /{{kT}}_{\mathrm{ex}}})}.\end{eqnarray} \tag{ 6 }$$

Over a typical range of T_ex = 5–100 K, the following formula for ${{\mathcal{T}}}_{\mathrm{mod}}$ gives deviations of <1% from the exact expression (Equation (6)), without having to calculate the partition function Q explicitly:

$$\begin{eqnarray}&&{\mathrm{log}}^{2}({{\mathcal{T}}}_{\mathrm{mod}})=1+{\mathrm{log}}^{2}[0.053\;{T}_{\mathrm{ex}}^{1.97}].\end{eqnarray} \tag{ 7 }$$

Having computed the ${{\mathcal{T}}}_{\mathrm{mod}}$ cube from the ${T}_{\mathrm{ex}}$ cube using either Equations (6) or (7), we then obtain the column density cube from this and the τ cube,

$$\begin{eqnarray}&&{N}_{\mathrm{CO}}=5.57\times {10}^{18}\;{{\rm{m}}}^{-2}\;{{\mathcal{T}}}_{\mathrm{mod}}\int \tau {dV},\end{eqnarray} \tag{ 8 }$$

where the integral over V means per 1 km s⁻¹ channel. We emphasize, however, that we have actually done the full calculation from Equation (5)/(6) in this study.

A moment map of the column density cube from Equation (5) is shown in Figure 14 (bottom panel). This now embodies a much larger "conversion factor" than in Dame et al. (2001). In the cool and opaque case, we might have typical values τ ≳ 100 and T_ex ∼ 10 K. For a cloud with an effective line width ∼5 km s⁻¹, N_CO ∼ 4.1 × 10²² m⁻². In fact, the bottom panel of Figure 14 shows that column densities computed in this way can reach a peak value some 5× higher than this, or N_CO ∼ 2 × 10²³ m⁻². Assuming that the ¹²CO abundance is 10⁻⁴ (although see Section 4.3), this is equivalent to ${N}_{{{\rm{H}}}_{2}}$ ∼ 0.4–2 × 10²⁷ m⁻², or a mass column Σ ∼ 760–3700 M_⊙ pc⁻². Even the lower end of this range is close to the mass columns in parsec-scale "dense clumps" seen in previous surveys (e.g., Barnes et al. 2011), but is much higher than typical GMC mass columns. We can convert the above cool, opaque example to an equivalent X_CO ∼ (8–40) × 10²⁴ m⁻²/(K km s⁻¹), and so realize that the Milky Way's molecular gas mass may need to be revised upwards by a significant amount, at least in some locations; see the discussion in Section 4.3.

Equally remarkable as the actual numbers are the structural revelations in Figure 14. As has been seen in many previous studies, the integrated intensity iso-CO maps (Appendix, Figure 16, panles (a) to (f)) have a rather complex combination of structures. The respective data cubes reveal that this complexity breaks down into individual regions or clouds at various ${V}_{\mathrm{LSR}}$ , which may be clumpy and compact, or filamentary, or amorphous on larger angular scales. However, the cubes by themselves do not reveal any obvious systematic trend among these structures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In contrast, the τ, ${T}_{\mathrm{ex}}$, and N_CO cubes show systematically different structural characteristics. As already noted in Section 4.1, the τ cubes are clumpy and high-contrast, while the ${T}_{\mathrm{ex}}$ cubes are amorphous and low-contrast. In Figure 14, the N_CO panel (bottom) reveals the combination of the τ and ${T}_{\mathrm{ex}}$: a high-contrast and more filamentary medium than either of the above constituents. One is naturally drawn to compare such images with the almost ubiquitous filamentary structures evident in the Hi-GAL maps.

Based on the distribution of points in Figure 13 and the conversion embodied in Equation (5), we can in principle compute mass, density, column density, or any other distribution functions for the entire molecular cloud population of the 4Q, such as in Figure 15 (discussed next). For many derived parameters (such as mass), we also need distances to individual clouds to obtain them from the column density computed here. However, assignment of distances to all the structures we see in the cubes is far beyond the scope of this paper.

Despite this, our results have a number of significant implications for studies of molecular cloud populations in disk galaxies. For example, a variety of observational measurements and chemical models indicate that the fractional CO abundance relative to molecular hydrogen typically¹⁰ lies in the range 10^−3.5–10^−4.5; a mean value of 10⁻⁴ is often used in the literature. As discussed in the last section, this means that we should be able to simply convert our computed N_CO to ${N}_{{{\rm{H}}}_{2}}$ by such a factor. But we also have the standard conversion of I_CO to ${N}_{{{\rm{H}}}_{2}}$ using Equation (4). For these two conversions to be consistent, the various quantities must be related. To make these relationships clearer, we rewrite Equation (4) as

$$\begin{eqnarray}&&{N}_{{{\rm{H}}}_{2}}={X}_{\mathrm{CO}}{I}_{\mathrm{CO}}={R}_{12}{N}_{\mathrm{CO}},\end{eqnarray} \tag{ 9 }$$

where we define the ratio R₁₂ from the actual CO abundance in the gas, R₁₂ ≡ ${N}_{{{\rm{H}}}_{2}}$/N_CO, or typically 10^+3.5 to 10^+4.5 from above. Note that this definition makes the gas abundance ratio R₁₂ different to the empirically determined X_CO conversion factor. Then

$$\begin{eqnarray}&&{N}_{\mathrm{CO}}/{I}_{\mathrm{CO}}={X}_{\mathrm{CO}}/{R}_{12}.\end{eqnarray} \tag{ 10 }$$

With the standard values on the RHS of Equation (10), the ratio on the LHS should be a constant,

$$\begin{eqnarray}&&{N}_{\mathrm{CO}}/{I}_{\mathrm{CO}}=1.8\times {10}^{20}\;\mathrm{CO}\;\mathrm{molecules}\;{{\rm{m}}}^{-2}/({\rm{K}}\;\mathrm{km}\;{{\rm{s}}}^{-1}).\end{eqnarray} \tag{ 11 }$$

However, when we plot Equation (11) versus I_CO in Figure 15, we see a marked variation in the actual ratio X/R as shown in Equation (10), instead of the flat line that would be expected from Equation (11) with fixed values of X and R (shown as a dashed blue line in Figure 15).

Even at the lowest values of I_CO (≲6 K km s⁻¹ per 1 km s⁻¹-wide channel) in Figure 15, the modal fit (blue line) of X/R already begins above the standard value (dashed line), while the mean values (green line) lie around (3–4) × 10²⁰ CO molecules m⁻²/(K km s⁻¹), or twice the standard value. If one assumes X_CO is fixed to the standard value, this implies that R₁₂ = 5 × 10³, or half its standard value. On the other hand, if one assumes R₁₂ is fixed to its standard value, this implies that X_CO = (3–4) × 10²⁴ H₂ molecules m⁻²/(K km s⁻¹), or twice its standard value. As I_CO rises, the modal fit in X/R also rises further to almost 10²¹ CO molecules m⁻²/(K km s⁻¹) at I_CO = 40 K km s⁻¹ per 1 km s⁻¹-wide channel. This means that at these brightness levels, either X_CO is 4× higher than standard, or R_CO is 4× lower than standard, or some combination of these two options applies. At I_CO > 20 km s⁻¹, the X_CO/R₁₂ distribution actually bifurcates¹¹ into two regimes, one with the ratio up to 2× higher than the standard value, and another up to 10× higher. A weighted power-law fit to the upper trend (sloping blue line in the plot) gives

$$\begin{eqnarray}&&X/R\propto {{I}_{\mathrm{CO}}}^{0.38}\end{eqnarray} \tag{ 12 }$$

across the measured range of I_CO, or equivalently,

$$\begin{eqnarray}&&{N}_{\mathrm{CO}}=1.6\times {10}^{20}\;{{\rm{m}}}^{-2}\;{I}_{\mathrm{CO}}^{1.38}.\end{eqnarray} \tag{ 13 }$$

This is a remarkable result for several reasons, but mainly because in retrospect, it should have been more widely recognized before now. For example, the fact that Equation (13) is superlinear falls naturally out of our very simple radiative transfer assumptions, and theoretically was predicted on this basis (Krumholz & Thompson 2007; Narayanan et al. 2008). Also, because Equation (13) is superlinear, it is no longer as physically meaningful to talk about a single X-factor: the direct conversion from CO integrated intensity to hydrogen column density can be seen as an illusion that is correct only in the mean, and masks the very real effects of high optical depth in the ¹²CO line. Thus, the scale factor in Equation (13) has units of molecules m⁻²/(K km s⁻¹)^1.38, and so dimensionally is not a "conversion factor."

Nevertheless, we can compute average quantities from our data, which may still be useful in spatially under-resolved studies, such as extragalactic or high-redshift star formation settings. The equivalent median value of the X-factor for all our data is X_CO = 3.8 × 10²⁴ m⁻²/(K km s⁻¹) or 2.1× the standard value, while a column-weighted median conversion is even larger than this, X_CO = 5.1 × 10²⁴ m⁻²/(K km s⁻¹) or 2.8× the standard value.

However, Equations (12)–(13) represent just the modal trend in Figure 15. It should be clear that the range of the variations in X_CO/R₁₂ is far too large to be attributed to changes in just one of these parameters. Even the 1σ variation from the modal trend in Figure 15 is 0.3 in the log, or a factor of 2. But the extrema in Figure 15 show X/R ratios from about 30× larger to 100× smaller than standard. This represents >3 orders of magnitude empirical uncertainty in the conversion of I_CO to ${N}_{{{\rm{H}}}_{2}}$. This variation exceeds any uncertainty in the assumptions we have used in this analysis, such as the assumption of LTE at a single ${T}_{\mathrm{ex}}$, or values of the ratios R, or even of the constancy of X_CO in the face of the metallicity gradient in the Galactic disk (Heyer & Dame 2015).

Such uncertainties aside, in the mean it is likely that, as I_CO rises from 4 to 40 K km s⁻¹ in each 1 km s⁻¹-wide channel, X_CO is very likely rising and that R₁₂ is probably falling (meaning that the gas phase abundance of ¹²CO relative to H₂ is probably rising). Therefore, the conversion of I_CO to an equivalent mass surface density is not as simple as has often been assumed. While we cannot yet derive the exact corrections to cloud masses, and the Milky Way's total molecular mass, without clouds' individual distances, this result has obvious implications for studies of the molecular mass distribution in the Milky Way, star formation laws in disk galaxies such as the Kennicutt–Schmidt relations, and other topics (Leroy et al. 2008; Lada 2015). At the very least, it implies that such relations must be recalibrated for the highest optical depth clouds and regions, which in our data are widely distributed across the Galactic Plane, and seem to have not yet been adequately taken into account.