Relations Between Molecular Cloud Structure Sizes and Line Widths in the Large Magellanic Cloud

TL;DR

This study investigates how the size-line width relation in molecular clouds within the Large Magellanic Cloud varies with star formation activity and mass surface density, highlighting the roles of gravity and feedback in cloud dynamics.

Contribution

It provides new observational evidence linking line width variations to star formation activity and surface density in LMC molecular clouds, extending previous studies.

Findings

Line width increases with star formation activity.

Line width correlates with mass surface density.

Both gravity and feedback influence cloud dynamics.

Abstract

We present a comparative study of the size-line width relation for substructures within six molecular clouds in the Large Magellanic Cloud (LMC) mapped with the Atacama Large Millimeter/submillimeter Array (ALMA). Our sample extends our previous study, which compared a Planck detected cold cloud in the outskirts of the LMC with the 30 Doradus molecular cloud and found the typical line width for 1 pc radius structures to be 5 times larger in 30 Doradus. By observing clouds with intermediate levels of star formation activity, we find evidence that line width at a given size increases with increasing local and cloud-scale 8${\mu}$m intensity. At the same time, line width at a given size appears to independently correlate with measures of mass surface density. Our results suggest that both virial-like motions due to gravity and local energy injection by star formation feedback play important roles in determining intracloud dynamics.

Tables

Table 1. Table 1: Summary of Map Parameters at a Common Resolution of 3 . ′ ′ arcsecond \farcs 5

Region	Transition	Ref. R.A.	Ref. Dec.	Mosaic Size	$\Delta v_{\mathrm{ch}}$	$T_{\mathrm{rms},12}$aarms noise in a channel map of width $\Delta v_{\mathrm{ch}}$.	$T_{\mathrm{peak},12}$bbPeak brightness temperature in cube.	$T_{\mathrm{rms},13}$	$T_{\mathrm{peak},13}$
		(J2000)	(J2000)	$\left(″\times ″\right)$	(km s-1)	(K)	(K)	(K)	(K)
30 Dor	$J$=2–1	5h38m47$\stackrel{s}{.}$0	$-69$°04′36″	50 $\times$ 050	0.5	0.05	43.6	0.05	14.1
N59C	$J$=1–0	5h35m18$\stackrel{s}{.}$8	$-67$°36′12″	160 $\times$ 260	0.2	0.19	26.9	0.14	06.2
A439	$J$=1–0	5h47m26$\stackrel{s}{.}$1	$-69$°52′46″	150 $\times$ 165	0.2	0.17	13.2	0.11	04.8
GMC104	$J$=1–0	5h21m04$\stackrel{s}{.}$0	$-70$°13′29″	150 $\times$ 160	0.2	0.21	17.4	0.10	06.4
GMC1	$J$=1–0	4h47m30$\stackrel{s}{.}$8	$-69$°10′32″	320 $\times$ 130	0.2	0.19	14.6	0.12	05.2
PCC	$J$=2–1	5h24m09$\stackrel{s}{.}$2	$-71$°53′37″	80 $\times$ 220	0.2	0.13	10.4	0.14	03.9

Table 2. Table 2: Power Law Fit Parameters: log ⁡ σ v = a 1 ​ log ⁡ R + a 0 subscript 𝜎 𝑣 subscript 𝑎 1 𝑅 subscript 𝑎 0 \log\sigma_{v}=a_{1}\log R+a_{0}

	12CO										13CO
Cloud	$a_1$				$a_0$				${\chi }_{\nu }^2$	$\epsilon$aarms scatter in $\log {\sigma }_v$ relative to the best-fit line. Units are dex.	$a_1$		$a_0$	${\chi }_{\nu }^2$	$\epsilon$aarms scatter in $\log {\sigma }_v$ relative to the best-fit line. Units are dex.
30Dor	0\@alignment@align.60	0.04 $\pm$	0\@alignment@align.22	0.03	4\@alignment@align.30	0.10 $\pm$	0\@alignment@align.43	0.08	0.21	0.04	10\@alignment@align.86	0.14 $\pm$
N59C	0\@alignment@align.97	0.04 $\pm$	-0\@alignment@align.34	0.02	17\@alignment@align.63	0.19 $\pm$	1\@alignment@align.00	0.12	-0.35	0.04	11\@alignment@align.41	0.24 $\pm$
A439	0\@alignment@align.53	0.02 $\pm$	-0\@alignment@align.24	0.02	9\@alignment@align.09	0.15 $\pm$	0\@alignment@align.83	0.04	-0.38	0.02	3\@alignment@align.07	0.13 $\pm$
GMC104	0\@alignment@align.46	0.02 $\pm$	-0\@alignment@align.31	0.02	6\@alignment@align.77	0.12 $\pm$	0\@alignment@align.36	0.02	-0.34	0.02	4\@alignment@align.76	0.15 $\pm$
GMC1	0\@alignment@align.29	0.03 $\pm$	-0\@alignment@align.38	0.02	9\@alignment@align.59	0.16 $\pm$	0\@alignment@align.57	0.08	-0.50	0.04	5\@alignment@align.98	0.15 $\pm$
PCC	0\@alignment@align.32	0.03 $\pm$	-0\@alignment@align.34	0.03	6\@alignment@align.56	0.13 $\pm$	0\@alignment@align.44	0.10	-0.47	0.05	4\@alignment@align.19	0.14 $\pm$
All	0\@alignment@align.65	0.03 $\pm$	-0\@alignment@align.33	0.02	60\@alignment@align.28	0.25 $\pm$	0\@alignment@align.54	0.04	-0.31	0.03	38\@alignment@align.99	0.25 $\pm$

Equations

$$J(T_{\mathrm{ex}})=f_{\mathrm{bm}}^{-1}{T}_{12,\mathrm{pk}}+J(T_{\mathrm{cmb}}),$$

$$J(T)\equiv \frac{h\nu /k}{\exp (h\nu /kT)-1}.$$

$$T_{13}=f_{\mathrm{bm}}[J(T_{\mathrm{ex}})-J(T_{\mathrm{cmb}})][1-\exp (-{\tau }_{13})],$$

$$N{(}^{13}\mathrm{CO})=\frac{3k}{8{\pi }^{3}B{\mu }^2}\exp \left[\frac{hBJ(J+1)}{k{T}_{\mathrm{ex}}}\right]$$

$$\frac{T_{\mathrm{ex}}+hB/3k}{1-\exp (-h\nu /k{T}_{\mathrm{ex}})}\frac{\int {\tau }_{13}dv}{J+1},$$

$$\frac{N({\mathrm{H}}_2)}{N({}^{13}\mathrm{CO})}=3\times 10^6.$$

$$2\mathcal{T}+\mathcal{W}=2\left(\frac{3}{2}M{\sigma }_v^2\right)-\frac{3}{5}\frac{G{M}^2}{R}=0,$$

$$\frac{M_{\mathrm{CO}}}{M_{\odot }}=4.3X_2\frac{L_{\mathrm{CO}}}{\mathrm{K}\mathrm{km}{\mathrm{s}}^{-1}\mathrm{p}{\mathrm{c}}^2},$$

$$\log {\sigma }_v=a_1\log R+a_0,$$

$$\log {\sigma }_v[\mathrm{km}{\mathrm{s}}^{-1}]=0.65\log R[\mathrm{pc}]-0.33,$$

$${\Sigma }_{\mathrm{vir}}=\frac{M_{\mathrm{vir}}}{\pi R^2}\propto \frac{{\sigma }_v^2}{R}=v_0^2,$$

$${\Sigma }_{\mathrm{vir}}-\Sigma =\frac{20}{3\pi G}\frac{P_{\mathrm{ext}}}{\Sigma }.$$

Full text

galaxies: ISM — radio lines: ISM — ISM: molecules — Magellanic Clouds

††software: CASA (McMullin et al., 2007), astrodendro (http://www.dendrograms.org), Kapteyn (Terlouw & Vogelaar, 2015), Scipy (Jones et al., 2001), SCIMES (Colombo et al., 2015), Astropy (Astropy Collaboration et al., 2013).

1 Introduction

The physical conditions within giant molecular clouds establish the initial conditions for star formation, thus understanding the factors that determine molecular cloud properties is of major interest. A correlation between size and line width, of the form $\sigma_{v}\propto R^{\alpha}$ with $\alpha\approx 0.5$, has long been noted in samples of nearby molecular clouds (Larson, 1981; Solomon et al., 1987, hereafter S87). This correlation, hereafter referred to as the $R$–$\sigma_{v}$ relation, is usually interpreted as the result of turbulent motions in the interstellar medium on all scales (Mac Low & Klessen, 2004; Falgarone et al., 2009). It closely resembles the turbulent cascade with a power-law slope falling between the Kolmogorov (1941) and Burgers (1939) values for incompressible and highly supersonic turbulence, respectively (see also Falgarone et al., 1994; Brunt & Heyer, 2002; Kritsuk et al., 2013; Federrath, 2013). At the same time, a study of 13CO emission in the Boston University-FCRAO Galactic Ring Survey by Heyer et al. (2009) showed that the normalization of the relation, $v_{0}=\sigma_{v}/R^{1/2}$, exhibits a linear correlation with mass surface density, $\Sigma=M/\pi R^{2}$, across more than an order of magnitude in $\Sigma$. The data are consistent with a state of virial balance between gravity and turbulent motions, except that the observed normalization $v_{0}$ is about a factor of 2 too large. Subsequent work by Field et al. (2011) has suggested that the larger than expected $v_{0}$ may result from external pressure confinement, although to a lesser extent than has been inferred for clouds in the outer Galaxy (Heyer et al., 2001) or near the Galactic Center (Oka et al., 2001; Shetty et al., 2012). On the other hand, Ballesteros-Paredes et al. (2011) have interpreted the Heyer et al. (2009) result in terms of gravitational collapse near free-fall, which differs from the virial equilibrium prediction by a factor of $\sqrt{2}$ in $v_{0}$, and is thus roughly consistent with the GRS data. A third possibility is that errors in the measured or inferred cloud properties create the appearance of excess kinetic energy when in fact clouds are close to being virialized.

Since interstellar turbulence transfers energy across spatial scales, the $R$–$\sigma_{v}$ relation can be studied on scales much smaller than the full extent of molecular clouds (e.g., Myers, 1983), all the way down to the $\sim$0.1 pc scales at which the thermal contribution to the line width becomes significant (Goodman et al., 1998). Kinetic energy spectra derived from techniques such as principal component analysis (PCA; e.g. Brunt & Heyer 2002; Brunt 2003) or wavelet transforms (e.g., Ossenkopf & Mac Low, 2002) have been interpreted as requiring turbulence to be driven on large ($>$10 pc) scales. The energy can be provided by gas accretion onto the disk of the Galaxy (Klessen & Hennebelle, 2010) or to some degree by large-scale instabilies and spiral waves (Wada et al., 2002). Stellar feedback in the form of supernovae, winds or expanding H II regions also produce enough energy to explain the observed energy budget (for reviews, see e.g. Mac Low & Klessen, 2004; Krumholz, 2014; Klessen & Glover, 2016), and so there is ongoing debate about the astrophysical origin of the observed turbulent motions. Progress is slowed by the fact that a priori it can be difficult to identify the driving scale of turbulence because energy can cascade to both larger and smaller scales (Vestuto et al., 2003). Another complication arises from the heterogeneous nature of Galactic data sets: high spatial resolution data are generally limited to the nearest clouds, and for most larger samples of clouds, distance uncertainties and distance-related selection effects are quite significant (e.g., Traficante et al., 2018).

The Large Magellanic Cloud (LMC) serves as an ideal laboratory to study molecular cloud turbulence on a range of scales and across a range of galaxy environments. The spatial dynamic range achievable from ground-based observations is rivalled only by observations of the Milky Way and M31, whereas the proximity ($d\approx 50$ kpc; Pietrzyński et al. 2019) and low inclination ($i\approx 34\arcdeg$; van der Marel & Kallivayalil 2014) of the galaxy enables clouds and H II regions to be easily identified. With the resolution provided by the Atacama Large Millimeter/submillimeter Array (ALMA), we have been investigating molecular cloud substructure across the LMC. Contrasting two clouds with very different star formation activity at the same spatial resolution, Wong et al. (2017, hereafter Paper I) found that the GMC associated with the actively star-forming 30 Doradus region shows a factor of $\sim$5 higher line width at a given spatial scale than a quiescent GMC in the outskirts of the galaxy. In this paper, we extend our previous study by analyzing four additional LMC clouds showing intermediate levels of star formation activity.

2 Observations

2.1 Sample Selection

Using data from the Magellanic Mopra Assessment (MAGMA) survey, Wong et al. (2011) identified 450 regions of contiguous CO emission (“islands”) in the LMC via the CPROPS emission segmentation software (Rosolowsky & Leroy, 2006). Additional observations conducted in 2012–3 increased the MAGMA coverage of the LMC by $\sim$20%, and the MAGMA Data Release 3 (DR3) contains 472 “islands” defined using a slightly more stringent set of CPROPS parameters (requiring a 3.5$\sigma$ rather than 3$\sigma$ peak). To select the cloud sample for ALMA mapping, we examined the joint distribution of CO and 8$\mu$m intensities for the islands, as measured by the MAGMA and Spitzer SAGE (Meixner et al., 2006) programs respectively, considering only the 65 clouds with angular areas between 6.25 and 12.5 arcmin2 (Figure 1, left). This range of cloud sizes was chosen to maximize the achievable spatial dynamic range while still allowing the mosaic to be completed within a single ALMA observation. We aimed for a sample exhibiting high CO intensity**, spanning a wide range in 8$\mu$m intensity, and distributed broadly across the galaxy, while also being** complementary to already approved ALMA observations. Our final sample consisted of four clouds, which we label as GMC1 and GMC104 (designations from the NANTEN cloud catalogue of Fukui et al. 2008), A439 (designation from the island catalogue of Wong et al. 2011), and N59C (designation based on the associated H II region catalogued by Henize 1956). To summarize, our sample consists of bright CO-emitting clouds selected to have a radius of $\sim$25 pc and to exhibit a wide range of mid-infrared surface brightness. The locations of our sample clouds within the LMC, including the two additional clouds studied in Paper I, are indicated in Figure 1 (right).

2.2 ALMA Cycle 4 Data

Observations toward four molecular clouds in the LMC were obtained in ALMA Cycle 4 under project code 2016.1.00193.S (PI: Wong). Each cloud was observed in two frequency settings, the first covering the 12CO(1–0) line and the second covering 13CO(1–0), C18O(1–0), CS(2–1), and C34S(2–1). The velocity resolution after Hanning smoothing was 61 kHz ($\sim$0.16 km s*-1*) across a bandwidth of 59 or 117 MHz. Although the CS(2–1) line was detected in all four clouds, in this paper we consider only the 12CO and 13CO data. Observations of 12CO were conducted in both the 12m and 7m arrays, providing sensitivity to structures up to 45″ in size, while observations of 13CO were conducted in the 12m array only, which is sensitive to structures up to 30″ in size at the observing frequency of 110 GHz.

While we did not revisit or modify the system calibration, we re-imaged the calibrated visibilities using the CASA package (McMullin et al., 2007). For each spectral line, the visibility data were imaged together using the tclean task in CASA 5.0.0. The imaging grid was set to have square pixels of width 0.5″ and velocity channels spaced by 0.2 km s*-1*. The size of the imaging grid was 800 pixels square for all clouds besides GMC1, for which extending the grid to 1000 pixels in R.A. was necessary. The flexible visibility weighting scheme of Briggs (1995) was used. For most clouds we obtained satisfactory imaging results using the auto-multithresh procedure in tclean with the clark deconvolver. In brief, the procedure continuously updates the deconvolution mask based on the current residual image. The initial masking of the image is controlled by the parameters pbmask, sidelobethreshold and noisethreshold, which we set to 0.2, 2, and 3 respectively. Regions smaller than the beam are pruned from the mask (minbeamfrac=1) and the resulting mask is convolved with a Gaussian 4 times the synthesized beam size (smoothfactor=4) and clipped at 10% of the smoothed peak (cutthreshold=0.1). We also expand the mask to an enclosing low-level contour of 1.5$\sigma$ (lownoisethreshold=1.5) to allow recovery of fainter extended emission. We found slightly different tclean parameters worked best for generating the deconvolution mask for the 7m data (noisethreshold=5, minbeamfrac=0.5, smoothfactor=1; all other parameters left unchanged). Deconvolution proceeded down to a 1$\sigma$ stopping threshold.

For combining the 12m and 7m data for 12CO, a hybrid imaging procedure was adopted, in order to ensure that the larger field-of-view 7m data were properly deconvolved. The 12m and 7m data were separately deconvolved using the auto-multithresh procedure in tclean, as described above, and then the combined visibilities were imaged and deconvolved using the union of the 12m and 7m masks, with no further mask adjustment. This final deconvolution was performed using the multiscale deconvolver (Cornwell, 2008) with scales of 0, 4, and 12 pixels, a stopping threshold of 2$\sigma$, and a smallscalebias parameter of 0.6 (the default value).

For one cloud (GMC104), significant sidelobes persisted after applying our standard imaging procedure, and interactive masking needed to be applied to finalize the imaging. We found that this procedure resulted in a $\sim$50% increase in the total integrated CO flux, but a negligible impact on the properties of the structures that are measured following decomposition (§3.3). This is not surprising, given that the structure properties are primarily sensitive to the brightest emission peaks, where the deconvolved flux is well-constrained.

For purposes of comparison, we smoothed the cubes to a common resolution of 3$\farcs$5 (0.8 pc). The rms noise and peak detected signal within the maps at this resolution are summarized in Table 1. To reduce noise when computing flux spectra or moment maps, we apply a signal mask to the cube prior to integration in the spatial or velocity dimensions. The mask was obtained by starting at the 3.5$\sigma$ contour and expanding to the surrounding 2$\sigma$ contour, then extending the mask by 1 channel towards the blue- and redshifted ends of each spectrum. We refer to this signal mask as the “dilated” mask. We compare the ALMA-derived 12CO flux spectra with those obtained from the single-dish MAGMA survey in Figure 2. The MAGMA spectra were obtained by applying the ALMA mosaic gain to the MAGMA map and employing the same masking procedure that was applied to the ALMA cubes. Since the differences are comparable to the $\sim$20% calibration uncertainties in the much lower sensitivity MAGMA maps, we have not attempted to merge the interferometer and single-dish data. Note that we lack single-dish or ACA 7m data for the 13CO line, so the spectra compared in the lower panels of Figure 2 are both from the ALMA cube, but using different signal masks (a dilated mask derived from the 12CO cube and then transferred to the 13CO, and one derived from the 13CO cube directly). The agreement indicates that the achieved signal-to-noise ratio is high enough that the measured fluxes are not sensitive to the details of the masking approach.

2.3 Archival Data for 30 Dor and PCC

We include archival ALMA data presented in Paper I for comparison with the data from our new ALMA Cycle 4 observations. The archival data were obtained over somewhat smaller fields of view in the $J=2\rightarrow 1$ transition of 12CO and 13CO. They include total power (TP) but not 7m data.

Data for the 30 Dor-10 molecular cloud (hereafter “30 Dor”) were collected in ALMA Cycle 0 under project code 2011.0.00471.S (Indebetouw et al., 2013). The ALMA data, covering a field 50″ $\times$ 50″ in size, were combined with total power data from APEX using the feather task in CASA to recover large-scale flux. The native resolution of the datacubes was 2$\farcs$38 $\times$ 1$\farcs$54 for 12CO and 2$\farcs$47 $\times$ 1$\farcs$59 for 13CO with 0.5 km s*-1* channels. The cubes were then smoothed to a 3$\farcs$5 circular beam for comparison with the other clouds.

For the quiescent cloud PGCC G282.98$-$32.40 (Planck Collaboration et al., 2016), which we refer to as the “Planck Cold Cloud (PCC)”, ALMA Cycle 2 observations were obtained in 2014 and 2015 under project code 2013.1.00832.S (Paper I). The observations cover a region of 220″ $\times$ 80″ and include total power data from ALMA merged using the feather task in CASA. The native resolution of the datacubes was 1$\farcs$72 $\times$ 1$\farcs$19 for 12CO and 1$\farcs$81 $\times$ 1$\farcs$24 for 13CO with 0.2 km s*-1* channels. Again, the cubes were then smoothed to a 3$\farcs$5 circular beam for comparison with the other clouds.

By including both Band 3 and Band 6 observations within our sample, we are assuming that the $J=1\rightarrow 0$ and $J=2\rightarrow 1$ lines trace similar structures within GMCs, and thus a direct comparison of results obtained from the different lines is possible. This is reasonable given the relatively modest excitation requirements for both lines ($E_{10}=5.5$ K, $E_{21}=11$ K) and the relatively small (factor of $\lesssim$2) departures of the $I_{\rm CO(2-1)}/I_{\rm CO(1-0)}$ ratio from unity in both Galactic (e.g. Sakamoto et al., 1995; Nishimura et al., 2015) and extragalactic (e.g. Leroy et al., 2009) studies. In particular, Sorai et al. (2001) measure a luminosity-weighted CO(2–1)/CO(1–0) brightness temperature ratio of $\sim$0.9 across the LMC. We note, however, that variations in the 13CO(2–1)/13CO(1–0) intensity ratio are expected to be larger than variations in the CO(2–1)/CO(1–0) ratio, given the lower optical depth of the 13CO lines (e.g. Nishimura et al., 2015). To more rigorously test our assumption that $J$ level does not strongly affect structure properties would require a comparative structural analysis of $J=1\rightarrow 0$ and $J=2\rightarrow 1$ data for the same cloud, which we defer to a future paper.

2.4 IR and Dust Comparison Images

We use published mid-infrared and far-infrared images of the LMC to measure global characteristics of the GMCs in our sample. The Spitzer 8$\mu$m and 24$\mu$m images from the SAGE legacy program (Meixner et al., 2006) are employed as indicators of star formation activity, and longer wavelength Herschel imaging from the HERITAGE key program (Meixner et al., 2013) are employed to trace dust temperature and mass. Specifically, for 8$\mu$m emission we use the 2″-pixel mosaic of point source subtracted residual images (SAGE_LMC_IRAC8.0_2_resid.fits), while for 24$\mu$m emission we use the 2$\farcs$49-pixel mosaic (SAGE_LMC_MIPS24_E12.fits). We use the dust temperature ($T_{\rm dust}$) and dust column density ($N_{\rm dust}$) maps generated by Utomo et al. (2019) at 13 pc (53″) resolution from the HERITAGE observations. Because of their coarser resolutions, the $\lambda\geq 24\,\mu$m images are used only to describe the overall GMC properties presented in §3.5. To do this, the boundary of each GMC is defined by a dilated mask derived from the MAGMA data starting from 3.5$\sigma$ peaks and extended to the 2$\sigma$ edge, regridded to match the SAGE pixel grid. Within each GMC boundary we obtain the median 8$\mu$m and 24$\mu$m intensity and a mean value for $N_{\rm dust}$ and $T_{\rm dust}$.

As the resolution of the 8$\mu$m images is similar to that of the ALMA data, we also measure the mean 8$\mu$m intensity in CO-emitting structures at a resolution matched to the CO data. For this purpose, the SAGE residual image is convolved to 3$\farcs$5 assuming a native resolution of 2″, and then regridded to match each ALMA image.

3 Results

3.1 Line Intensity Maps

We present the integrated intensity maps of 12CO and 13CO emission for the full sample of six clouds in Figures 3 and 4. These were obtained by applying the dilated masks described in §2.2 to the cubes and summing in velocity. Since the mosaic field of view (indicated by the red dashed line) differs across the sample, the FWHM size of the common-resolution beam (3$\farcs$5, corresponding to 0.8 pc) is indicated in the lower left of each plot. In addition to matching the spatial resolution of the cubes, we have added Gaussian noise at the beam scale to reach a uniform 1$\sigma$ noise of 0.25 K for 12CO and 0.15 K for 13CO. These values were scaled down by $\sqrt{2.5}$ for 30 Dor, for which the channel map spacing is 0.5 km s*-1* rather than 0.2 km s*-1*. Since line widths in the 30 Dor region are substantially larger ($\gtrsim$1 km s*-1*) than in other parts of the LMC (Indebetouw et al. 2013; Paper I), the mismatch in channel width and noise should not have a significant impact on our results.

The ALMA maps reveal a wealth of structure, which we plan to investigate more deeply in future papers, but what is immediately apparent is that the emission structures are nested rather than discrete, and that much of the brightest emission occurs in filaments. In these and subsequent plots, we will present the results for the six clouds in the following order: 30 Dor, N59C, A439, GMC104, GMC1, and PCC. This ordering is by decreasing median 8$\mu$m intensity (Figure 1, left), and reflects a sequence of decreasing star formation activity.

3.2 LTE Analysis

Following Paper I, we conduct a simple local thermodynamic equilibrium (LTE) analysis to infer the 13CO column density from the observed 12CO and 13CO emission. The analysis assumes that both lines share a common excitation temperature (e.g., Nishimura et al., 2015). We assume the 12CO(2–1) line is optically thick at line center and not subject to beam dilution, so that for a given line of sight the excitation temperature is uniform and given by (e.g., Bourke et al., 1997):

[TABLE]

where the beam filling fraction $f_{\rm bm}$ is assumed to be 1, $T_{\rm 12,pk}$ is the peak temperature of the 12CO line profile, and

[TABLE]

The beam-averaged 13CO optical depth is then calculated from the brightness temperature $T_{13}$ at each position and velocity in the cube by solving

[TABLE]

again assuming $f_{\rm bm}=1$. Since $\tau_{13}$ varies linearly with $T_{13}$ in the optically thin limit (for a given value of $T_{\rm ex}$), we allow negative values of $\tau_{13}$ due to noise. Given our assumption of a single $T_{\rm ex}$ at each sky position, and because of the limited range of $T_{\rm ex}$ for most clouds (as discussed below), these noise values tend to average out when integrated in the cube.

For any given $T_{\rm ex}$, there is a maximum allowed value of $T_{13}$ beyond which $\tau_{13}$ becomes undefined. For $T_{\rm ex}=6$ K we require $T_{13}<2.8$ K for $J=1\rightarrow 0$ and $T_{13}<2$ K for $J=2\rightarrow 1$. To reduce the number of undefined values we impose a minimum value on $T_{\rm ex}$ of $T_{\rm floor}$ = 6 K under the assumption that lower inferred values reflect beam dilution of 12CO (i.e., $f_{\rm bm}<1$). We chose $T_{\rm floor}$ based on an examination of the distribution of inferred $T_{\rm ex}$ for pixels detected at high signal-to-noise ratio (SNR$>$5) in 13CO. Pixels with $T_{\rm ex}<6$ K are extremely rare and account for 0.2% or less of the high SNR pixels in any given cloud. As a result, imposition of the floor has little effect on the high SNR column densities, while avoiding problems with artificially high column densities due to low nominal $T_{\rm ex}$ values in low SNR regions.

Following determination of $T_{\rm ex}$ and $\tau_{13}$, the total 13CO column density in cm*-2*, summed over all rotational levels, is calculating using (Garden et al., 1991):

[TABLE]

where $J$ is the rotational quantum number of the lower state, $B$ is the rotational constant for 13CO (55.1 GHz), and $\mu$ is the dipole moment of 13CO (0.112 debye).

The calculated uncertainties in $N(^{13}{\rm CO})$ mainly reflect the uncertainties due to map noise. The assumption of a single $T_{\rm ex}$ that can be derived from $T_{\rm 12,pk}$, as well as the definition of the 3D masks used for integration, introduce additional systematic uncertainties. The most important effect is beam dilution, which reduces $T_{\rm 12,pk}$ when smoothing to a common resolution of 3$\farcs$5; our resulting underestimate of $T_{\rm ex}$ would lead us to overestimate $\tau_{13}$ (Equation 3), although this is mitigated somewhat by the fact that we would have also overestimated $f_{\rm bm}$. The net effect on $N$(13CO) of lowering $\tau_{13}$ while raising $T_{\rm ex}$ is difficult to assess, and would be better constrained by measuring the 13CO excitation with additional $J$ transitions.

We derive an LTE-based estimate of clump masses and column densities by scaling $N$(13CO) to $N$(H2) using a fixed abundance ratio of

[TABLE]

Although the 13CO abundance in the LMC is not well-constrained by observations, and may be subject to spatial variations, our adopted value is consistent with values inferred or adopted in previous work (Heikkilä et al., 1999; Mizuno et al., 2010; Fujii et al., 2014). Henceforth we denote LTE-based mass surface densities as $\Sigma_{\rm LTE}$.

3.3 Structural Decomposition

To identify significant emission structures in the datacubes we used the Python package astrodendro111http://www.dendrograms.org, which decomposes emission into a hierarchy of structures (Rosolowsky et al., 2008; Shetty et al., 2012; Colombo et al., 2015). Our procedure follows closely that used in our previous analysis of the PCC and 30 Dor clouds (Paper I). The algorithm identifies local maxima in the cube above the 3$\sigma_{\rm rms}$ level that are also at least 2.5$\sigma_{\rm rms}$ above the merge level with adjacent structures. Each local maximum is required to span at least two synthesized beams in area. Isosurfaces surrounding the local maxima are categorized as trunks, branches, or leaves according to whether they are the largest contiguous structures (trunks), are intermediate in scale (branches), or have no resolved substructure (leaves). For a given cloud, the trunks do not overlap other trunks and leaves do not overlap other leaves, but every trunk can be decomposed into leaves (and usually branches). The resulting dendrogram for CO in the A439 cloud is shown in Figure 5, with leaves, branches, and trunks colored green, black, and maroon, respectively.

The basic properties of the identified structures are also determined by astrodendro, including their spatial and velocity centroids ($\bar{x},\bar{y},\bar{v}$), the integrated flux $S$, rms line width $\sigma_{v}$ (defined as the intensity-weighted second moment of the structure along the velocity axis), the position angle of the major axis (as determined by principal component analysis) $\phi$, and the rms sizes along the major and minor axes, $\sigma_{\rm maj}$ and $\sigma_{\rm min}$. All properties are determined using the “bijection” approach discussed by Rosolowsky et al. (2008), which associates all emission bounded by an isosurface with the identified structure. From these basic properties we have calculated additional properties, including the effective rms spatial size, $\sigma_{r}=\sqrt{\sigma_{\rm maj}\sigma_{\rm min}}$, the spherical radius $R=1.91\sigma_{r}$, following S87, the luminosity $L=Sd^{2}$, adopting $d=50$ kpc (Pietrzyński et al., 2019), the virial mass $M_{\rm vir}=5\sigma_{v}^{2}R/G$, derived from solving the equilibrium condition

[TABLE]

and the luminosity-based mass (from 12CO)

[TABLE]

where $X_{2}=1$ for a standard (Galactic) CO to H2 conversion factor (Bolatto et al., 2013). In this paper we have adopted $X_{2}=2.4$ for the CO(1–0) line based on the virial analysis of the MAGMA GMC catalog by Hughes et al. (2010), and $X_{2}=3$ for the CO(2–1) line assuming a CO(2–1)/CO(1–0) line ratio of 0.8. The line ratio is known to vary with cloud conditions (Sorai et al. 2001; Paper I), with values $\sim$0.5 for clouds in the outskirts of the LMC and rising to $\sim$1.2 near 30 Dor, so our adoption of a constant value is only approximate.

Our expression for virial mass is not fully self-consistent given that S87 assumed a truncated $\rho\propto r^{-1}$ density profile in deriving $R=1.91\sigma_{r}$, whereas we have assumed a constant density sphere in deriving the potential energy ${\cal W}$. For $\sigma_{v}$ in km s*-1* and $\sigma_{r}$ in pc, our expression for $M_{\rm vir}$ simplifies to $2.22\times 10^{3}\,\sigma_{v}^{2}\sigma_{r}$ $M_{\odot}$, which is about 10% higher than the equivalent expression in S87. For simplicity, and consistent with the crudeness of our treatment of the virial theorem, we choose not to correct for this offset.

As with any emission segmentation approach, one can justifiably question the reality of the recovered structures in the dendrogram (see discussion in Shetty et al., 2010; Beaumont et al., 2013). Since we are principally interested in recovering the line widths of structures spanning a range of different sizes, our results should be consistent with other approaches that measure emission properties within bounded regions of the data cube—i.e., that avoid extrapolation of properties beyond the structure boundaries. Any approach that seeks to recover line widths at different size scales must ultimately compare nested rather than disjoint structures, and so is unlikely to yield dramatically different results. Note that we do not attempt to generate histograms (e.g., clump mass spectra) of the structure properties. Even if limited to a set of disjoint structures (e.g., the subset of dendrogram leaves), such histograms would be difficult to interpret given the continuous nature of the CO emission (Figure 3 and 4) and the fact that leaf structures will tend to be similar by construction. Uniformly sampled property measurements (e.g., Hughes et al., 2013; Leroy et al., 2016) are probably better suited for obtaining parameter distributions.

3.4 Size vs. Line Width Relations

Figures 6 and 7 show the individual $R$–$\sigma_{v}$ relations for the six clouds in 12CO and 13CO respectively. Each relation is fit by a power-law model of the form

[TABLE]

with the slope ($a_{1}$) and intercept ($a_{0}$) of the fitted line reported in Table 2. Following Paper I, fitting is performed using the kmpfit module of the Python package Kapteyn (Terlouw & Vogelaar, 2015), which treats errors in both axes using the effective variance method. Points in the gray shaded regions are excluded from fitting due to resolution limitations. However, because a sharp truncation of the data can skew the fit results, we have also repeated the fitting without excluding the shaded regions. The resulting slopes and intercepts are usually within the quoted 1$\sigma$ uncertainties, but are sometimes discrepant by nearly 3$\sigma$. We have therefore chosen to plot a (conservative) 3$\sigma$ confidence interval as blue shading in Figures 6 and 7, while still listing the 1$\sigma$ uncertainties in Table 2. We note that the fit parameters are quite sensitive to the selection of data points to be fitted and their relative uncertainties, and thus fits to particular subsets of the data, or to differently weighted data, could differ from the reported fits by much more than 1$\sigma$. Moreover, a power law is generally a poor fit to the data, as reflected by the reduced $\chi^{2}$ values all being substantially greater than 1 (Table 2).

In Figures 6 and 7 the fiducial relation of S87 for Galactic clouds is shown in pink. The LMC clouds appear roughly consistent with this relation, although smaller structures tend to exhibit a wide range of line widths. Comparing the 12CO and 13CO relations for a given cloud, we find no systematic differences between the fitted slopes or intercepts. Particularly for the 13CO line, however, there appears to be a progression from large to small line widths at a given size when comparing clouds with decreasing star formation activity. This is consistent with the discrepancy in the $R$–$\sigma_{v}$ relations for 30 Dor and PCC previously noted in Paper I, but the four additional clouds clearly demonstrate a continuous variation across the sample. Variations in the fitted slope, in particular a fairly steep slope in the case of N59C, make this trend less obvious when examining only the fitted values of $a_{0}$; the trend is revealed more clearly when referencing the $R$–$\sigma_{v}$ relation to a common slope (§3.6).

Figure 8 provides a summary view of the $R$–$\sigma_{v}$ relation across the six clouds, with points for each cloud colored according to the median 8$\mu$m intensity within the cloud. The aggregate relation for 12CO can be described by

[TABLE]

with an rms scatter of $\varepsilon=0.25$ dex in $\sigma_{v}$ (see Table 2). In these and subsequent plots, the fitting is still performed on individual points weighted by their uncertainties, with binned values plotted to illustrate overall trends. The aggregate scatter is substantially larger than the scatter for any single cloud, indicating that differences between clouds are substantial. Again, these differences are in the sense of higher line widths for 8$\mu$m-bright clouds (blue and magenta colors) and smaller line widths for 8$\mu$m-faint clouds (orange and red colors). This is seen for the 12CO-identified structures (upper left panel), the 13CO-identified structures (upper right panel), and persists also when the marker colors are based on the local rather than cloud-scale 8$\mu$m intensity (lower left panel). These first three panels of Figure 8 provide direct evidence for a coupling between the IR brightness and the CO velocity dispersion.

3.5 IR Brightness or Cloud Surface Density?

We have seen that the six clouds in our sample all show a consistent trend of increased line width at a given size with increased IR surface brightness. A natural interpretation of this trend is that IR brightness tracks star formation activity and that higher star formation activity results in stronger stirring of the ISM due to feedback. Here we consider an alternative interpretation, which posits that molecular clouds and their substructures lie close to simple virial equilibrium as defined by Equation (6). Increased line width then reflects higher mass surface densities, which also tend to correlate with higher star formation activity because the molecular gas is more susceptible to star formation or simply more abundant. Figure 8 (bottom right) shows that when structures are color coded by the local surface density derived from LTE analysis, there is shift in surface density from the bottom to top of the plot which closely resembles the shift in IR brightness seen in the other panels.

The normalization of the $R$–$\sigma_{v}$ relation, $v_{0}=\sigma_{v}/R^{1/2}$, is related to a virial mass surface density,

[TABLE]

so the plot of virial vs. luminous surface density (the so-called boundedness plot, Figure 9) provides another view of the trends seen in the size-line width relation. As expected, the 8$\mu$m-bright clouds are seen at higher $\Sigma_{\rm vir}$ in Figure 9, i.e. toward the top of the figure. It is clear, however, that $\Sigma_{\rm vir}$ (and thus $v_{0}$) is also correlated with observational measures of surface density, such as $\Sigma_{\rm CO}$ and $\Sigma_{\rm LTE}$; indeed, most of the 13CO structures lie close to the line of simple virial equilibrium ($\Sigma_{\rm vir}=\Sigma_{\rm obs}$, the diagonal dashed lines in each panel of Figure 9). We have also drawn lines of constant external pressure $P_{\rm ext}$, derived from setting the right-hand side of Equation (6) equal to $4\pi R^{3}P_{\rm ext}$ (Field et al., 2011):

[TABLE]

We note that the observed 12CO and 13CO structures show considerable scatter around the virial equilibrium line, with no clear “threshold” density above which structures tend to be virialized. Furthermore, they do not appear consistent with confinement by a single value of $P_{\rm ext}$, although a role for a variable external pressure cannot be excluded.

If the IR brightness is largely responding to $\Sigma$, as the virial interpretation would imply, we should find a good local correlation between 8$\mu$m intensity and $\Sigma$. To investigate this we compare in Figure 10 the 8$\mu$m intensity averaged within each 13CO dendrogram structure with the surface density of the structure measured in three different ways: from the virial theorem (left panel), CO brightness (middle panel), or LTE analysis (right panel). While a correlation is apparent in all three panels, it becomes progessively weaker for mass estimators that are less sensitive to $\sigma_{v}$ ($\Sigma_{\rm vir}\propto\sigma_{v}^{2}$, and $\Sigma_{\rm CO}\propto\sigma_{v}$ for an optically thick CO line with a constant saturation brightness temperature.) This suggests that the correlation between 8$\mu$m intensity and $\Sigma$, while significant, is driven in part by sensitivity of common mass estimators to the CO line width.

On the other hand, 8$\mu$m intensity alone does not appear to fully account for changes in $v_{0}$. This is apparent from Figure 11, which shows that at fixed 8$\mu$m intensity, $\Sigma_{\rm vir}\propto v_{0}^{2}$ (indicated by the marker colors) increases with independent measures of $\Sigma$ (in this case, the LTE-based value, but the CO-based surface densities show the same pattern). Conversely, at fixed $\Sigma$, $v_{0}$ increases with 8$\mu$m intensity. Although our crude LTE treatment of CO excitation may introduce intrinsic correlations between the two axes (e.g., by underestimating the column density in cold regions), we interpret this result as pointing to significant roles for both feedback (represented by the abscissa) and self-gravity (represented by the ordinate) in setting the normalization of the $R$–$\sigma_{v}$ relation (represented by the the marker color). With a slope significantly below unity, the data are also consistent with a steep dependence of the SFR surface density on $\Sigma$, a point we return to in §4.

3.6 Cloud-averaged correlations

We can also examine the relationship between $v_{0}$, IR emission, and mass surface density on cloud-averaged scales, as seen in Figure 12. Here the error bars are conservatively drawn to indicate the standard deviation of the individual values of $v_{0}$ across the set of dendrogram structures, or across the cloud in the case of the abscissa. We see again that the cloud-averaged $v_{0}$ (shown as the vertical axis on each panel) correlates well with the median 8$\mu$m intensity; the upper left panel of Figure 12 shows this is consistent across all clouds and both line tracers (12CO and 13CO). The observed trend is driven in part by the much higher $v_{0}$ values in 30 Dor, but including values for the N113 cloud, based on a preliminary analysis of ALMA 12CO (2–1) and 13CO (2–1) maps (Sewiło et al., in preparation), further strengthens the case for a general trend across clouds.

Considering that 8$\mu$m is a complex tracer (attributed to FUV-irradiated PAH molecules) that could be intrinsically related to molecular gas column density, we can check whether this correlation persists with other tracers of star formation activity. As the two panels on the right side of Figure 12 show, a comparably good correlation is seen with 24$\mu$m intensity and dust temperature, both of which should be responsive to the FUV radiation field in ways that are independent of the PAH emission. At the same time, $v_{0}$ also shows a good correlation with dust column density inferred from FIR emission (bottom left panel). This provides further evidence that both star formation activity and mass surface density contribute to setting the characteristic $v_{0}$ for a cloud.

4 Discussion and Conclusions

The main results of our study can be summarized as follows. Firstly, while there is substantial scatter in the line width at a given size within a cloud, the scatter among clouds is even larger, allowing one to meaningfully assign a mean value of the normalization $v_{0}$ to each cloud in our sample. Secondly, the characteristic $v_{0}$ for a cloud correlates independently with both its mean IR brightness and its mean surface density—although we emphasize that the surface density is inferred from CO or 13CO line intensity, which may be sensitive to variations in radiative excitation and gas velocity dispersion. Our results are consistent with at least two interpretations. The first is that higher line widths are due to energetic feedback from recent star formation, which tends to occur in high density regions. The second is that the higher line widths are due to gravity-driven motions in high column density structures, with some of these structures collapsing to form stars.

The ambiguity between these two interpretations reflects the ongoing debate over which of the possible sources of turbulence in the interstellar medium dominates: gravitational collapse, stellar feedback, or large-scale galactic dynamics (e.g., Padoan et al., 2016; Ballesteros-Paredes et al., 2011; Krumholz & Burkhart, 2016). It is reasonable to expect a combination of these factors to be at work, with the balance between them shifting as a function of spatial scale. Our ALMA observations probe scales of 1–50 pc, substantially smaller than the diameters of ionized superbubbles (Chu et al., 1995) or the inferred disk thickness (Elmegreen et al., 2001). Thus we lack sensitivity to energy injection on larger scales, although large-scale flows can still generate a turbulent cascade to the smaller scales that we probe (Klessen & Hennebelle, 2010). A model by Krumholz & Burkhart (2016) offers a simple prediction for distinguishing feedback from gravitationally dominated turbulence on kiloparsec scales. Relying on a state of vertical dynamical equilibrium and marginal Toomre stability, they predict ${\rm SFR}\propto\sigma_{v}^{2}$ for a feedback-driven model and ${\rm SFR}\propto f_{g}^{2}\sigma_{v}$ for a gravity driven model, where $f_{g}$ is the gas fraction.

However, these scalings may not apply on the scales of individual molecular clouds, where equilibrium may not hold and the gas density can deviate significantly from the critical Toomre value. On the scale of individual molecular clouds, feedback simulations predict that the integrated star formation efficiency increases with $\Sigma$, since more feedback is required to halt collapse, and the ratio of the self-gravity force to the rate at which feedback injects momentum depends on surface density (e.g., Fall et al., 2010; Raskutti et al., 2016; Grudić et al., 2018). This is qualitatively consistent with our finding of a steep dependence of star formation rate on $\Sigma$ (Figure 11), although a quantitative comparison will require taking a detailed census of recent star formation in these clouds.

In the absence of a simple analytic prediction to distinguish feedback-driven from gravitationally-driven turbulence, we can examine whether the highest dispersion molecular gas is associated with small or with large structures. We noted previously in Paper I that in the PCC, and to some extent in the 30 Dor cloud, some of the highest dispersions are found in the smallest resolved structures (the dendrogram leaves, green symbols in Figure 6). This is most easily apparent in the $R$–$\sigma_{v}$ diagram, or equivalently a plot of specific kinetic energy ($\sigma_{v}^{2}$) as a function of size scale. We present such a plot in Figure 13, now including all six clouds. We note that the highest specific K.E. is generally associated with the largest structures, especially when using the 13CO line, which is less sensitive to opacity broadening (see below). The most quiescent clouds, PCC and GMC1, show a somewhat flatter upper envelope in the distribution of points, indicating significant energy injection on small scales. The other clouds, especially GMC104 and A439, seem more consistent with a simple energy cascade from large to small scales. It therefore appears difficult to ascribe a single dominant mechanism for turbulent energy injection across the entire sample.

We note that interpretations of $v_{0}$ in terms of feedback or gravity rely on our ability to infer physical quantities from measured observables. In doing so we should be aware of the following caveats.

Opacity broadening: Hacar et al. (2016) have analyzed the overestimate of the intrinsic velocity dispersion that results from high line opacity. The overestimate amounts to up to a factor of 2–3 in $\sigma_{v}$ (and thus up to an order of magnitude in $\Sigma_{\rm vir}$) in the case of the 12CO line, and may be responsible for many of the points which scatter above the virial equilibrium line in Figure 9 (left). We note, however, that we obtain consistent results for $v_{0}$ from both the 12CO and 13CO analyses, while the latter should be less susceptible to (though not immune from) opacity effects. 2. 2.

Mass uncertainties: We have estimated structure masses using either a constant $X$-factor ($\Sigma_{\rm CO}$) or a simple LTE analysis ($\Sigma_{\rm LTE}$), both of which involve simplifying assumptions and require adopting highly uncertain abundance and line ratios. In particular, $\Sigma_{\rm LTE}$ is based on an assumed $T_{\rm ex}$ that is likely underestimated due to beam dilution (§3.2). Although a systematic error in mass would not affect our characterization of the size-line width relation, it would affect our ability to disentangle the influences of line width and surface density on 8$\mu$m emission.

Although we are unable to identify the dominant driving mechanism for turbulence in our clouds, there are some obvious next steps that will bring us closer to doing so. Direct comparison with numerical simulations, translated into the observational domain using radiative transfer modeling, will help interpret the scatter in line width, particularly on small scales. Ongoing wide-field mapping of the LMC with the Morita 7-m array of ALMA will extend our analysis to larger scales and thus better constrain the form of the $R$–$\sigma_{v}$ correlation. In addition, improved characterization of the dust mass, ionized gas properties, and the young stellar population, made possible with near-infrared and optical surveys and (in the near future) with the James Webb Space Telescope (JWST), will constrain more tightly the available energy from stellar feedback and provide mass estimates that are independent of CO emission. These advances should enable substantial progress is characterizing the energy flow within turbulent molecular clouds.

Images and analyses presented in this paper are available for download at the URL https://mmwave.astro.illinois.edu/almalmc/, or an updated URL can be found at the arxiv.org listing for this paper.

We thank Mike Grudić, Michael Fall, and the anonymous referee for helpful comments on an earlier draft of this paper. We thank Dyas Utomo for assistance with the dust maps. The Mopra radio telescope is part of the Australia Telescope National Facility which is funded by the Australian Government for operation as a National Facility managed by CSIRO. This paper makes use of the following ALMA data: ADS/JAO.ALMA #2011.0.00471.S, ADS/JAO.ALMA #2013.1.00832.S, and ADS/JAO.ALMA #2016.1.00193.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. R.I. acknowledges support from NSF AAG award 1313276. A.H. acknowledges support from the Centre National d’Etudes Spatiales (CNES). K.T. was supported by NAOJ ALMA Scientific Research Grant Number 2016-03B. The work of M. S. was supported by NASA under award number 80GSFC17M0002. Part of this research was conducted at the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. R.S.K. and S.C.O.G. acknowledge financial support from the German Science Foundation (DFG) in the Collaborative Research Center SFB 881 The Milky Way System (subprojects B1, B2, and B8). This research made use of astrodendro, a Python package to compute dendrograms of astronomical data, SCIMES, a Python package to find relevant structures in dendrograms of molecular gas emission using the spectral clustering approach, and Astropy, a community-developed core Python package for astronomy. Mopra, ALMA, Herschel, Spitzer