ABSTRACT
We study the probability distribution function (PDF) of mass surface densities, Σ, of infrared dark cloud (IRDC) G028.37+00.07 and its surrounding giant molecular cloud. This PDF constrains the physical processes, such as turbulence, magnetic fields, and self-gravity, that are expected to be controlling cloud structure and star formation activity. The chosen IRDC is of particular interest since it has almost 100,000 solar masses within a radius of 8 pc, making it one of the most massive, dense molecular structures known and is thus a potential site for the formation of a "super star cluster." We study Σ in two ways. First, we use a combination of NIR and MIR extinction maps that are able to probe the bulk of the cloud structure up to Σ ∼ 1 g cm−2(AV ≃ 200 mag). Second, we study the FIR and submillimeter dust continuum emission from the cloud utilizing Herschel-PACS and SPIRE images and paying careful attention to the effects of foreground and background contamination. We find that the PDFs from both methods, applied over a ∼20'(30 pc)-sized region that contains ≃1.5 × 105 M⊙ and enclosing a minimum closed contour with Σ ≃ 0.013 g cm−2 (AV ≃ 3 mag), shows a lognormal shape with the peak measured at Σ ≃ 0.021 g cm−2 (AV ≃ 4.7 mag). There is tentative evidence for the presence of a high-Σ power-law tail that contains from ∼3% to 8% of the mass of the cloud material. We discuss the implications of these results for the physical processes occurring in this cloud.
Export citation and abstract BibTeX RIS
1. INTRODUCTION
The probability distribution function (PDF) of mass surface density, Σ, is one of the simplest metrics of interstellar cloud structure. This Σ-PDF is, in principle, much easier to observe than other distributions, such as volume density, thus making it a convenient metric with which to compare observed and simulated clouds. The Σ-PDF shape should be sensitive to physical processes occurring in the clouds. For example, simulations of driven supersonic hydrodynamic (and if including magnetic fields, super-Alfvénic) turbulence of non-self-gravitating gas in periodic boxes yield lognormal Σ-PDFs (e.g., Federrath 2013; Padoan et al. 2014), i.e., the area-weighted PDF, pA, can be well fit by a lognormal:
where 
is mean-normalized Σ. The lognormal width, σlnΣ, grows as turbulent Mach number increases. In simulations with self-gravity, Σ-PDFs are seen to develop high-end power-law tails, perhaps tracing regions undergoing free-fall collapse (Cho & Kim 2011; Collins et al. 2011; Kritsuk et al. 2011; Federrath & Klessen 2013). However, simulations of self-gravitating, strongly magnetized (trans-Alfvénic), turbulent clouds with non-periodic boundary conditions are also needed for comparison with observed Σ-PDFs. Such clouds are expected to have smaller star formation efficiencies per mean free-fall time, 
ff, likely implying they would have smaller mass fractions in any high-Σ power-law tail. Accurate quantification of the Σ-PDF in real star-forming clouds is needed to constrain theoretical models.
Observationally, Kainulainen et al. (2009, hereafter K09) performed NIR extinction mapping of ∼20 nearby clouds, ranging from "quiescent," non-star-forming clouds to more active clouds. They found quiescent cloud Σ-PDFs are well described by lognormals, while star-forming clouds have high-Σ power-law tails. Note that in practice a lognormal is fit to the observed 
, which is then used to derive 
(over the considered range of Σ), which then defines the mean 
, where σlnΣ is standard deviation of lnΣ (Butler et al. 2014, hereafter BTK14).
However, the ability of this observational method to accurately measure the position of the Σ-PDF peak, typically at Σ ∼ 0.01 g cm−2 (i.e., AV ∼ 2 mag; we adopt conversion Σ/(1 g cm−2) ≡ Av/(224.8 mag); Kainulainen & Tan 2013, hereafter KT13) in the K09 clouds, has been questioned by Schneider et al. (2015b) and Lombardi et al. (2015) due to difficulties of disentangling foreground and background contributions. Lombardi et al. argued that low-Σ PDF uncertainties are so large that observed PDFs are all consistent with power laws.
Extending these studies to denser, higher-Σ clouds, perhaps more typical of most Galactic star formation, KT13 used combined NIR (Kainulainen et al. 2011) + MIR (Butler & Tan 2012, hereafter BT12) extinction maps to study Σ-PDFs of 10 infrared dark clouds (IRDCs). These dense structures, typically several kiloparsecs away, have Σ-PDFs extending to ∼1 g cm−2 (e.g., Tan et al. 2014). One uncertainty of these maps is choice of opacity (i.e., at 8 μm) per unit Σ, with variations ∼30% expected for different dust models, e.g., moderately coagulated thin and thick ice mantle models of Ossenkopf & Henning (1994, hereafter OH94). KT13 considered simple, rectangular regions enclosing contours of Σ = 0.03 g cm−2 (AV ≃ 7 mag), i.e., "complete" for Σ above this value. However, two caveats limit this completeness: first, MIR-bright regions are not treated by extinction mapping, so are excluded from the PDF; second, extinction mapping saturates at high-Σs, depending on MIR image depth, typically at ∼0.3–0.5 g cm−2 for Spitzer GLIMPSE (Churchwell et al. 2009) images (BT12). At low-Σs, better probed by NIR extinction mapping, systematic uncertainties are ∼0.01 g cm−2. The KT13 Σ-PDFs, extending down to AV = 7 mag, did not cover PDF peaks well, so provided weak constraints on true shapes: PDFs could be fit with lognormals or power laws.
BTK14 presented a higher dynamic range Σ-map of IRDC G028.37+00.07, both to higher Σ ≃ 0.7 g cm−2 and lower Σ ≃ 0.013 g cm−2 (AV ≃ 3 mag). From a 20' × 19' region around the IRDC (Figure 1(a)), BTK14 derived a Σ-PDF where the peak was observed at Σpeak ≃ 0.03 g cm−2. Furthermore, the entire PDF from 0.013 to 0.4 g cm−2 was well described by a lognormal with 
(
) and σlnΣ' = 1.4. There appears to be relatively little mass in a high-Σ power-law tail, surprising given the IRDC (KT13) and GMC (Hernandez & Tan 2015, hereafter HT15) are both self-gravitating with virial parameters close to unity and some star formation has already started.
Figure 1. (a) Top: 160 μm image of IRDC G028.37+00.07 (ellipse) and surroundings. Smaller rectangle shows 20' × 19' area analyzed by BTK14 and in this Letter. Larger rectangle indicates approximate area, ∼1° × 0
7, studied by S15b. (b) Bottom left: Galactic background model for 160 μm intensity for the BTK14 region evaluated with the Small Median Filter (SMF) method. (c) Bottom right: same as (b), but with the Galactic Gaussian (GG) method. The seven vertical lines are locations of profiles shown in Figure 2.
Download figure:
Standard image High-resolution imageAnother method to measure Σ is via submillimeter dust continuum emission. However, this also depends on dust temperature, T, so multiwavelength studies are needed to probe the spectral energy distribution (SED) peak. The best data for this come from Herschel-PACS and SPIRE observations probing 70–500 μm. However, derived maps have angular resolution ∼20'', i.e., ∼10× worse than the NIR+MIR extinction maps.
Schneider et al. (2015b, hereafter S15b) utilized Herschel-derived Σ-maps to study the same IRDC/GMC examined by BTK14. They derived Σ-PDFs for the IRDC ellipse region (Figure 1(a)) and a surrounding "GMC" region defined by a 13CO(1-0) emission contour equivalent to AV ∼ 2 mag, extending approximately over the larger rectangle shown in Figure 1. Note that this is significantly larger than the BTK14 region. S15b found their IRDC Σ-PDF was well fit by a single power law for AV ≳ 30 mag. The GMC region could also be fit with a power law, especially for AV ≳ 40 mag. Below this, S15b claimed to detect a peak in the Σ-PDF at AV ∼ 20 mag. S15b proposed that the presence of power-law tails indicated the cloud was undergoing multi-scale, including global, quasi-free-fall collapse.
Here, we revisit the Σ-PDF toward IRDC/GMC G028.37+00.07, especially comparing PDFs derived from dust extinction and emission methods. We examine reasons for the different results of BTK14 and S15b and derive the most accurate Σ-PDF of this massive protocluster.
2. METHODS
2.1. NIR+MIR Extinction Derived Σ Map
We utilize the NIR+MIR extinction map from BTK14 with one modification. The NIR extinction map, based on statistical estimates of stellar extinctions in ∼30' regions, requires choosing an "off-position" of negligible local IRDC/GMC extinction. BTK14 utilized an off-position at l = 284, b = −004. However, S15b noted this location may be too close to the GMC, which is confirmed in the 13CO(1-0) map of HT15. We therefore choose a new off-position at l = 283, b = +03 with AV that is 3 mag smaller than the previous position. The net effect is to add an offset of AV = 3 mag to the BTK14 map. We will see this has only a minor effect on the Σ-PDF: in particular, the mean extinction remains close to 
.
Conversion of the extinction map into a Σ-map requires an assumption about opacity per unit mass at a given wavelength. Here, for the Spitzer-IRAC ∼ 8 μm map we adopt 7.5 cm2 g−1 (BT12). The NIR+MIR combination assumes τ8 μm = 0.29τK (KT13). For conversion to AV we follow KT13, adopting Σ/(1 g cm−2) ≡ AV/(224.8 mag). Overall, from consideration of different dust models (BT12), we estimate ∼30% systematic uncertainties due to these opacity choices.
2.2. Submillimeter Emission Derived Σ Map
We use PACS and SPIRE images from the Herschel Infrared GALactic plane survey (Hi-GAL; Molinari et al. 2010). Derivation of Σ and T is performed via pixel-by-pixel graybody fits to the 160, 250, 350, and 500 μm data (e.g., Battersby et al. 2011; S15b), first re-gridded at the 500 μm image resolution of 36''. Specifically, Σ and T are derived via
where Iν is the observed intensity of the corresponding band, Bν is the filter-weighted value of the Planck function, τν is the optical depth, and κν is the filter-weighted opacity.
However, to derive local cloud properties, the images need correction for foreground (FG) and background (BG) diffuse ISM emission along the line of sight. Assuming emission from the cloud in these submillimeter bands is optically thin, one method is to estimate FG+BG emission as one combined, constant intensity. S15b adopted this method, selecting a region outside the GMC of interest for FG+BG column density corresponding to AV ≃ 2 mag, which was then subtracted as a constant offset.
We first assess the FG based on the ISM model derived from observed 24 μm intensities toward "saturated" dark regions of the IRDC by Lim & Tan (2014, hereafter LT14) and Lim et al. (2015, hereafter LCT15). These saturated regions, also seen at shorter wavelengths and possibly at 70 μm, are where observed intensities are similar to within instrumental noise in spatially independent locations. This is assumed to be caused by the IRDC blocking essentially all BG light, so the observed intensity is that of a spatially smooth FG. Several 24 μm FG measurements across the IRDC are made and then averaged to estimate a mean. Then the Draine & Li (2007) diffuse ISM SED model is normalized to this value to predict 160, 250, 350, and 500 μm filter-weighted FG values. These are subtracted from the submillimeter images, i.e., the FG is first assumed to be constant across the IRDC/GMC.
Once we have FG-subtracted images, we next assess the BG in two ways. First, we assess BG intensity at each wavelength adapting the "Small Median Filter" (SMF) method (Butler & Tan 2009, hereafter BT09), which is applied outside the IRDC ellipse with square filter of 35' size. BG emission behind the IRDC ellipse is estimated via interpolation from the surrounding regions (Figure 1(b)), following the BT09 weighting scheme.
As a second, "Galactic Gaussian" (GG), method, we follow Battersby et al. (2011) and assume that the Galactic BG follows Gaussian profiles in latitude. We fit Gaussians to the minimum intensities exhibited along strips with longitude width ∼2' (see Figures 1(c) and 2 for examples at 160 μm). The method involves clipping higher intensity values arising from discrete clouds and iteratively converges on a final result. Note that here we also assume the FG intensity is a Gaussian of the same width and subtract that off the images before deriving the final BG model.
Figure 2. Seven example profiles of the Galactic Gaussian (GG) background model (dashed lines), evaluated along l values shown in Figure 1(c). Solid lines show the total observed intensity (after FG subtraction), dashed lines show the fitted Gaussian function to this profile, and the dotted line is the GG FG model. Vertical dot-dashed lines show the Galactic latitude range of the BTK14 region.
Download figure:
Standard image High-resolution imageThe above methods result in FG+BG-subtracted images at 160, 250, 350, and 500 μm. Then, at each pixel we fit the graybody function (Equation (2)), including its filter response weighting, to derive Σ and T. This fitting requires an assumed form of κν. This has been studied via MIR to FIR extinction by LT14 and LCT15, who find evidence of generally flatter extinction laws over the ∼8–70 μm range as Σ increases, consistent with OH94 and Ormel et al. (2011) dust models that include grain growth via ice mantle formation and coagulation. For consistency with these extinction results, we adopt the OH94 thin ice mantle model with 105 years of coagulation at a density of nH = 106 cm−3 as our fiducial. At submillimeter wavelengths, this model exhibits κν = 0.1(ν/1000 GHz)β cm2 g−1 with β ≃ 1.8. We will also explore the effects of varying β from 1.5 to 2.
Finally, we generate higher-resolution (HiRes) Σ maps by re-gridding to the 250 μm image pixels (18'' resolution; 6'' pixels) and then repeating the above analysis, but now fixing temperatures from the lower angular resolution maps. These HiRes Σ maps are better able to probe smaller, higher Σ structures.
3. RESULTS
3.1. Mass Surface Density and Temperature Maps
Figure 3 shows the submillimeter emission derived Σ (HiRes) and T maps of the IRDC, starting from maps derived with no FG and BG subtraction, and then showing the effects of the two background estimation methods (SMF and GG). The overall result of FG+BG subtraction on the Σ map is to effectively remove ∼0.1 g cm−2 across the cloud. As we will see this has a major effect on the shape of the Σ-PDF. Note that for study of the BTK14 region, we consider the GG method to be superior to SMF as it allows for the large-scale structure of the Galactic plane. Also, the SMF method tends to underestimate mass surface densities in the structures just outside the IRDC ellipse.
Figure 3. Submillimeter emission derived HiRes Σ (top row) and T (middle row) maps based on no FG+BG subtraction (left), SMF (middle), and GG (right) methods. Bottom row shows images from Spitzer-IRAC 8 μm, -MIPS 24 μm, and Herschel-PACS 70 μm.
Download figure:
Standard image High-resolution imageThe temperature maps are also strongly affected by FG+BG subtraction, which leads to a lowering of the temperatures measured in the IRDC, as well as revealing warmer localized patches in the surroundings. Some of these appear to correspond to 8, 24, and 70 μm brighter regions (Figure 3, lower panels).
The top row of Figure 4 compares the FG+BG-subtracted Σ maps derived from NIR+MIR extinction and from submillimeter emission. White patches in the extinction map are locations of bright MIR emission that prevent an absorption measurement against the Galactic background. White patches in the submillimeter emission derived map are locations where the background subtraction has caused the estimated flux from the cloud to become negative in at least one wavelength. The second row of Figure 4 presents the same information, but now with a simplified color scheme for the Σ scalebar, which can be compared to regions of the Σ-PDF, as discussed below.
Figure 4. Top row, left to right: NIR+MIR extinction derived Σ map, submillimeter emission derived Σ map using SMF BG subtraction, submillimeter emission derived Σ map using GG BG subtraction. Middle row, left to right: as the top row, but with a simplified color scheme for the Σ scalebar to highlight values around the peak of the Σ-PDFs. Bottom row, left to right: submillimeter emission derived Σ map (GG BG subtraction), with β = 1.5, 2.0, and pixel-by-pixel Σ comparison between NIR+MIR extinction and fiducial submillimeter emission (GG) methods (red lines show median and 1σ range in data binned perpendicular to the one-to-one line).
Download figure:
Standard image High-resolution imageBroadly similar morphologies are seen in these maps, but with the submillimeter emission derived maps tending to find moderately higher values in the IRDC, although the differences between the SMF and GG background subtraction methods are comparable to the differences between the Submillimeter Em. (GG) model and the extinction map. Note also that the superior resolution of the extinction map.
The third row of Figure 4 shows the effect of applying different values of β for deriving Σ (GG case). Note that β = 1.8 is closest to the behavior of the OH94 thin ice mantle model. Slightly lower values of β ≃ 1.5 (which lead to warmer derived temperatures and thus lower values of Σ; e.g., Guzmán et al. 2015) are one way of reconciling differences between the extinction and emission derived maps.
Finally, Figure 4 also shows a pixel-by-pixel comparison of Σ's derived by NIR+MIR extinction and the fiducial submillimeter emission (GG) method. The fraction of the pixels with both Σs > 0.013 g cm−2 with values within 30% of each other is 0.608.
3.2. Σ Probability Distribution Functions
Figure 5(a) shows the area-weighted Σ-PDF of IRDC/GMC G028.37+00.07. Note that PDF normalization is with respect to the total "BTK14" area (see Figure 1). The result from the submillimeter dust emission map with no FG and BG subtracted (green solid line) is similar to S15b's results, being intermediate between their GMC and IRDC PDFs (normalized to their particular areas), as expected given the region geometries (Figure 1). FG-only subtraction (green dotted) has a modest effect: the PDF peak moves to lower Σ by ≃0.02 g cm−2.
Figure 5. Top row: area-weighted Σ-PDFs of IRDC/GMC G028.37+00.07. (a) Left: effects of foreground/background subtraction. Green solid line shows the Submillimeter Em. case with no FG/BG subtracted. Green dotted line shows only the FG subtracted. Brown/black solid lines show the SMF/GG BG subtraction, respectively. Blue solid line shows the NIR+MIR extinction result, including FG+BG subtraction. (b) Middle: effects of β variation. Black solid line shows the fiducial Submillimeter Em. GG case, which has effective β ≃ 1.8. The β = 1.5, 2 cases are shown with dotted purple and orange lines, respectively. Variable β model is shown with a solid purple line. (c) Right: lognormal fits to fiducial Σ-PDFs. Shaded regions correspond to the colors shown in Figure 4 (middle and bottom rows). Note that the gray region shows values below the minimum closed contour. Bottom row (d)–(f): same as top row, but for mass-weighted Σ-PDFs.
Download figure:
Standard image High-resolution image
Background subtraction (SMF: brown solid; GG: black solid) leads to larger shifts of the PDF peak, reducing its value by several factors. Overall the shapes of the two submillimeter emission derived Σ-PDFs are quite similar. The NIR+MIR derived PDF (blue solid) is also similar, especially to the PDF using GG background estimation.
Figure 5(b) illustrates effects of varying β from 1.5 to 2 on the Submillimeter Em. (GG) PDF. Smaller values of β reduce the amount of inferred high-Σ material (thus boosting the low-Σ distribution). Dust properties may vary systematically with Σ due to grain growth (LCT15). Figure 5(b) also shows a simple "Variable-β" model where β = 1.8 − 0.3(Σ/[0.5 g cm−2]), i.e., lower β in high-Σ regions broadly consistent with models of grain growth (e.g., Ormel et al. 2011), applied iteratively until convergence is achieved. The resulting Σ-PDF is quite similar to the fiducial Submillimeter Em. (GG) case.
Figure 5(c) shows lognormal fits to the Submillimeter Em. (GG) and NIR+MIR Ext. derived PDFs, with Σpeak ≃ 0.023, 0.019 g cm−2 (AV ≃ 5.1, 4.3 mag), respectively. Note that fitting is done only for Σ ≥ 0.013 g cm−2 (AV ≃ 3 mag), i.e., above the minimum closed contour in the NIR+MIR extinction map. With Σpeak greater than Σ of the minimum closed contour, the lognormal fitting is well constrained by the data. Values of 
are 0.044, 0.038 g cm−2, while σlnΣ' = 1.35, 1.15, respectively, compared to 1.4 reported by BTK14.
The Σ-PDFs are well fit by single lognormals. The fraction of mass above Σ = 0.013 g cm−2 that is in excess of the lognormal fits in their high-Σ (>0.1 g cm−2) tails is pl ∼ 0.08, 0.03 for Suubmillimeter Em. (GG) and NIR+MIR Ext. PDFs, respectively. Such small fractions may be consistent with similarly low values of ff. Krumholz & Tan (2007) estimated ff ∼ 0.02, including results from observed IRDCs, while Da Rio et al. (2014) estimated ff ≃ 0.04 in the Orion Nebula Cluster. However, detailed study of numerical simulations to link the mass fraction in these "tails," i.e., pl, with ff is still needed for trans-Alfvénic, turbulent, global clouds.
Note that these high-Σ excess fractions in the PDFs are in fact not particularly well described with power-law tails to the lognormals. Given the small size of the excess fractions, the limited dynamic range of Σ where they appear, and the potential systematic errors that enter at high Σs, we do not fit power-law functions, but rather focus on pl as our metric for deviation of the PDF from a lognormal shape.
Note also that while we consider the GG method preferable to SMF, lognormal fitting results are not too sensitive to this choice: with SMF, Σpeak decreases by 25%, σlnΣ' increases by 25%, and pl decreases by 20%.
Figures 5(d)–(f) mirror Figures 5(a)–(c), but now for mass-weighted PDFs. Values of Σpeak of the lognormal fits are 0.060, 0.051 g cm−2 for the Submillimeter Em. (GG) and NIR+MIR Ext. PDFs, respectively. Values of 
are 0.067, 0.054 g cm−2, while σlnΣ' ≃ 1.10, 0.96, respectively. These mass-weighted Σ-PDFs are also well fit by single lognormals, with the PDF peak being significantly above the minimum closed contour level.
4. DISCUSSION
We have measured the Σ-PDF from NIR+MIR extinction and submillimeter dust emission in a contiguous 30 pc-scale region centered on a dense, massive IRDC that extends to its surrounding GMC. This material is likely to eventually form a massive star cluster. The two methods give similar results, especially detecting the area-weighted Σ-PDF peak at ≃0.021 g cm−2 and the mass-weighted Σ-PDF peak at ≃0.055 g cm−2, both significantly higher than the minimum closed contour at ≃0.013 g cm−2. Comparison of extinction and emission-based methods is important to assess systematic uncertainties. The consistency of these results increases our confidence in the reliability of the Σ maps and their PDFs.
Some differences may result from dust opacity uncertainties, including potential systematic variations with density due to grain growth, which have greater influence on the submillimeter emission method. Angular resolution also leads to differences: the NIR+MIR extinction map has ∼2'' resolution, while the submillimeter emission derived map has 18'' resolution, which will tend to smooth out high-Σ peaks, thus limiting its ability to see highest-Σ regions. NIR+MIR extinction mapping faces problems of saturation at high Σ, but this should only become important at Σ ≳ 0.65 g cm−2 (BTK14): i.e., most of the range of Σ of Figure 5 is unaffected. NIR+MIR extinction mapping also fails in MIR-bright regions. Figure 4 (middle row) and Figure 5(c) show that a significant reason for difference in the amount of Σ > 0.4 g cm−2 material in the NIR+MIR extinction map compared to the Submillimeter Em. (GG) map results from the material around the central MIR-bright source in the IRDC. However, this is not enough to explain the claimed power-law tail of S15b's analysis. Rather, most of the difference of our results from S15b's are caused by our higher estimate of the diffuse Galactic plane background subtraction level.
The Σ-PDFs are well fit by single lognormals, even though this IRDC and GMC region is gravitationally bound with virial parameters of about unity (KT13; HT15). The Σ-PDF peak is greater than the minimum closed contour: i.e., peak position is not too sensitive to the choice of map boundary. If we analyze a smaller 15' × 15' region centered on the IRDC, then Σpeak ≃ 0.023 g cm−2, σlnΣ' = 1.22, and pl ≃ 0.07 (averaging results from Submillimeter Em. (GG) and NIR+MIR extinction methods), similar to the results for the 20' × 19' region. The position of this peak likely has physical significance, e.g., depending on properties of protocluster turbulence, magnetic fields, and/or feedback, e.g., protostellar outflows, and thus constrains theoretical models of star cluster formation.
Another important result is the mass fraction in the high-Σ power-law tail (or lognormal excess), pl. There is tentative evidence for a small tail being present at Σ ≳ 0.2 g cm−2 and containing ∼3%–8% of the total mass, but subject to the systematic uncertainties discussed above. Still, we consider this to be the most accurate measurement of this high-Σ lognormal excess mass fraction since we have measured Σ with two independent methods, which both detect the lognormal peak. This mass fraction also constrains theoretical models, especially protocluster star formation rate and thus duration of star cluster formation. Relatively small pl may imply small ff (Krumholz & McKee 2005; Kritsuk et al. 2011; Federrath & Klessen 2012), and thus an extended duration for star cluster formation (Tan et al. 2006; Da Rio et al. 2014). Better quantification of the relation between pl and ff should be an additional goal of star cluster formation simulations.
We thank E. Schisano for Herschel data and discussions, and S. J. Carey, C. Federrath, J. Kim, N. Schneider, and B. Wu for discussions. The comments of an anonymous referee helped to improve the article. J.C.T. acknowledges support from NASA ADAP grant NNX15AF21G.