THE PHYSICAL PROPERTIES OF HIGH-MASS STAR-FORMING CLUMPS: A SYSTEMATIC COMPARISON OF MOLECULAR TRACERS

Twenty-seven sources were mapped with BEARS on April 22 through 26, 2004 using the Nobeyama Radio Observatory 45 m telescope. BEARS is a 25 element focal plane array of SIS mixers (Sunada et al. 2000). The beams are spaced at 411 on the sky with a beamsize of approximately 178 (see Tatematsu et al. 2004). To map a source at nearly full beam resolution, the array must be dithered in 4  ×  4 grid with a spacing of 2055. Also, to minimize systematic calibration difference between the individual beams, the source was mapped with the array rotated by +90° and −90°. Each source was observed in position-switched mode with OFF positions that were determined to be clean of emission in CS J = 5 → 4 (Shirley et al. 2003). The final map covers approximately 3' × 3' with 100 spectra.

The BEARS backend is a digital autocorrelator with 1024 channels. The spectral resolution is 37.8 kHz resulting in 0.12 km s⁻¹ resolution at 93.7 GHz. Each spectrum was centered at v_LSR for the strongest hyperfine component of N₂H⁺ J = 1 → 0 (FF₁ = 23 → 12 at 93.1737767 GHz). During data reduction, the spectra were smoothed by four channels to improve the signal-to-noise; the final spectral resolution is 0.49 km s⁻¹.

The spectra were placed onto the T*_A scale by comparing observations of each beam in the array with the Nobeyama S100 receiver. T_mb is calculated by using a main-beam efficiency of the telescope with S100 of 0.515. The S100 receiver uses the standard chopper-wheel calibration method (Penzias & Burrus 1973). The bright source, G40.50, was observed with each beam in the array on two separate days. The peak of the N₂H⁺ J = 1 → 0 emission was used to compute the calibration corrections for each beam. We estimate that the calibration is consistent with T*_A determined by S100 to within 15%. Weather conditions were good to excellent during the observations. The typical wind speed was below 10 m s⁻¹.

The spectra were reduced using the NEWSTAR software package of the Nobeyama Radio Observatory. A linear baseline was removed from each spectrum and the spectra at the same offset position were averaged. Since the actual beam centers on the sky do not form a perfectly regular grid, the spectra were interpolated onto a regular grid using a Gaussian smoothing function. The final spectra were written out as an SPFITS file (Greisen & Harten 1981) and were read into AIPS++ where the individual spectral arrays could be manipulated. Gaussian fits were performed using the profile fitter in AIPS++.

Observations of N₂H⁺ J = 3 → 2, HCO⁺ J = 3 → 2, and HCS⁺ J = 6 → 5 of the 27 sources were performed at the Caltech Submillimeter Observatory 10.4 m telescope on 10 nights between 2002 June and 2004 July. The 230 GHz single-pixel receiver (Kooi et al. 1992, 1998) with a 50 MHz AOS backend was used for all observations. The average spectral resolution was 0.15 km s⁻¹; however, the spectra were smoothed to a velocity resolution of 0.21 km s⁻¹ (N₂H⁺, HCO⁺) and 0.14 km s⁻¹ (HCS⁺). We made oversampled, boustrophedonic (alternating scan direction between left-to-right and right-to-left lines) On-The-Fly maps of N₂H⁺ and HCO⁺ on a grid with 10'' resolution (Mangum et al. 2007) while HCS⁺ observations were single position-switched spectra toward the 350 μm dust continuum peak position.

The standard chopper calibration-wheel method was used to calibrate the spectra on the T*_A scale (Penzias & Burrus 1973). Calibration on planets (Mars, Jupiter, and Venus) was used to convert to the T_mb scale assuming that the sidelobes do not contribute significantly (Kutner & Ulich 1981). Pointing was performed every hour on planets or bright high-mass star-forming cores (W28A2(1)) when planets were not available. We estimate our pointing is better than 10'' when no planets were available and is better than 5'' with planets.

All 27 sources were mapped using the 1 mm ALMA prototype sideband-separating receiver on the 10 m Heinrich Hertz Telescope (HHT) on Mt. Graham, Arizona over six nights between 2008 October and 2009 June. We observed in dual polarization four IF modes (5.0 GHz IF), with HNC J = 3 → 2 (271.9811 GHz) centered in the upper sideband and SO N_J = 6₇ → 5₆, C₂H N_J = 3_J → 2_J, and CH₃OH-E J_K = 2₁ → 1₀ in the lower sideband. We also observed the isotopologues HN¹³C J = 3 → 2 (261.263 GHz) centered in the upper sideband and ³⁴SO N = 6₇ → 5₆ (256.877 GHz) also centered in the upper sideband between 2009 November and 2010 April. The backend filter banks have 256 channels with 1 MHz width and separation. At these frequencies, we have velocity resolution of ∼1 km s⁻¹ and ∼30'' resolution.

For each source, we made two arcminute square OTF maps with 10'' spacing. All sources were scanned in R.A. and decl. and most were also scanned again in R.A. and decl. with the observational scheme rotated by an arbitrary position angle or offset by 5'' from the original map center. Pointing was performed using Jupiter, Saturn, or Mars and the main beam efficiency was calibrated with Jupiter or bootstrapped from the intensity of S140 if Jupiter was not available. Pointing accuracy for the HHT is typically 5'' rms.

Spectra were reduced using GILDAS CLASS reduction software. The ALMA prototype receiver is a dual polarization receiver. We scale the horizontal polarization (H_pol) data to the vertical polarization (V_pol) data using the main beam efficiency determined for each feed. After a linear baseline was removed, the spectra were convolved onto a 11 × 11 10'' grid using a Gaussian-tapered Bessel function of the form $\frac{J_1 ({\frac{r}{a}})}{({\frac{r}{a}})} \mathrm{exp}- (\frac{r}{b})^2$ where $a = 1.55 (\frac{\mathrm{\theta _{mb}}}{3})$ and $b = 2.52 (\frac{\mathrm{\theta _{mb}}}{3})$ (Mangum et al. 2007). A Bessel function is a more accurate representation of the spatial frequency response of a single dish. Tapering with a Gaussian is a compromise to improve computation time, although the spatial frequencies lost in doing so are not significantly different from those probably already tapered by the response of the receiver feed horns (Mangum et al. 2007).

We also downloaded SCUBA 850 μm archived data (Di Francesco et al. 2008) for the 23 of our sources observed with SCUBA. We display a portion of the image centered on the continuum peak for each available source (see the Appendix).

We present maps of the integrated intensity of two wavelengths of dust continuum emission and seven molecule transitions. We calculate the integrated intensity from

$$\begin{equation} I(T_{\rm mb}) = \int T_{\rm mb} dv \;\; \pm \sqrt{\delta v_{{\rm line}} \delta v_{{\rm chan}}} \sigma _{T_{\rm mb}}, \end{equation} \tag{ 1 }$$

where δv_line is the full velocity width of the line, the width of the line above $\sigma _{T_{\rm mb}}$, and δv_chan is the channel width of the spectrometer. Emission maps of integrated intensity, I(T_mb), are shown in the Appendix. All 27 sources are detected in the mapped transitions of HCO⁺, N₂H⁺, HNC, CCH, and SO but not every one is detected in the rarer isotopologues (H¹³CO⁺, H¹³NC, ³⁴SO) nor the transitions of HCS⁺, CH₃OH, and SO₂. The HCO⁺ J = 3 → 2 transition typically has the highest integrated intensity with a median intensity, I(T_mb) = 61.5 K km s⁻¹. Figure 2 shows histograms of the integrated intensity for every transition we observed. Observed properties for the mapped transitions are listed in Tables 3–8.

Zoom In Zoom Out Reset image size

Most transitions trace a single emission peak centered near the dust continuum peak position. The emission maps of N₂H⁺ are significant exceptions. For 10 of our sources (W3(OH), RCW142, M8E, G8.67, W43S, W44, ON2S, W75(OH), S140, and CepA), we find that N₂H⁺ emission shows significant spatial differentiation that is not reflected in the dust emission or emission from other molecules. These 10 sources not only have multiple intensity peaks in the N₂H⁺ emission, but these peaks are also spatially separated from the peak of the dust and other molecule emission. In general, the dust and other molecules peak near the water maser position while the N₂H⁺ emission in the differentiated sources peaks several arcseconds (∼15'') away from the water maser position—more than can be accounted for with pointing errors alone. Both N₂H⁺ J = 3 → 2 and J = 1 → 0 show this differentiation, indicating that this is a spatial variation in the abundance of N₂H⁺ and not a radiative transfer effect.

The spatial variation in the N₂H⁺ emission may be understood from the dependence of N₂H⁺ chemistry on temperature. N₂H⁺ chemistry and HCO⁺ chemistry are directly related by the gas phase abundance of CO. HCO⁺ is created by CO while N₂H⁺ is more abundant when CO is frozen out onto dust grains (see Jørgensen et al. 2004). Each source with spatial variations in the N₂H⁺ emission also has centrally peaked HCO⁺ emission, supporting the conclusion that these sources are chemically differentiated. Warmer temperatures (T > 20 K) are required to return CO to the gas phase, suggesting that N₂H⁺ is preferentially differentiated toward warmer, more evolved sources.

To test this hypothesis, we looked at several physical properties that might indicate the presence of at least one evolved protostar embedded in the clump. The isothermal dust temperature, T_iso, defined as the temperature required to match the mass of the best-fit dust continuum radiative transfer model (Mueller et al. 2002),

$$\begin{equation} T_{{\rm iso}} = 41 ({\rm K}) \left[ \ln \left(\frac{M_{{\rm iso}}}{5.09 \times 10^{-8} \,M_{\odot } S_{\nu } ({\rm Jy}) D^2 ({\rm pc})} + 1 \right) \right] ^{-1} \;\;\;, \end{equation} \tag{ 2 }$$

should be higher for warmer, more evolved clumps. The spread in T_iso is within a few K of the median T_iso (22 K), reflecting the unresolved temperature gradients found for these sources by Mueller et al. (2002; see Table 2). Warmer clumps will also be more luminous, although more massive clumps will have larger luminosities even at lower temperatures. We make a simple correction for clump mass by taking the ratio of clump luminosity to the dust-determined isothermal mass, M_iso, defined as

$$\begin{equation} M_{{\rm iso}} = \frac{S_{\nu } D^2}{B_{\nu } \kappa _{\nu }} = 5.09 \times 10^{-8} \,M_{\odot } S_{\nu }({\rm Jy}) D^2({\rm pc})\big(e^{41K/T_{{\rm iso}}} - 1\big), \end{equation} \tag{ 3 }$$

where S_ν is the source dust flux density in a 120'' aperture (Mueller et al. 2002). The median L_bol/M_iso value for the differentiated sources is somewhat larger than for the remaining sample ($110 \,{\pm}\, 58 \frac{L_{\odot }}{M_{\odot }}$ compared with $50 \,{\pm}\, 23 \frac{L_{\odot }}{M_{\odot }}$). The difference suggests that the differentiated clumps are warmer on average, although this conclusion is not robust.

If the hypothesis that N₂H⁺ is differentiated with respect to the dust in more evolved clumps is correct, then the ionizing radiation we expect from more evolved protostars should lead to radio continuum (H ii) emission from the clump, provided there is a protostar of early enough type (>B0). Using the Shirley et al. (2003) clump classifications (see Table 2), we find that eight of the strongly differentiated sources are associated with H ii regions, indicating the presence of a more evolved massive protostar (see Table 2). Care must be taken in interpreting the evolutionary state of the clump as a whole as recent interferometric work has shown that many of these regions have multiple cores in different evolutionary states that are unresolved in the single dish beam (e.g., Hunter et al. 2008; Brogan et al. 2009).

Strong outflows in sources with N₂H⁺ differentiation may shock heat gas in the clump, leading to increased N₂H⁺ destruction. Most observations of massive star-forming regions, including those studied here, are observed with resolution insufficient to determine the outflow structure in the clump. It is possible that a single protostellar object powers some outflows, though others may reflect the combined effect of outflows from multiple protostars in the clump (Bally 2008). Interferometric follow-up is required to directly test the effect of outflows on N₂H⁺ chemistry.

The majority of our sources show N₂H⁺ peak emission that is spatially coincident with the dust peak at 30'' resolution, indicating that emission within the central beam is dominated by cold (T < 25 K), dense (n > 10⁴ cm⁻³) gas. However, this does not mean that N₂H⁺ is not chemically differentiated in these sources. Unresolved chemical differentiation may be reflected in the poor correlation between N₂H⁺ emission and emission from the dust and other molecules. We explore these correlations in Sections 3.2 and 3.3.

Spatial differentiation in N₂H⁺ emission has been observed toward other high-mass clumps. Pirogov et al. (2003, 2007) studied 35 dense molecular cloud cores known to be forming massive stars or clusters of stars and found N₂H⁺ differentiation in ${\sim} \frac{1}{3}$ of their sources. This differentiation is not what we would expect based on studies of N₂H⁺ emission in less evolved, low-mass star-forming regions, where N₂H⁺ emission is well correlated with the dust emission (Lee et al. 2001; Crapsi et al. 2005; Emprechtinger et al. 2009). For a single core, N₂H⁺ differentiation indicates a warmer source that is likely in a more advanced evolutionary state. However, as Pirogov et al. point out, in high-mass star-forming regions, emission from cores in several different evolutionary states overlap in the single dish beam, leading to a more complicated emission profile. In addition, more massive stars evolve faster, so emission from the hot region surrounding a massive young stellar object (YSO) may drive the observed differentiation while we simply cannot spatially resolve the colder, less evolved regions where N₂H⁺ is still well coupled to the dust. Dynamical activity in high-mass clumps may generate additional pressure on an individual collapsing core causing the collapse to proceed faster than free fall, thus producing high densities before the molecules responsible for N₂H⁺ destruction (e.g., CO) are frozen out (Lintott et al. 2005). Lintott et al. (2005) suggest that preventing freeze out facilitates N₂H⁺ destruction near massive YSOs and enhances emission from other species that freeze out in more quiescent low mass cores (Myers & Benson 1983; Evans et al. 2001).

Most of our sources do not show large-scale spatial differentiation in any transitions other than the two N₂H⁺ transitions, with the notable exception of ON2S. The 850 μm continuum emission for ON2S shows two peaks in the larger, more southern source. The northeast peak is brighter than the southwest peak in both 350 μm and 850 μm emission indicating that the northwest source is warmer. HCO⁺, CS, and HNC emission appears to trace the more northern of the two peaks in the 850 μm continuum, while the position of the N₂H⁺, SO, and CCH emission is more consistent with the position of the more southern peak (see Figure 12 in the Appendix). A colder southern source with strong N₂H⁺ emission is consistent with the predictions for N₂H⁺ chemical differentiation, however the chemical differentiation in SO and CCH is not easily explained. Higher spatial resolution interferometric observations are required to resolve the individual northeastern and southern sources. In a future work, we will use the Submillimeter Array to study the temperature and the chemistry of these sources individually in order to better understand what may be causing the unusual SO spatial differentiation in this source (M. Reiter et al. 2011, in preparation).

Two of the mapped molecules in this survey, CCH and SO, are molecules with unpaired electrons and therefore sensitive to the Zeeman effect and potential tracers of the magnetic fields (e.g., Bel & Leroy 1989, 1998). Zeeman splitting is too small to be observed with our spectral resolution for the transitions considered. The similarity of the CCH and SO emission to the dust continuum emission in this survey indicates that observations of these Zeeman tracers at higher resolution should typically be centered within 15'' of the dust continuum position. Ultimately, observing magnetic fields in these clumps via Zeeman splitting in lower frequency transitions of these molecules will best be done with the sensitivity and high spatial resolution available with ALMA (e.g., 3 mm band for CCH and 9 mm band for SO).

3.1.1. Molecular Tracers with Non-detections

We quantify how well the molecules trace the dust (and each other) by calculating the correlation between peak molecular emission and the dust emission at both 350 μm and 850 μm (correlations with 350 μm dust emission are plotted in Figure 3). We use the Bayesian IDL routine LINMIX_ERR (Kelly 2007) to perform a linear regression to find the slope of the best fit line. The relationship between dust emission and molecule emission may not be strictly linear, however, given the scatter in the correlation plots (see Figure 3), we are not justified in fitting a more complex function. We also calculate the Spearman rank coefficient for every combination of dust and molecule intensity. The Spearman rank analysis determines the degree of correlation between two distributions nonparametrically, and is therefore less sensitive to outliers in the data. Larger values of the correlation coefficient (r_corr ⩾ 0.7) indicate a more robust correlation.

Zoom In Zoom Out Reset image size

Correlations of each molecule with the 350 μm and the 850 μm dust emission are in excellent agreement (see Columns 1 and 2, respectively, of Table 9), which is encouraging given that the dust observations have different resolutions (15'' versus 10''), and were taken with different telescopes using different observing techniques (jiggle mapping versus scanning). High n_eff tracers CCH, HNC, CS, and SO are strongly correlated (r_corr ⩾ 0.7) with both 350 μm and 850 μm dust emission. In general, transitions with higher n_eff do a good job of tracing the dust continuum flux density.

In contrast, both transitions of N₂H⁺ show the poorest correlation with the dust emission (r_corr ⩽ 0.5), in agreement with the Zinchenko et al. (2009) result that N₂H⁺ does not trace the dust emission in high-mass star-forming regions. As discussed in the previous section, this is likely due to the destruction of N₂H⁺ in the warm, CO-rich gas in the center of these sources, unresolved within the 30'' beam. This divergent N₂H⁺ chemical behavior is more apparent when compared to the other molecules.

We also inter-compare the integrated intensities for each molecule transition (Columns 4 through 10 in Table 9). While a detailed understanding of these molecule–molecule correlations requires multi-dimensional, clumpy radiative transfer modeling with realistic (and probably complicated) molecular abundance profiles, several striking trends are apparent among the correlations. Molecular tracers with n_eff >10⁵ cm⁻³ (CCH, HNC, CS, and SO) have integrated intensities that are generally well correlated (r_corr ⩾ 0.7) with each other. CCH and HNC (n_eff = 1.5 × 10⁵ cm⁻³ and n_eff = 1.7 × 10⁵ cm⁻³, respectively) are strongly correlated with every transition. The two densest gas tracers, CS and SO, are poorly correlated with N₂H⁺, again supporting the hypothesis that N₂H⁺ is destroyed in the warm gas in the center of these clumps (see Section 3.1).

CS and SO have somewhat weaker correlation coefficients overall compared to CCH and HNC even though CS and SO are the highest n_eff tracers in this survey. This may be due to sulfur's complicated chemistry, which, in part, remains poorly understood due to the difficulty observing S⁺ or S ii. Beuther et al. (2009) studied the chemistry of sulfur in molecules in their study of high-mass protostellar chemical evolution. As discussed in that work, CS is released from grains early (CS desorbs from grains at temperatures of ∼30 K), then reacts with OH to form SO and eventually SO₂. It is possible that shock heating from strong outflows may also drive up gas temperatures, facilitating this chemical pathway. With higher spatial resolution, such as that used in Beuther et al. (2009), we may be able to use the relative abundances of CS, SO, and SO₂ to determine the chemical evolutionary state of the individual cores in our sources. Without resolving individual cores, we can say that this kind of chemical evolution likely contributes to the weaker correlations with S-bearing molecules. CS- and SO-integrated intensities are only modestly correlated (r_corr = 0.67) as we would expect for unresolved chemical differentiation due to the destruction of CS in the production of SO.

Tracers with n_eff ≈10⁴ cm⁻³ (HCO⁺ and N₂H⁺ J = 3 → 2) also show relatively weak correlations with CS and SO (except for HCO⁺ with SO), though they are well correlated with the other transitions. Given these poor correlations and the fact that N₂H⁺ is the only molecule to show substantial differentiation at 30'' resolution, we conclude that the N₂H⁺ abundance structure within the central telescope beam is significantly different than for the other molecule species. These poor correlations suggest that N₂H⁺ is chemically differentiated with respect to the other molecules in many of our sources, but this differentiation is spatially unresolved for the majority of the sources in this survey.

To compare the behavior of the molecules, we calculate the column density. Molecules with similar n_eff will be excited in the same gas, so we expect that as one transition gets brighter, so should other transitions with similar n_eff, in the absence of strong abundance gradients. Poor correlations may indicate chemical differentiation or the development of substructure (cores within the clump), even if the scales on which these things are occurring are too small to be resolved with the single dish beam.

We calculate the column density from

$$\begin{equation} N_{ul} = \frac{4\pi Q(T_{\rm ex}) e^{E_u/kT_{\rm ex}} B_{\nu }(T_{\rm ex})}{hcg_u A_{ul} [J_{\nu }(T_{\rm ex}) - J_{\nu }(T_{{\rm cmb}})] } I(T_{\rm mb}) \frac{\tau }{1-e^{-\tau }} \;\; \rm {cm}^{-2}, \end{equation} \tag{ 4 }$$

where g_u is the statistical weight for an energy level E_u, A_ul is the transition Einstein A, B_ν is the Planck function, J_ν(T) = (hν/k)/[exp (hν/kT) − 1] is the Planck function in temperature units (K), Q is the partition function for the molecule in LTE with excitation temperature T_ex, and τ/(1 − e^−τ) is the correction for non-zero optical depth. We approximate the average excitation temperature along the line of sight with the isothermal dust temperature T_iso (see Section 3.1). T_iso averages over significant temperature gradients in the clump, and may not represent the temperature of the gas where a specific transition is excited. However, for the range of T_iso for these clumps, the column density is not very sensitive to the differences in temperature (see Figure 4), changing by at most a factor of two.

Zoom In Zoom Out Reset image size

We attempt to correct for optical depth, τ, for those molecules for which we have observed a rarer isotopologue. With a value for τ, we can correct the column density for optical depth by multiplying the optically thin column density with a correction factor, τ/(1 − e^−τ) (Goldsmith & Langer 1999). If we assume that the isotopologue (H¹³CO⁺ for example) is optically thin, then we can use the ratio of the peak temperature of the rarer isotopologue (H¹³CO⁺) to the peak temperature of the main isotopologue (HCO⁺) to solve for the optical depth. To do so, we assume that the molecular abundance of H¹³CO⁺ relative to HCO⁺ is equal to the ISM abundance of atomic ¹³C relative to ¹²C ([¹³C]/[¹²C] ${=}\, \frac{1}{77}$; Wilson & Rood 1994). With this assumption and ignoring background radiation, we can calculate the optical depth τ as

$$\begin{equation} \frac{[T_{\mathrm{H^{13}CO^+}}]}{[T_{\mathrm{HCO^+}}]} = \frac{1-e^{-\tau _{\mathrm{H^{13}CO^+}}}}{1 - e^{-\tau _{\mathrm{HCO^+}}}} = \frac{1 - e^{-X \tau }}{1 - e^{-\tau }}, \end{equation} \tag{ 5 }$$

where X is the ratio of [¹³C]/[¹²C], τ is the optical depth of the HCO⁺ line, and H¹³CO⁺ is assumed to be thin. For HNC, where both isotopologues show linewidths that are approximately Gaussian, we find a median τ = 3.26 ± 0.72. We know that HCO⁺ J = 3 → 2 must be optically thick in at least some of our sources because we see self-absorption in the line profiles. Measuring the temperature of a self-absorbed line is difficult and leads to considerable uncertainty in the temperature of the main isotopologue. ISM abundance ratios for [¹³C]/[¹²C] and [³⁴S]/[³²S] are also uncertain to at least 10%, further increasing our τ uncertainty. Not surprisingly, our highly uncertain optical depth correction dramatically increases the scatter in our column densities, and thus decreases the correlation in many cases. We have calculated correlation coefficients for both the optically thin and the optically thick column density and report both in Tables 10 and 11, respectively.

We calculate total H₂ column density from the 350 μm dust emission or the 850 μm dust emission for sources without 350 μm data from

$$\begin{equation} N_{\mathrm{H}_2} = \frac{S_{\nu }^{{\rm obs}}}{B_{\nu }(T_d) \kappa _{\nu } \mu m_{{\rm H}} \Omega _{{\rm beam}}}, \end{equation} \tag{ 6 }$$

where S^obs_ν is the observed flux in a 30'' aperture and κ_ν = 0.10 cm² g⁻¹ is the opacity of gas and dust at 350 μm, and κ_ν = 0.0197 cm² g⁻¹ at 850 μm (Ossenkopf & Henning 1994). We use T_iso as the dust temperature and assume a mean molecular weight of μ = 2.34 as is appropriate for dense molecular gas. We plot the correlation between the molecular column densities and H₂ column density in Figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Correlations between the H₂ column density and the optically thick molecule column densities follow the same general trends as correlations of the molecule-integrated intensities with the dust emission, although this is expected since N_ul ∼ I(T_mb) in the optically thin limit. Indeed, the most pronounced differences are for correlations of optical-depth corrected column densities with the H₂ column density. These are made notably worse by the optical depth correction, dropping from strongly correlated (r_corr = 0.78 for the integrated intensity correlation) to relatively poorly correlated (r_corr = 0.54 for the column densities) in the case of HNC.

Zoom In Zoom Out Reset image size

Our suite of molecular transitions represents a range of excitation conditions, probed by n_eff, allowing us to examine how the derived clump properties change as a function of the gas that a specific transition is tracing. In this section, we examine how size R, linewidth Δv, virial mass M_vir, mass surface density Σ, and average density 〈n〉 vary with n_eff (see Figure 8 and Table 12).

Zoom In Zoom Out Reset image size

Table 12. Median Source Properties

Transition	log(neff)	R	Σ	log 〈n〉	Mvir	Δv	L'	$\frac{L_{{\rm bol}}}{M_{{\rm vir}}}$
	(cm−3)	(pc)	($\frac{\mathrm{g}}{\mathrm{cm^2}}$)	($\frac{\mathrm{g}}{\mathrm{cm^3}}$)	(M☉)	($\frac{\mathrm{km}}{\mathrm{s}}$)	($\frac{\mathrm{Jy}}{\mathrm{pc}^{2}} \frac{\mathrm{km}}{\mathrm{s}}$)	($\frac{L_{\odot }}{M_{\odot }}$)
N2H+ 1–0	3.30	0.29	0.17	4.66	308	2.77	58.30	168.42
HCO+ 3–2	4.36	0.25	0.36	5.12	640	6.50	24.35	74.57
N2H+ 3–2	4.43	0.26	0.73	5.28	639	4.14	11.41	51.08
CCH 3J–2J	5.18	0.25	0.55	5.28	524	3.93	11.51	96.06
HNC 3–2	5.23	0.24	0.35	5.19	446	4.86	9.23	61.06
CS 5–4	5.57	0.20	0.53	5.23	700	4.81	11.32	154.05
SO 67–56	6.38	0.17	1.43	5.82	689	5.09	2.37	102.09

Download table as:  ASCIITypeset image

Our sample of high-mass star-forming regions represents a range of masses, chemistries, and evolutionary states. These source-to-source variations are evident in the scatter in the relationships between molecule emission and dust emission (see Sections 3.2 and 3.3 and Figures 3–7). These variations also lead to a broad range of values for the physical properties derived from this sample (though they tend to be positively skewed). We therefore take the median value of each physical property and examine its behavior as a function of n_eff and characterize the spread of values around that median with the median absolute deviation. Because we analyze trends observed in a population of 27 sources, we emphasize that these trends are primarily qualitative as the uncertainties on the median values are large. However, we note that trends in the median value of a given physical property that persist despite the substantial scatter in the values for individual sources are remarkable given the diversity of our sources and relatively small sample size.

4.1.1. Size

Clump size, as determined from a molecular intensity map, depends on the integrated intensity threshold chosen. In this paper, we use the half-peak intensity contour though other definitions that depend on a fixed intensity level are possible. Comparisons between the different definitions generally give the same qualitative conclusions (see Shirley et al. 2003). We find the deconvolved angular radius of a circle with the same area as is enclosed within the half-peak intensity contour,

$$\begin{equation} \theta = \sqrt{ \frac{A}{\pi } - \frac{\theta ^2_{\mathrm{beam}}}{4}}, \end{equation} \tag{ 7 }$$

where A is the area enclosed in the contour in arcsec² and θ_beam is the beam size in arcsec. The angular radius can be converted to a linear radius using the distance, R = θD.

The observed clump size also depends on the underlying density structure of the clump. The large-scale structure of the clumps can be described, on average, with a power-law density distribution (n∝r^−p; Mueller et al. 2002; see also Shirley et al. 2003). The convolution of a Gaussian telescope beam pattern with a power-law intensity distribution results in an FWHM core size slightly larger than the beamsize with the observed FWHM larger for flatter power-law indices, p (see Shirley et al. 2003). To account for this, we convolved the N₂H⁺ J = 1 → 0 data to 30'' resolution (originally 18'') in order to study all the mapped transitions at the same resolution. With these caveats in mind, we analyze the distribution of sizes for each core in the mapped transitions.

An anti-correlation is observed between the median clump size and the effective density of the tracer (Figure 8). This is the expected trend for a power-law density distribution; the higher density tracers are preferentially excited in a smaller emitting region within the clumps. A similar trend was observed in the HCN and CS multi-transition study of Wu et al. (2010); however, their study contained observations with very different spatial resolutions, from 60'' to 30'', while the observations in this paper are all similar (30'').

We also see a trend when we look at clump sizes as a function of n_eff for each evolutionary state. If we separate sources into two categories, one for less evolved sources (weak/no radio continuum and UCH ii) and one for more evolved sources (CH ii and H ii), we find that more evolved sources are ∼0.1 pc larger than less evolved sources, although the difference in size as a function of evolutionary state becomes less pronounced for higher n_eff where source sizes get closer to the beam size. We note that this is a qualitative trend as this analysis requires separating 27 sources into four categories. With so few sources entering into the medians, the median absolute deviation is dominated by small number statistics.

4.1.2. Masses

Assuming these clumps are self-gravitating, we calculate the virial mass from

$$\begin{equation} M_{{\rm vir}} = \frac{5 R \Delta v^2}{8 a_1 a_2 G \ln (2)} \approx 209 \,{M}_{\odot } \;\; \frac{(R / 1 \,\mathrm{pc}) (\Delta v / 1 \,\mathrm{km \,s^{-1}})^2 }{ a_1 a_2 } \end{equation} \tag{ 8 }$$

$$\begin{equation} \begin{array}{cc}a_1 = \displaystyle\frac{1 - p/3}{1 - 2p/5}, & p < 2.5, \\ \end{array} \end{equation} \tag{ 9 }$$

where a₁ is the correction for power-law density distribution and a₂ is the correction for non-spherical shape (Bertoldi & McKee 1992). For aspect ratios less than 2, a₂ ≈ 1, so we do not include its effects here. We use power-law indices p found in Mueller et al. (2002) to calculate a₁. We assume a value of a₁ = 1.39 (corresponding to the median power-law index, p = 1.75) for sources without 350 μm modeling.

In the M_vir calculation, the linewidth, Δv², reflects the thermal energy and turbulent pressure supporting the cloud against gravity. Linewidths may be broadened by optical depth and systematic collapse motions (for example), leading to an artificially large estimate of the mass. Observing a rarer isotopologue of the molecule mitigates the effects of optical depth in some cases, as the optically thin linewidth from the rarer isotopologue may be substituted for the optically thick linewidth of the main isotopologue. Other systematic motions can still broaden the linewidth and, in many cases, we cannot correct for those effects. Thus, any mass estimate dependent on the linewidth is only an upper bound.

We have observed isotopologues of four species in this survey, H¹³CO⁺, HN¹³C, C³⁴S, and ³⁴SO, and have derived reliable linewidths for the first three. Comparing the median M_vir calculated with the main and the rarer isotopologue linewidths, we find that the broader linewidths of the main isotopologues lead to a systematic overestimation of the mass by a factor of up to 1.83 (see Figure 9).

Zoom In Zoom Out Reset image size

Δv is also not correlated with n_eff (see Figure 8). For every source with a reliable fit for the linewidth, the optically thin isotopologue has a smaller linewidth than the main isotopologue, although this difference is most pronounced for HCO⁺, which we expect to be optically thick (see Figure 10). No correlation emerges when we compare the optically thin isotopologues with n_eff, we find only that the scatter in the lack of correlation is somewhat decreased.

Zoom In Zoom Out Reset image size

4.1.3. Mass Surface Density and Volume Density

McKee & Tan (2003) first noted that almost all regions forming massive stars have mass surface densities Σ ∼ 1 g cm⁻² (equivalent to a column density of >2 × 10²³ cm⁻²). Recent numerical simulations by Krumholz et al. (2007a) and Krumholz & McKee (2008) have shown that Σ ⩾ 1 g cm⁻² is a minimum threshold to prevent fragmentation of a massive clump into several smaller objects. Σ of order 1 corresponds to an opacity high enough for gas to be heated substantially (hundreds of K) and to substantial distances (⩾1000 AU) from the protostar (Krumholz et al. 2007a). Warmer temperatures are crucial to prevent the clump from fragmenting into smaller objects, as the Jeans Mass decreases for higher densities for a fixed temperature (M_J ∼ T^3/2).

We calculate the mass surface density of our sources using

$$\begin{equation} \Sigma = \frac{M_{{\rm vir}}}{\pi R^2} \approx 0.665 \frac{(M_{{\rm vir}} / 1.0 \times 10^4 \,M_{\odot })}{(R / 1 \,\mathrm{pc})^2} \,\mathrm{g \,cm^{-2}}. \end{equation} \tag{ 10 }$$

The median surface density for each transition is positively correlated with n_eff. This is as expected; higher n_eff transitions trace denser regions of the clump, corresponding to larger surface densities. Surface densities are >1 for ${\sim} \frac{1}{3}$ of our sources when calculated from transitions with n_eff ≈ 10⁴ cm⁻³. When measured with a tracers with n_eff above a few 10⁶ cm⁻³, ${\sim} \frac{2}{3}$ of the clumps have Σ ⩾ 1.0 g cm⁻². Unlike median size, median surface density shows no apparent trend with evolutionary state.

Krumholz & McKee (2008) point out that even though suppression of fragmentation allows massive stars to form from massive cores, some low-mass stars will form before radiation halts fragmentation. Therefore, massive stars are expected to be found in clusters. Volume-averaged density, 〈n〉, for a clump with clustered star formation embedded within will be dominated by the diffuse gas between the high density cores, leading to smaller values of 〈n〉 than would be expected for an individual star-forming core (Dunham et al. 2010). 〈n〉 is then a valuable tool for determining whether objects discovered in new and upcoming surveys (e.g., BGPS) are individual star-forming cores or clumps with cores embedded within. For this survey, all of our sources are clumps rather than individual cores, so we focus on characterizing the variation in 〈n〉 with n_eff.

Assuming these clumps are spherical, we can calculate the volume-averaged density,

$$\begin{equation} \mbox{$\langle n\rangle $} = \frac{M_{{\rm vir}}}{\frac{4}{3} \pi R^3 \mu m_{{\rm H}}}, \end{equation} \tag{ 11 }$$

where μ is the mean molecular weight of molecular gas (μ = 2.34), m_H is the mass of hydrogen, and R is the radius of the half-peak intensity contour. We find that 〈n〉 increases with n_eff, exactly as we expect for transitions with higher n_eff tracing denser gas (see Figure 7). As a function of evolutionary state, we find that the 〈n〉 values for sources associated with no H ii emission or (U)CH ii emission are essentially the same, while 〈n〉 for sources associated with H ii regions tends to be higher. The difference is small and may simply reflect small number statistics.

We find poor correlations (r_corr < 0.7) between column density and the physical variables L_bol, M_iso, 〈n〉, and Σ. Several factors contribute to these poor correlations. A specific molecular transition will trace out the parts of a clump with conditions sufficient to excite that transition. Clumps likely contain unresolved cores that, if they are forming stars, have strong gradients with temperature and density increasing toward the center of the core. Chemical variability from source to source as well as within our beam increases the scatter in comparisons of the column densities.

Qualitatively, the column density of higher n_eff tracers tend to be larger for larger L_bol and M_iso, though the correlation never becomes robust (r_corr > 0.7). Both volume-averaged density, 〈n〉, and mass surface density, Σ, are completely uncorrelated with molecular column densities regardless of n_eff. Ratios of column densities also fail to correlate with any of the physical variables. As discussed in the previous section, the various tracers in this study are sensitive to different excitation conditions and we believe that analyzing the median value of the clump properties is the best way to understand how derived clump properties depend on the tracer chosen.

Understanding clustered, massive star formation in our own galaxy is an essential step to understanding star formation in other galaxies. Except for the nearest galaxies, star formation can only be assessed via the unresolved star formation indicators across a galaxy. Beginning with Schmidt (1959), star formation prescriptions in other galaxies have relied on the correlation between the galaxy's gas content and its star formation rate (SFR). Kennicutt (1998) refined this relationship using H i and CO emission to measure the amount of atomic and molecular gas available for star formation, yielding a power-law index of 1.4–1.5. Molecular gas is essential to star formation, but much of the molecular gas traced by CO (n_crit = 2.8 × 10³ cm⁻³, n_eff = 30–50 cm⁻³) is probably diffuse, and therefore may not be directly involved in the star formation process (Gao & Solomon 2004).

More recent work has focused on higher density tracers such as HCN J = 1 → 0 (n_crit = 1.2 × 10⁶ cm⁻³, n_eff = 1.4 × 10⁴ cm⁻³) and HCN J = 3 → 2 (n_crit = 6.8 × 10⁷ cm⁻³, n_eff = 1.5 × 10⁵ cm⁻³) to study the dense molecular gas that is much more likely to be associated with star formation (Gao & Solomon 2004; Bussmann et al. 2008). The star formation law derived from studying dense molecular gas yields a power-law index much closer to unity. Wu et al. (2010) extended the HCN J = 1 → 0 sample to galactic cloud cores and found that the data were still well fit with a power-law index of 1. This led to the interpretation that the linear correlation between dense molecular gas and the SFR may be understood as the cumulative contribution of individual star-forming "units."

Extragalactic star formation studies have used the molecular line luminosity, L', to estimate the mass of dense gas available for star formation (see Juneau et al. 2009 for a summary of dense gas tracers). L' is defined as

$$\begin{equation} L^{\prime } = 23.504 \pi \Omega _{s \ast b} (1+z)^{-3} D_L^2 I(T_{\rm mb}), \end{equation} \tag{ 12 }$$

where D_L is the luminosity distance in Mpc, z is the redshift, and Ω_s*b is the solid angle of the source convolved with the beam in arcsec². As described in Mangum et al. (2008), two lines with the same brightness temperature and spatial extent will have the same L', such that ratios of L' indicate the physical conditions in the gas.

Empirically, L' is found to trace the mass of molecular gas,

$$\begin{equation} M_{{\rm dense}} = \alpha L^{\prime }_{{\rm mol}}, \end{equation} \tag{ 13 }$$

where α is a scaling factor that has traditionally been based on the H₂-mass-to-line-luminosity found in the Milky Way (Downes et al. 1993). A direct test of how well L' actually traces the mass of dense gas is not possible without an independent determination of the mass. For our galactic clumps, we use M_iso (see Section 3.1; Mueller et al. 2002), which is determined from the optically thin 350 μm dust flux in a 120'' aperture. M_iso is independent of the integrated intensity of line and independent of the source size determined from the map of that line, making this method of determining mass completely independent from L'.

Comparing L' with M_iso, we find that the two are highly correlated (r > 0.7; see Figure 11) for every transition. Using LINMIX_ERR to do a linear regression, we find power-law indices that are slightly less than 1 in every case. There is no trend in power-law index with n_eff. Even N₂H⁺, which shows strong chemical differentiation in the integrated intensity maps and the poorest correlations with the dust flux, has an L' that is well correlated with M_iso (with somewhat more scatter than the other transitions). All of the dense gas tracers we have observed, despite chemical and kinematical behavior that affects the other derived properties of these clumps, are linearly related to the mass of dense gas available for star formation. Power-law indices slightly smaller than 1 indicate that L' will tend to underestimate the mass, and we show that using higher n_eff tracers will not resolve this issue.

Zoom In Zoom Out Reset image size

Looking at L' versus M_vir for multiple transitions of HCN and CS, Wu et al. (2010) also found that L' is an excellent tracer of the dense gas. As those authors discuss, questions have been raised concerning how reliably HCN, CS, and other high n_eff molecules trace the mass, given the high densities required for excitation. The power indices near unity we find support the Wu et al. (2010) conclusion that L' is linearly related to the mass regardless of n_eff of the tracers in these two studies because they trace only the dense gas most relevant to star formation.