OBSERVATIONAL EVIDENCE OF DYNAMIC STAR FORMATION RATE IN MILKY WAY GIANT MOLECULAR CLOUDS

We present a brief description of the SFC catalog here (see Lee et al. 2012, for more detail). Lee et al. (2012), following the approach of Murray & Rahman (2010) and Rahman & Murray (2010), identified 280 SFCs in and near the Galactic plane.

Ionizing photons from young clusters can travel tens or even hundreds of parsecs through the interstellar medium (ISM) before being absorbed. It follows, and it is observed, that free–free emission regions can be much larger (∼1°–2° in diameter) than the star clusters that power them. In addition, star clusters tend to form in associations (just like stars tend to form in clusters), so that the free–free emission seen in a given direction may be powered by more than one star cluster. SFCs represent systems of clusters that are able to carve out bubbles of size ∼10–100 pc in their host GMCs. We assign a free–free flux f_ν to each SFC by dividing up the free–free flux of their host WMAP sources (the 1°–2° wide regions of free–free emission) in proportion to the 8 μm flux of each SFC, as seen by Spitzer. This is motivated by the linear correlation between free–free and 8 μm flux seen on scales of a few parsecs.

Lee et al. (2012) calculated the distance to each SFC by fitting the galactic longitude and the central local standard at rest velocity v_lsr to the Clemens (1985, hereafter C85) rotation curve. For SFCs that reside inside the solar circle (R_⊙ ≤ 8.5 kpc), they use the radial velocities of the absorption lines along the line of sight to disambiguate between near (absorption line velocity is less than that of the SFC) and far (absorption line velocity can be as large as the tangent point velocity) distances.

The rotation curve calculated by C85 is known to introduce significant errors in the distance to the objects in the outer galaxy due to Perseus arm streaming motions (Fich et al. 1989). We, therefore, update the Lee et al. (2012) distances using the rotation curve of Brand & Blitz (1993, hereafter BB93) at R₀ = 8.5 kpc and Θ₀ = 220 $\,\mathrm{km}\,{{\rm{s}}}^{-1}$.

Figure 1 shows that the distances given by the two rotation curves are generally in good agreement. The four outliers with C85 "near" distances that are close to zero have C85 solutions that are slightly inside (by ∼0.03 kpc) the solar circle, while the BB93 solutions for these objects are slightly outside (by ∼0.12 kpc) the solar circle.

Zoom In Zoom Out Reset image size

We reject solutions with galactocentric radius greater than 16 kpc, i.e., we are using a prior that there is no significant massive (O star) star formation at such large radii. The typical error in distance is ∼35%, the dominant sources of error being the streaming velocity of ionized gas, the expansion velocities of bubbles carved out by SFCs, and the motions of spiral arms (Lee et al. 2012).

We use the new all-sky catlaog of Galactic molecular clouds presented by MML16, where we describe in detail how clouds are identified. MML16 presents a number of correlations between cloud properties. We provide a brief summary here.

Clouds are identified as coherent molecular structures from the ¹²CO (J1-0) survey of Dame et al. (2001). This survey has a modest angular resolution (7.5 arcmin) compared to other observations, e.g., Jackson et al. (2006), but it has the advantage of providing a uniform data set covering the entire Galactic plane in longitude. We limit our study to −5° < b < 5°.

The identification of clouds from position–position–velocity (PPV) (or [l, b, v]) cubes is challenging, as the ISM has a fractal structure, and because unrelated clouds along a given line of sight can have similar projected velocities. Given this difficulty, it is striking that studies using different data sets and different structure identification techniques find relatively consistent scaling relations between various cloud properties, such as size, velocity dispersion, mass and surface density (see the review by Hennebelle & Falgarone 2012). MML16 finds scaling relations similar to those reported in the earlier catalogs.

MML16 employs a combination of Gaussian decomposition and hierarchical clustering analysis to identify clouds. First, every spectrum of the the entire CO PPV cube is decomposed into a sum of Gaussian components. Next, they build a cube of integrated emission W_CO where each grid cell (l, b, v) is assigned W_CO integrated over Gaussian components whose central velocities equal v. This cube is much more sparse than the original brightness cube as the integrated emission of each Gaussian component is concentrated in a single cell (l, b, v) and not spread out in velocity. This new cube facilitates the identification of coherent structures down to the noise level of the data. To do so, they used a classical threshold descent. First, they identified islands in PPV space as neighboring cells with W_CO higher than some high threshold value. The threshold is progressively lowered to 0.5 K $\,\mathrm{km}\,{{\rm{s}}}^{-1}$. At each step of this descent, new cells are revealed. They are either attached to previously identified structures or classified as new structures. This watershed method is similar in spirit to clumpfind (Williams et al. 1994) or to dendrogrammes (Rosolowsky et al. 2008). A total of 8107 coherent structures/clouds are identified over the whole Milky Way disk, recovering 98% of the total CO emission.

As for the SFCs in this paper, the distance to each GMC is estimated assuming that its average velocity follows the Galactic rotation curve of Brand & Blitz (1993). In the inner Galaxy, there is an ambiguity as two distances along the line of sight (dubbed near and far) have the same observed velocity. MML16 relied on a statistical method that has been used by many other studies (e.g., Dame et al. 1986; Solomon et al. 1987; Grabelsky et al. 1988; García et al. 2014). The original idea was to choose the distance that provides a cloud physical size, R, that matches more closely the size-linewidth σ_v ∝ R^α relation, established using clouds not subject to the distance ambiguity.

The size-linewidth relation appears to vary with the cloud column density Σ (see e.g., Heyer et al. 2009). To account for this dependence on Σ, MML16 selects the distance that best matches the relation $\sigma \propto {(R{\rm{\Sigma }})}^{0.42}$. The cloud distances chosen in this manner agree well with the distances of associated SFCs, as shown in Section 2.4.

The cloud catalog provides the position in (l, b, v), the size, the velocity dispersion, and the distance to each cloud. It also provides the gas mass estimated as

$$\begin{eqnarray}&&{M}_{{\rm{g}}}={W}_{\mathrm{CO}}^{\mathrm{tot}}\,{X}_{\mathrm{CO}}\,2\mu {{\rm{m}}}_{{\rm{H}}}\,{D}^{2}\tan {(\delta )}^{2}\end{eqnarray} \tag{ 3 }$$

where ${W}_{\mathrm{CO}}^{\mathrm{tot}}$ is the total CO emission of all the Gaussian components associated to the cloud, X_CO = 2 × 10²⁰ cm⁻² K⁻¹ km⁻¹ s (Bolatto et al. 2013), μ = 1.36 takes into account the contribution from helium, m_H is the mass of hydrogen, D is the distance to the cloud and δ = 0125 is the angular size of the pixel. We also provide the WMAP free–free and IRAS 100 μm flux of clouds, measured by integrating the emission over pixels associated with each cloud.

We identify matches between SFCs and GMCs if they meet the following criteria:

$$\begin{eqnarray}&&\sqrt{(\delta {l}^{2}+\delta {b}^{2}+\delta {v}^{2})}\leqslant 1\end{eqnarray} \tag{ 4 }$$

where

$$\begin{eqnarray}\delta l & = & ({l}_{\mathrm{GMC}}-{l}_{\mathrm{SFC}})/\max ({\sigma }_{{\rm{l}},\mathrm{GMC}},{R}_{\mathrm{SFC}}/D,0\buildrel{\circ}\over{.} 5)\\ \delta b & = & ({b}_{\mathrm{GMC}}-{b}_{\mathrm{SFC}})/\max ({\sigma }_{{\rm{b}},\mathrm{GMC}},{R}_{\mathrm{SFC}}/D,0\buildrel{\circ}\over{.} 5)\\ \delta v & = & ({v}_{\mathrm{GMC}}-{v}_{\mathrm{SFC}})/\max ({\sigma }_{{\rm{v}},\mathrm{GMC}},{\sigma }_{{\rm{v}},\mathrm{SFC}},7\,\mathrm{km}\,{{\rm{s}}}^{-1}).\end{eqnarray} \tag{ 5 }$$

By R_SFC/D we mean the angular size of the SFC measured in degrees rather than radians; ${\sigma }_{v,\mathrm{GMC}}$ is the rms velocity dispersion of a GMC, while ${\sigma }_{v,\mathrm{SFC}}$ is the half-spread velocity of SFCs (i.e., $({v}_{\max }-{v}_{\min })/2$ along the bubble walls; see Lee et al. 2012 for more detail). Out of 280 SFCs, 256 have measured median velocities v_SFC. Since the host WMAP sources are $\sim 2^\circ$ in size, we only allow SFC-GMC matches to the clouds that are large enough to encompass more than 10 pixels but smaller than $2^\circ$ in their longest axis; larger objects are likely not isolated, self-gravitating clouds.⁷ The selection criteria on the GMC size limits the cloud count to 5469. In summary, we perform a cross-correlation between 256 SFCs and 5469 GMCs.

Using the criteria given by Equations (4) and (5), we find that approximately half of SFCs are matched with multiple GMCs. Visual inspection reveals that these clouds—which are often smaller in size than SFCs—trace the outer rim of bubbles blown by their SFC counterpart when we overlay them on GLIMPSE 8 μm images. It is likely that the clouds originate from a single massive cloud that was disrupted by stellar winds and radiation pressure from its massive star clusters. These SFC-GMC "complexes" are often found inside the solar circle, so that even though multiple objects (SFCs or GMCs) coincide in the (l, b, v)-space, they can be at vastly different distances. We reject any GMC that is not within 0.3 d_SFC from the centroid of the host SFC, where d_SFC is the heliocentric distance of the SFC. For SFCs without a distance measurement, we use the distance of the GMC that is most closely matched in the (l, b, v)-space.

Last, we make a visual inspection for the luminous SFCs (those with expected M_⋆ ≥ 10⁴ M_⊙) to ensure that we have made a sensible cross-correlation. Only ∼10% of the initial SFC-GMC complexes are mismatched. Most of the mismatches stem from conflicting distance ambiguity resolution between the SFC catalog and the GMC catalog.

The fact that most of the luminous SFCs are matched with GMCs suggests a good agreement in distance ambiguity resolution between the SFC sample (resolution by absorption line velocities) and the GMC sample (resolution by a fit to a size-linewidth column-density relation), enhancing our confidence in the use of the latter technique.

We identify a known cluster from Morales et al. (2013) in the direction of SFC-GMC complexes to reassign SFC Nos. 27, 28, 93 to near distances and GMC Nos. 136, 171, 388, 398, 446, 523, 625, 812, 1312, 1656, 2733, 2734, 2816, 3422 to far distances.⁸ Distances to SFC Nos. 36 and 252 are changed to far distances after correcting for a bug in a code used to produce the catalog in Lee et al. (2012). We also manually match SFC No. 202 with GMC No. 2420, SFC No. 65 with GMC No. 2071 (associated with the W40 cluster (Mallick et al. 2013)), and SFC No. 251 with GMC No. 482 (associated with RCW 120 (Anderson et al. 2010)). Four SFCs (Nos. 31, 35, 245, 248) had no GMC match, because their distances were truncated to 12 kpc (because of their likely association with a ∼3 kpc ring; see Lee et al. 2012); we recalculated their distances.

We present the results of our cross-correlation in Table 1.

The result is 191 unique SFCs matched to 389 unique GMCs, recovering 93.5% of the total SFC free–free luminosity, 83.8% of the total SFC free–free flux, and 9% of the total GMC gas mass. All of the top 24 most luminous SFCs are matched to at least one GMC. For the rest of the paper, we identify these 389 GMCs matched to 191 SFCs as 191 SFC-GMC complexes.

Each complex inherits the sum of the CO fluxes and gas masses of the matched GMCs, as well as the mass-weighted mean Galactic coordinates l, b, velocity v, heliocentric distance d, and galactocentric radius R_gal. Following MML16, we calculate the angular size and surface area of each SFC-GMC complex by solving for eigenvalues of the moment of inertia matrix:

$$\begin{eqnarray}\phi =\left[\begin{array}{cc}{\sigma }_{l}^{2} & {\sigma }_{{lb}}^{2}\\ {\sigma }_{{lb}}^{2} & {\sigma }_{b}^{2}\end{array}\right]\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}{\sigma }_{l}^{2} & = & \displaystyle \frac{{\displaystyle \sum }_{i}{M}_{i}{\displaystyle \sum }_{c=1}^{4}{({l}_{c,i}-{l}_{\mathrm{SG}})}^{2}}{4{\displaystyle \sum }_{i}{M}_{i}}\\ {\sigma }_{b}^{2} & = & \displaystyle \frac{{\displaystyle \sum }_{i}{M}_{i}{\displaystyle \sum }_{c=1}^{4}{({b}_{c,i}-{b}_{\mathrm{SG}})}^{2}}{4{\displaystyle \sum }_{i}{M}_{i}}\\ {\sigma }_{{lb}}^{2} & = & \displaystyle \frac{{\displaystyle \sum }_{i}{M}_{i}{\displaystyle \sum }_{c=1}^{4}({l}_{c,i}-{l}_{\mathrm{SG}})({b}_{c,i}-{b}_{\mathrm{SG}})}{4{\displaystyle \sum }_{i}{M}_{i}}.\end{eqnarray} \tag{ 7 }$$

Here, i denotes each constituent MML16 cloud, M_i is the individual cloud mass, and (l_SG, b_SG) are the mass-weighted mean Galactic coordinates of the host SFC-GMC complex. The quantities $({l}_{c,i},{b}_{c,i})$ are the Galactic coordinates of the two semi-major and two semi-minor vertices (c = 1–4) of each cloud, whose semi-major axis we define as 3σ_l and semi-minor axis as 3σ_b.

Using the two eigenvalues ${R}_{\max }$ and R_min of ϕ, we define the angular radius of each SFC-GMC complex as ${R}_{\mathrm{ang}}\equiv {({R}_{\max }{R}_{\min }^{2})}^{1/3}.$ We define the velocity dispersion of the SFC-GMC complex as

$$\begin{eqnarray}&&{\sigma }_{v,\mathrm{SG}}^{2}=\displaystyle \frac{{\sum }_{i}{M}_{i}(\delta {v}_{i}^{2}+{\delta }_{+}{v}_{i}^{2}+{\delta }_{-}{v}_{i}^{2})}{{\sum }_{i}{M}_{i}},\end{eqnarray} \tag{ 8 }$$

where $\delta {v}_{i}={v}_{\mathrm{cent},{\rm{i}}}-{v}_{\mathrm{SG}}$, ${\delta }_{+}{v}_{i}={v}_{\mathrm{cent},{\rm{i}}}+{\sigma }_{v,i}-{v}_{\mathrm{SG}}$, and ${\delta }_{-}{v}_{i}={v}_{\mathrm{cent},{\rm{i}}}-{\sigma }_{v,i}-{v}_{\mathrm{SG}}$ with ${v}_{\mathrm{cent},i}$ as the central velocity of each cloud i, v_SG the mass-weighted velocity of the host SFC-GMC complex, and ${\sigma }_{v,i}$ the velocity dispersion of each cloud i. We have verified that the SFC-GMC complexes follow the ${\sigma }_{v}-{M}_{g}$ and ${\alpha }_{\mathrm{vir}}-{M}_{g}$ (where ${\alpha }_{\mathrm{vir}}\equiv 5{\sigma }_{v}^{2}{R}_{g}/{{GM}}_{g}$ is the cloud virial parameter in which ${R}_{g}=d\tan ({R}_{\mathrm{ang}})$) relations reported in MML16. The properties of the SFC-GMC complexes are presented in Table 2.

Table 2.  SFC-GMC Complexes (Sorted by Luminosities)

SFC	l	b	v	σv	Rang	Rmax	Rmin	Rgal	d	WCO	Mg
No.	(°)	(°)	(km s−1)	(km s−1)	(°)	(°)	(°)	(kpc)	(kpc)	(K $\,\mathrm{km}\,{{\rm{s}}}^{-1}$)	(M⊙)
227	336.30	−0.05	−75.34	10.80	0.35	0.41	0.33	4.56	10.79	3.94e+02	9.47e+05
228	337.90	−0.02	−55.00	7.65	0.58	0.81	0.49	5.08	11.82	2.93e+03	8.49e+06
68	30.01	−0.25	98.41	5.34	0.42	0.52	0.38	4.44	8.42	2.08e+03	2.97e+06
111	79.18	0.45	−3.21	3.48	0.95	1.00	0.93	8.77	4.25	9.64e+03	3.51e+06
274	358.63	−0.30	−1.41	9.37	0.24	0.47	0.17	6.75	15.24	4.83e+02	2.33e+06
2	0.14	−0.64	15.44	3.30	0.29	0.44	0.24	0.25	8.25	5.44e+02	7.69e+05
249	347.78	0.08	−95.96	6.88	0.33	0.42	0.30	2.70	10.32	4.23e+02	9.33e+05
110	76.56	0.25	−1.04	3.75	0.91	1.07	0.84	8.68	4.59	7.99e+03	3.43e+06
72	31.02	−0.09	102.25	5.59	0.45	0.57	0.40	4.42	6.69	1.34e+03	1.24e+06
191	311.67	0.10	−50.97	5.53	0.47	0.56	0.43	6.52	7.08	2.60e+03	2.69e+06

Note. The order of SFC Nos. appear shuffled compared to Table 1, because we adopt mass-weighted distances of matched GMCs for each SFC here.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

KM05 present a semi-analytic model of ${\epsilon }_{\mathrm{ff}}$ based on the idea of turbulence-regulated star formation. They posit that within a turbulent cloud, whose density distribution is well characterized by a log-normal distribution—with its width governed by the Mach number—only the regions of densities larger than some critical value will collapse to form stars. This critical value depends on both the virial parameter and the Mach number. Clouds with larger α_vir and ${ \mathcal M }$ are harder to collapse so the critical ρ for such clouds will be larger, leading to a smaller star formation rate per freefall time _ff.

Assuming that the log-normal density distribution established by turbulence is maintained over the cloud lifetime (but see e.g., Vázquez-Semadeni et al. 2008; Cho & Kim 2011; Kritsuk et al. 2011; Collins et al. 2012; Federrath & Klessen 2013 for the evidence of dynamic density distributions in star-forming regions), KM05 finds

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{ff},0}=0.014{\left(\displaystyle \frac{{\alpha }_{\mathrm{vir}}}{1.3}\right)}^{-0.68}{\left(\displaystyle \frac{{ \mathcal M }}{100}\right)}^{-0.32},\end{eqnarray} \tag{ 16 }$$

(their Equation (30)). We append a subscript 0 to _ff to emphasize that their model assumes that the star formation rate is time independent. We see in Figure 4 that the KM05 model expects a scatter in _ff, σ = 0.24 dex, significantly smaller than the observed σ = 0.91 dex.

The KM05 model has been criticized for characterizing an entire cloud and its collapsing substructure by a single freefall timescale. Star-forming substructures have varying densities and, therefore, different freefall times (see e.g., Hennebelle & Chabrier 2011; Federrath & Klessen 2012, and references therein). Hennebelle & Chabrier (2011) present an analytic model of _ff, which takes into account not only the different freefall time for each collapsing structure, but the recycling of turbulent flow over a dynamical time (see also Federrath & Klessen 2012; Hennebelle & Chabrier 2013):

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{ff}}\propto {e}^{3{\sigma }_{\rho }^{2}/8}\left[1+\mathrm{erf}\left(\displaystyle \frac{{\sigma }_{\rho }^{2}-\mathrm{ln}{x}_{\mathrm{crit}}}{{2}^{1/2}{\sigma }_{\rho }}\right)\right]\end{eqnarray} \tag{ 17 }$$

where σ_ρ is the spread in the density distribution, set by the cloud Mach number:

$$\begin{eqnarray}&&{\sigma }_{\rho }^{2}=\mathrm{ln}(1+{b}^{2}{{ \mathcal M }}^{2})\end{eqnarray} \tag{ 18 }$$

(b = 0.25 for purely solenoidal driving and b = 1.0 for purely compressive driving; we take b = 0.5) and x_crit ≡ ρ_crit/ρ₀ is the critical density over which substructure begins to collapse, normalized by the average bulk density of the cloud ρ₀. We experiment with two criteria for collapsing: if the local Jeans length becomes comparable to the local shock width (Padoan & Nordlund 2011, PN criterion),

$$\begin{eqnarray}&&{x}_{\mathrm{crit}}\simeq 0.067{\theta }^{-2}{\alpha }_{\mathrm{vir}}{{ \mathcal M }}^{2}\end{eqnarray} \tag{ 19 }$$

with θ = 0.35 and if the local Jeans length becomes comparable to some prescribed fraction y_cut of the cloud size (Hennebelle & Chabrier 2011, y_cut criterion),

$$\begin{eqnarray}&&{x}_{\mathrm{crit}}=\displaystyle \frac{1}{5{y}_{\mathrm{cut}}}\displaystyle \frac{{\alpha }_{\mathrm{vir}}}{{{ \mathcal M }}^{2}}(1+{ \mathcal M }{y}_{\mathrm{cut}}^{p}).\end{eqnarray} \tag{ 20 }$$

Here, p is the power-law scale from the size-linewidth relationship σ ∝ R^p; we adopt p = 0.5. The normalization of Equation (17) represents the efficiency at which a clump gas converts to a star, and we adjust it so that the median model-inferred _ff matches that of the observed. The match warrants an unusually small normalization: 0.009 for the PN criterion and 0.002 for the y_cut criterion, compared to the usual 0.02–0.05.

These multi-freefall models of star formation show stronger dependence of _ff on the Mach number compared to the KM05 model. It has been argued that the observed scatter in _ff of Milky Way clouds can be explained if the cloud Mach number or the scale of turbulence driving vary by more than two orders of magnitude (e.g., Federrath 2013; Chabrier et al. 2014). Figure 4 demonstrates that the velocity dispersion or the size of the clouds simply do not vary enough to reconcile even the multi-freefall model with the observed scatter in _ff of the MML16 clouds. No reasonable choice of y_cut can reproduce the large σ = 0.91 dex. Furthermore, HC11 models predict a strong positive correlation between _ff and the total mass, in contrast to a weak negative correlation seen in the data.

Three facts emerge from our analysis that hint at the controlling parameters of star formation other than α_vir, ${ \mathcal M }$, and the scale of turbulence driving. First, the measured scatter is significantly larger than any theoretically expected scatter. Second, there is little systematic offset between the observed and the KM05 distribution of _ff. Third, an unusually small gas-to-core efficiency is required to match the median _ff expected from the multi-freefall models to that observed. All three facts suggest the distribution of Milky Way cloud ${ \mathcal M }$ and size is far too narrow to explain the observed scatter in _ff with what the multi-freefall models predict.

Both the larger scatter and the lack of systematic offset between the observed and theoretical distribution of _ff is also evident in Figure 5. The figure shows the star formation rate per unit area ${\dot{{\rm{\Sigma }}}}_{\star }$ plotted against the ratio between the gas surface density and the freefall time, Σ/τ_ff, a comparison advocated by Krumholz et al. (2012).

Zoom In Zoom Out Reset image size

Like Krumholz et al. (2012), we find a statistically significant correlation between ${\dot{{\rm{\Sigma }}}}_{\star }$ and Σ/τ_ff. A least-squares fit produces ${\dot{{\rm{\Sigma }}}}_{\star }\propto {({\rm{\Sigma }}/{\tau }_{\mathrm{ff}})}^{0.3}$, compared to the linear relationship (the volumetric star formation law) reported by Krumholz et al. (2012).

We interpret the difference in the power-law index as the result of the large scatter in the observed value of _ff. As Figure 5 attests, the data points span more than an order of magnitude both above and below _ff = 0.01, well beyond the expected uncertainty of three predicted by Krumholz et al. (2012). We find a dispersion about the volumetric law of 0.86 dex. Similarly, large scatters in star formation rates are observed in local molecular gas from Gould's belt (Evans et al. 2014) and in ATLASGAL clumps (Heyer et al. 2016). The sample of clouds assembled by Evans et al. (2014) falls squarely within the range _ff = 0.1 and _ff = 0.001. A few ATLASGAL clumps analyzed by Heyer et al. (2016) fall below _ff = 0.001, but no clumps lie beyond _ff = 0.1. Our GMC sample is more complete than the clump samples in either the Gould's belt or ATLASGAL clumps, which may explain part of the difference in the dispersions found in the different surveys.

An alternate form of star formation law envisions cloud–cloud collision to regulate the star formation rate on approximately orbital timescales (see e.g., Tan 2000; Suwannajak et al. 2014, and references therein):

$$\begin{eqnarray}&&{\dot{{\rm{\Sigma }}}}_{\star }\simeq {\epsilon }_{\mathrm{col}}{Q}^{-1}(1-0.7\beta ){\rm{\Sigma }}/{\tau }_{\mathrm{orb}}\end{eqnarray} \tag{ 21 }$$

where _orb = 0.1 is the rate at which gas is converted to stars for each collision, $Q\sim \sqrt{{\alpha }_{\mathrm{vir}}}$ is the Toomre parameter, τ_orb the orbital time, and β is the logarithmic derivative of the velocity profile. Assuming a flat rotation curve (β = 0), we show in Figure 6 that our star-forming clouds scatter about the collision-induced star formation law with a dispersion of σ = 0.94 dex. As in Figure 5, the full min-mid or mid-max scatter is more than an order of magnitude; it is unlikely that the shear-velocity-dependent term $(1-0.7\beta )$ will vary by more than a factor of order unity.

Zoom In Zoom Out Reset image size

As with the volumetric star formation law, the large scatter about the collision-induced star formation law suggests there are extra parameters that control the rate at which gas is converted into stars.

Padoan et al. (2012) studied magnetohydrodynamic simulations of turbulent, star-forming gas, finding that _ff depends exponentially on the ratio between the freefall and the dynamical time (which is equivalent to the square root of the virial parameter):

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{ff}}\simeq {\epsilon }_{w}{e}^{(-1.6{\tau }_{\mathrm{ff}}/{\tau }_{\mathrm{dyn}})},\end{eqnarray} \tag{ 22 }$$

where _w accounts for the mass loss due to protostellar winds and outflows, and τ_dyn is the dynamical time in their simulations.

Figure 7 shows the star formation rate per freefall time for star-forming clouds, plotted as a function of τ_ff/τ_dyn, where ${\tau }_{\mathrm{dyn}}={R}_{g}/{\sigma }_{\mathrm{GMC}}$ is the dynamical time of the host cloud. We do see a hint of a decrease in the upper envelope of the distribution of _ff with increasing τ_ff/τ_dyn. However, many GMCs are more efficient at producing stars than Equation (22) would suggest. More strikingly, the _ff values of GMCs span three orders of magnitude below the prediction by Padoan et al. (2012). Semenov et al. (2016) find that the distribution of cloud virial parameter in their Galactic scale simulation is broad enough to reproduce the observed scatter in _ff using the star formation prescription of Padoan et al. (2012). It is unlikely that the scatter in _ff we find for a narrow range of the virial parameter arises from observational uncertainties. The observed virial parameter needs to scatter by extra factors of ∼30 to match our data to Equation (22). The errors in cloud radius or velocity dispersion is less than a decade (Miville-Deschenes et al. 2016). The X_CO factor presents the biggest uncertainty in the cloud mass but the expected error in X_CO is approximately 30% (e.g., Bolatto et al. 2013). We may achieve better agreement between the observed and simulated distribution of _ff once Galactic-scale simulations acquires finer resolution.

Zoom In Zoom Out Reset image size

Padoan et al. (2012) noted that _ff depended on the Alfvénic Mach number ${{ \mathcal M }}_{a}$, with the lower Alfvénic Mach number (stronger magnetic field for a given strength of turbulence) yielding smaller _ff, but only for ${{ \mathcal M }}_{a}\gtrsim 5;$ below this value, _ff will increase again (see their Figure 3). Their simple fitting formula (our Equation (22)) corresponds to the minimum star formation rates in their models with ${{ \mathcal M }}_{a}\approx 5$. Can the variation in the magnetic field strength explain the large scatter in _ff?

Variations in ${{ \mathcal M }}_{a}$ may explain those points that lie above the red dashed line (corresponding to Equation (22)) in Figure 7, but the extremely low values seen below the red dashed line cannot be explained by variations in the magnetic field strength. The Alfvénic Mach number can be estimated as ${{ \mathcal M }}_{a}\sim {\alpha }_{\mathrm{vir}}({M}_{g}/{M}_{{\rm{\Phi }}})$ where α_vir is the virial parameter and ${M}_{{\rm{\Phi }}}=0.12\pi B{({R}_{g}/2)}^{2}/\sqrt{G}$ is the magnetic critical mass where B is the magnetic field strength. SFC-GMC clouds typically have volumetric number densities of n_H ∼ 10–200 g cm⁻³ so their B ∼ 1–10 μG (Crutcher 2012); we find that these clouds have M_g/M_ϕ ∼ 10–100. The median α_vir ∼ 0.76 with a scatter of 0.32 dex. We estimate ${{ \mathcal M }}_{a}$ ∼ 4–200.

The simulations of Padoan et al. (2012) show that an order of magnitude change in ${{ \mathcal M }}_{a}$ results in a factor of ∼3 change in _ff. Taking the variations in ${{ \mathcal M }}_{a}$ of the star-forming clouds into account, we expect an upward scatter from the reference value given by Equation (22) by factors of ∼10.

All four models described in Section 6—single-freefall turbulence-regulated, multi-freefall turbulence-regulated, collision-induced, and parametrization by virial parameter—assume a star formation rate that has no explicit dependence on time. These models predict that any observed scatter in _ff should arise from variations in the internal properties of clouds (e.g., virial parameter, Mach number, Alfvénic Mach number, or freefall time) or large-scale motions. We have demonstrated above that the distribution in the observed _ff of Milky Way GMCs is far too broad to be explained by any of these models (Equations (16), (20), (21), and (22)).

A number of authors have suggested that the star formation rate on scales of GMCs and smaller increases systematically with time (see e.g., Palla & Stahler 1999; Gutermuth et al. 2011 and Murray 2011 for observational evidence or see e.g., Lee et al. 2015; Murray & Chang 2015, and Murray et al. 2015 for theoretical studies). Sampling clouds at different evolutionary stages with a time-varying star formation rate may give rise to the broad distributions in and _ff. In this section, we allow _ff to be time dependent and compare the expected width in the distribution of _ff to that observed. We note that the explicit time dependence of _ff arises from the time-varying turbulent structure (i.e., the size-linewidth relation itself changes with time due to the interplay between turbulence and gravity; see Goldbaum et al. 2011 and Murray & Chang 2015). This explicit time dependence should not be confused with the implicit time dependence portrayed by Hennebelle & Chabrier (2011), whose model accounts for the reassembly of turbulent structure that is static (i.e., the turbulent flows follow the same size-linewidth relation each time they are recycled).

We proceed to test the effect of explicit time evolution of _ff on the distribution of _ff. First, we generalize the Krumholz & McKee (2005) model to allow for time-variable _ff,

$$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{* }}{{dt}}={\epsilon }_{\mathrm{ff},0}{\left(\displaystyle \frac{t}{{\tau }_{\mathrm{ff}}}\right)}^{\delta }\displaystyle \frac{{M}_{g}}{{\tau }_{\mathrm{ff}}},\end{eqnarray} \tag{ 23 }$$

where t = 0 corresponds to the time at which the first star forms, and δ = 0 corresponds to the constant _ff model.

Next, we use this prescription for _ff in the model of Feldmann & Gnedin (2011), which describes the concomitant evolution of the stellar mass and the gas mass of a star-forming cloud:

$$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{g}}{{dt}}=-{\epsilon }_{\mathrm{ff},0}{\left(\displaystyle \frac{t}{{\tau }_{\mathrm{ff}}}\right)}^{\delta }\displaystyle \frac{{M}_{g}}{{\tau }_{\mathrm{ff}}}-\alpha {M}_{* }+\gamma ,\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{* }}{{dt}}={\epsilon }_{\mathrm{ff},0}{\left(\displaystyle \frac{t}{{\tau }_{\mathrm{ff}}}\right)}^{\delta }\displaystyle \frac{{M}_{g}}{{\tau }_{\mathrm{ff}}}.\end{eqnarray} \tag{ 25 }$$

The quantity α parametrizes the rate at which stellar feedback disrupts the host GMC, while γ is the (possibly time dependent) rate of gas accretion onto the GMC.

Equations (24) and (25) are coupled ordinary differential equations (ODEs), whose solutions are M_g(t) and M_⋆(t). Massive stars that contribute most to the ionizing radiation typically live for $\langle {t}_{\mathrm{ms},{\rm{q}}}\rangle \simeq 3.9\,\mathrm{Myr}$ so we define the mass in these "live" stars as ${M}_{\star ,\mathrm{live}}(t)={M}_{\star }(t)-{M}_{\star }(t-\langle {t}_{\mathrm{ms},{\rm{q}}}\rangle )(t)$. We can then write the star formation efficiency as $\epsilon (t)\equiv {M}_{\star ,\mathrm{live}}(t)/({M}_{g}(t)+{M}_{\star ,\mathrm{live}}(t))$ and the star formation rate per freefall time as ${\eta }_{\mathrm{ff}}(t)\equiv \epsilon (t){\tau }_{\mathrm{ff}}/\langle {t}_{\mathrm{ms},{\rm{q}}}\rangle$. Note that η_ff is analogous to the observed _ff, and is not to be confused with ${\epsilon }_{\mathrm{ff},0}{(t/{\tau }_{\mathrm{ff}})}^{\delta }$.

We illustrate the predicted evolution and the distribution of η_ff of this model in Figure 8 for three different δ = 0, 1, 2, assuming no accretion (γ = 0). We set τ_ff to the median GMC freefall time ${\tau }_{\mathrm{ff}}\sim 6.7\,\mathrm{Myr}$ and set α = 3.5 so that the cloud disperses in ∼20 Myr e.g., Williams & McKee (1997) and Murray (2011). In Section 6.5 below, we obtain a mean cloud lifetime using MML16 clouds of 21–24 Myr. To isolate the effect of δ, we assume all the host GMCs have the same α_vir and ${ \mathcal M }$, or more precisely, the same intrinsic value of ${\epsilon }_{\mathrm{ff},0}$. We account for variations in these parameters after integrating the ODEs by adding the variations in quadrature.

Zoom In Zoom Out Reset image size

The distribution of η_ff that we calculate from the model with δ = 0 (the constant _ff model) features a very sharp peak. Because the lifetime of GMCs (∼20 Myr) is substantially longer than the lifetime of live stars (∼4 Myr), the distribution in η_ff is dominated by the period in which both the live stellar mass and the host GMC mass are roughly constant with $\epsilon \sim {\epsilon }_{\mathrm{ff},0}(4\,\mathrm{Myr}/{\tau }_{\mathrm{ff}})$ and, therefore, a constant ${\eta }_{\mathrm{ff}}\sim {\epsilon }_{\mathrm{ff},0}$. The small scatter (0.16 dex) arises partly from the initial rise in the live stellar mass, and partly from the final rapid decline in GMC mass around ∼20 Myr, when the stellar feedback disrupts the host GMC.

To make a fair comparison to the observed scatter in _ff, we add in quadrature the scatter from the initial and final transients, 0.16 dex, and the dispersion in ${\epsilon }_{\mathrm{ff},0},0.24$ dex, to find a total predicted dispersion of ${\sigma }_{\mathrm{log}{\eta }_{\mathrm{ff}}}=0.29$ dex. Compared to the observed σ_br = 0.91 ± 0.22 dex, the scatter predicted by the models of time-independent _ff is too small by 2.8 standard deviations.

Allowing for a monotonic increase in _ff with time significantly broadens the distribution in and ${\epsilon }_{\mathrm{ff}}$ (see second and third rows in Figure 8). Not only does the mass in live stars monotonically increase over most of the cloud lifetime, it increases more gradually from time zero. When δ = 1 (as predicted by Lee et al. 2015 and Murray & Chang 2015), we find σ = 0.54. Adding this in quadrature to the scatter in GMC properties, we find a total dispersion of ${\sigma }_{\mathrm{log}{\eta }_{\mathrm{ff}}}=0.59$ (1.4 sigma away from σ_br). Similar calculations show a total dispersion of ${\sigma }_{\mathrm{log}{\eta }_{\mathrm{ff}}}=0.93$ for δ = 2.0 (within one standard deviation from σ_br).

Models assuming a constant star formation rate per freefall time predict too narrow a distribution in _ff. Using the simple model of Feldmann & Gnedin (2011), _ff ∝ t²—equivalently, the stellar mass increases as ${M}_{* }(t)\sim {t}^{3}$—fits the data better. Our result demonstrates that a time-varying star formation rate is in better agreement with the data compared to time-independent models. Determining the exact form of _ff that best fits the observation will require a careful parameter study, which is beyond the scope of this paper.

We note that all predicted distributions of η_ff (whether the star formation rate is assumed to be time dependent or constant) feature a deficit toward the high end and an excess toward the low end compared to the observations. The simple model of Feldmann & Gnedin (2011) assumes that feedback from stars begins to destroy the host GMC as soon as the stars form. If the feedback takes the form of gas pressure in H ii regions, or of radiation pressure, the host GMC will not be significantly affected until either pressure (or their sum) overcomes the pressure associated with the self-gravity of the GMC. Altering the model to account for this threshold effect will extend the distribution of to lower values. We leave this and similar model building efforts to future work.

We estimate the average cloud lifetime using the star-forming clouds (191 SFC-GMC complexes) and non-star-forming clouds that are massive enough—and gravitationally bound—to birth stars. We show that clouds indeed live substantially longer than ∼4 Myr. This separation of timescales is crucial if we are to use the width of the distribution of (or _ff) to distinguish between constant star formation rate models and models which allow for variations of the star formation rate with time.

Clouds that harbor at least one SFC are likely near the end of their lives: their SFCs have already carved out bubbles of size 10–100 pc. The lifetime of clouds can then be estimated by multiplying the effective lifetime of live stars by the ratio of the total number of potentially star-forming clouds to the number of clouds that have SFCs. We define potentially star-forming clouds as those that are both massive enough to birth Orion Nebula Cluster (i.e., ${M}_{\mathrm{gas}}\geqslant {M}_{* }\sim 1000\,{M}_{\odot }$) and gravitationally bound (α_vir < 3.3; see Appendix of MML16); there are 1014 such clouds. The average cloud lifetime is then $((1014+191)/191)\times 4=24\,\mathrm{Myr}$. The lifetime reduces to 21 Myr if we place the cloud lower mass limit at 10⁴ M_⊙ instead.

If the rate of star formation in the most massive (10⁶ M_⊙ and higher) GMCs does accelerate, the process of star formation must halt before  ≳ 0.1, roughly the largest value we see in such clouds. Stellar feedback—in the form of stellar winds, radiation pressure, protostellar jets, and supernovae—from massive star clusters can disrupt the natal cloud and inhibit future star formation. The idea that stellar feedback destroys GMCs was proposed long ago by Larson (1981). Galactic-scale numerical simulations (e.g., Hopkins et al. 2011 and Faucher-Giguère et al. 2013) find that the effects of stellar feedback are necessary for regulating star formation rates to the Kennicutt-Schmidt value. Turbulence in the ISM alone cannot slow the rate of star formation.

The existence of large expanding bubbles—a few to ∼100 pc wide—associated with regions bright in free–free emission is evidence for the effects of stellar feedback from young clusters (referred to as SFCs by Rahman & Murray 2010 and Lee et al. 2012). These massive young clusters, which power most of the free–free emission in the Milky Way, are identified by clear cavities in the 8 μm emission. The clusters are found in massive GMCs, but the star clusters are not enshrouded in molecular gas—they are instead in regions of ionized gas, which are enshrouded in molecular gas. In some cases, the massive stars may still be accreting gas from their natal protostellar disks, but the disks are not accreting gas from the host GMC.

We have examined our sample to see if the feedback from stars is affecting their host GMCs globally. There was no clear correlation between the variations in α_vir and the ratio of disruptive forces (gas and radiation pressure) to the binding force (dynamical pressure) of the host GMC. We also found no clear correlation between the size of the bubbles blown by SFCs in the host clouds and α_vir. It may be that the effects of radiation and gas pressure will only become evident later in the evolutionary history of the GMC, or it may be that other forms of feedback (e.g., supernovae) are responsible for disrupting the clouds. Given that GMCs in the Milky Way are found near spiral arms, combined with our finding that a large fraction of 10⁶ M_⊙ GMCs are most likely gravitationally bound (${\alpha }_{\mathrm{vir}}\lt 1$), it seems unlikely that they simply disperse on their own.