The Single-Cloud Star Formation Relation

One of the most important and well-established empirical results in astronomy is the Kennicutt-Schmidt (KS) relation between the density of interstellar gas and the rate at which that gas forms stars. A tight correlation between these quantities has long been measured at galactic scales. More recently, using surveys of YSOs, a KS relationship has been found within molecular clouds relating the surface density of star formation to the surface density of gas; however, the scaling of these laws varies significantly from cloud to cloud. In this Letter, we use a recently developed, high-accuracy catalog of young stellar objects from $\textit{Spitzer}$ combined with high-dynamic-range gas column density maps of twelve nearby ($<$1.5 kpc) molecular clouds from $\textit{Herschel}$ to re-examine the KS relation within individual molecular clouds. We find a tight, linear correlation between clouds' star formation rate per unit area and their gas surface density normalized by the gas free-fall time. The measured intracloud KS relation, which relates star formation rate to the volume density, extends over more than two orders of magnitude within each cloud and is nearly identical in each of the twelve clouds, implying a constant star formation efficiency per free-fall time $\epsilon_{\rm ff}\approx 0.026$. The finding of a universal correlation within individual molecular clouds, including clouds that contain no massive stars or massive stellar feedback, favors models in which star formation is regulated by local processes such as turbulence or stellar feedback such as protostellar outflows, and disfavors models in which star formation is regulated only by galaxy properties or supernova feedback on galactic scales.


INTRODUCTION
In galactic disks, there is a well-established correlation between the gas mass and star formation rate per unit area, when both quantities are measured in kpc-scale or larger patches (e.g., Kennicutt 1998; Bigiel et al. 2008;Leroy et al. 2013); this correlation is known as the Kennicutt-Schmidt (KS) relation (Schmidt 1959). The correlation, however, worsens as one measures progressively smaller regions, and there is little correlation between the carbon monoxide and ionizing or far-infrared luminosities -standard proxies for gas mass and star formation rate, respectively -of individual molecular clouds or filaments 100 pc in size (Mooney & Solomon 1988;Schruba et al. 2010;Onodera et al. 2010;Kruijssen & Longmore 2014;Ochsendorf et al. 2017;Zhang et al. 2019). This apparent lack of a correlation can be the result of a true spread in the star formation rate per unit mass among clouds (Lee et al. 2016), or the failure of the proxies for mass and star formation rate. The latter is possible when we estimate star formation rate using the luminosity of massive stars, because this proxy may be under-sampled on small scales (Calzetti et al. 2012), and also underestimates the true star formation rate until the stellar population is old enough (∼5-10 Myr) to have reached a statistical steady state between the formation of new massive stars and the deaths of older ones (e.g., Krumholz & Tan 2007). Conversely, gas tracers also suffer from undersampling (Calzetti et al. 2012) and timescale issues: massive stars can rapidly disperse the gas from which they formed (e.g., Chevance et al. 2021), and if we observe a stellar population where dispersal is well underway, we will underestimate the mass of gas that was present when the stars formed. Thus analyses based on massive stars tend to underestimate the star formation rate per unit mass in young clouds and overestimate it in old clouds. When we measure the KS relation in kpc-scale patches, we average over large numbers of clouds at random ages, these errors cancel, and we recover the correct mean star formation rate per unit mass. However, the uncertainties for individual clouds might nonetheless be substantial, artificially creating scatter in the KS relation at smaller scales (Feldmann & Gnedin 2011;Kruijssen & Longmore 2014;Kreckel et al. 2018).
Whether the observed large scatter in the KS relation at small scales indicates a real scatter in star formation rate per unit mass, or whether it is simply an artifact of the observational errors described above, has profound implications for our understanding of the mechanisms by which star formation is regulated. If it is real, this suggests that the KS relation on galactic scales is due to feedback processes acting at similar scales, most likely the balance between gravity and supernovae (Ostriker & Shetty 2011;Hopkins et al. 2011;Faucher-Giguère et al. 2013), and individual clouds are free to collapse to stars with high efficiency; indeed, the lack of a KS relation within individual clouds is a direct prediction of such models (Lee et al. 2016). On the other hand, if a KS relation does hold within single clouds, particularly those containing no stars massive enough to produce supernovae, this implies that some smaller-scale or more universal mechanism inhibits star formation within individual molecular clouds. These mechanisms include turbulence, magnetic fields (Krumholz & McKee 2005;Federrath & Klessen 2012) or stellar feedback in the form of protostellar outflows, stellar winds and ionizing radiation (Krumholz et al. 2012b;Federrath 2015;Xu et al. 2020;Guszejnov et al. 2021).
A natural experiment for deciding between these possibilities is to search for a KS relation within individual molecular clouds using counts of the recently formed stars or protostars identified by their bright infrared emission from circumstellar dust. The rarity of massive stars and their disruptive effect on their host cloud means they are poor tracers on cloud scales. Protostars, by contrast, have the advantage that they sample a much shorter time interval and therefore provide a much better estimate of the "instantaneous" star formation rate (SFR), and they allow measurements of the SFR even in clouds that lack massive stars and have not been significantly affected by feedback. Studies based on this method generally do find a reasonable correlation between the number of young stellar objects (YSOs) in a cloud and its gas mass above a certain density, or its gas mass divided by its mean-density free-fall time (Krumholz et al. 2012a;Lada et al. 2012;Heyer et al. 2016;Krumholz et al. 2019). Within molecular clouds, several studies have found a power-law correlation between the surface densities of YSOs and gas (Gutermuth et al. 2011;Lada et al. 2013;Willis et al. 2015). Most recently, Pokhrel et al. (2020) used high accuracy YSO catalogs and high dynamic range gas column densities to show the presence of these laws in twelve nearby clouds. This correlation is consistent with a star formation surface density being proportional to the gas surface density squared. The scaling of this law, however, varies significantly between clouds. Moreover, their analysis technique examines the density of gas around stars on a star-by-star basis, and therefore cannot easily determine whether there is a KS relation based on the volume density of gas. Using the same data, we apply a different approach to determine the star formation law that includes a dependence on the volume density of the gas. We find that the star formation law can be recast as an effectively universal linear dependence of the surface densities between star formation rate and gas mass per free-fall time, with a very less scatter between clouds.

OBSERVATIONS
The input data for our study consist of a matched set of protostellar catalogs and cloud column density maps. We use such matched catalogs and maps for the star-forming regions Ophiuchus, Perseus, Orion-A, Orion-B, Aquila-North, Aquila-South, NGC 2264, S140, AFGL 490, Cep OB3, Mon R2, and Cygnus-X. For H 2 column density maps, we use Herschel -derived column densities. For the clouds that are <500 pc distance, we used the column density maps from the Herschel Gould Belt Survey (André et al. 2010). Full details of the data reduction procedure for the clouds that are >500 pc away are provided in Pokhrel et al. (2020), but we summarize here for reader convenience. We construct the column density maps using Herschel/SPIRE and Herschel/PACS imaging at 160 µm, 250 µm, 350 µm, and 500 µm, convolved to a common resolution. In each pixel, we fit the observed spectrum using a model for dust emission in which the free parameters are the gas column density and the temperature; in these fits the dust opacity per unit mass at 500 µm is fixed to κ 500 µm = 2.90 cm 2 g −1 based on the OH4 dust model of Ossenkopf & Henning (1994). Our column density maps are the results of these fits, and can be expressed equivalently in column of H 2 molecules, N (H 2 ), or gas mass column Σ gas ; the two are related by where m H = 1.67 × 10 −24 g is the mass of a hydrogen atom and X = 0.71 is the hydrogen mass fraction of the local interstellar medium (Nieva & Przybilla 2012). We also mask pixels where the estimated dust temperature exceeds a threshold value that indicates a Rayleigh-Jeans limit, since in this regime the column density estimate becomes very uncertain -see Pokhrel et al. (2020) for details. In the highest density regions the dust emission can be optically thick even at 500 µm and our estimation of column densities may represent the lower limits. However, we are not probing gas beyond N(H 2 ) ∼ 10 23 cm −2 , and the effect on our results is minimal. To the extent that optical depth effects are significant, they would cause us to slightly underestimate the gas mass at the highest column densities. For protostars, we use the Spitzer Extended Solar Neighborhood Archive (SESNA) catalog compiled by R. Gutermuth et al. (in preparation). SESNA is constructed using combined Spitzer IRAC (Fazio et al. 2004) 3.6, 4.5, 5.8, 8.0 µm, MIPS (Rieke et al. 2004) 24µm, and near-IR (1.24, 1.67, 2.16 µm) photometry from the Two Micron All-Sky Survey (2MASS; Skrutskie et al. 2006) spanning ∼90 deg 2 . Near-IR photometry from the UK Infrared Deep Sky Survey Galactic Plane Survey (UKIDSS GPS Lucas et al. 2008) data was used exclusively for our most distant target, Cygnus-X. Sources with excess IR emission are distinguished from field stars and further subdivided into various YSO and contaminant classifications (e.g., background galaxies and unresolved molecular hydrogen shock emission) using a series of reddening-safe color and flux selections (Gutermuth et al. 2009). With a few exceptions (Gutermuth et al. 2011;Pokhrel et al. 2020), prior work on the intracloud KS relation employed protostar identifications that required 24 µm flux measurements (e.g., Heiderman et al. 2010;Evans et al. 2014). This requirement strongly limits protostar sensitivity due to confusion with resolved nebulosity and neighboring bright sources as are found in young stellar clusters (Kryukova et al. 2014;Megeath et al. 2016;Gutermuth & Heyer 2015). SESNA and related Spitzer censuses of YSOs make robust protostar identifications that do not require 24 µm photometry, improving protostar completeness under these circumstances (Gutermuth et al. 2009;Megeath et al. 2012). In addition, SESNA has a well-measured rate of contamination from extragalactic interlopers and edge-on disks (Gutermuth et al. 2008(Gutermuth et al. , 2009, and we can therefore correct statistically for these contaminating effects. The correction procedure is explained in detail in Pokhrel et al. (2020); all our analysis in this work makes use of the statistically-corrected data.

METHODS
Given the input catalogs, we construct a series of contours within which we measure the enclosed gas mass M gas , enclosed number of protostars N PS , and enclosed area A (measured in physical rather than angular units). Our approach is similar to that explored by Heiderman et al. (2010) and Lada et al. (2010). We place the lowest contour at the lowest value of N (H 2 ) such that the resulting contour is entirely enclosed by the footprints of the SESNA catalog and the column density map. We then place additional contours with a uniform spacing corresponding to 0.5 magnitudes in A V , where for our OH4 dust model 0.5 mag of extinction in V corresponds to a gas column N (H 2 ) ∼ 5 × 10 20 cm −2 until the smallest contour does not enclose any protostar. Our estimates for the minimum and the maximum A V for each cloud are given in Table 1. The result of this procedure is a set of (M gas , N PS , A) triples for each contour level in each cloud, which forms the basis for our analysis in this work. We show our data for one example cloud, Mon R2 GMC, in Figure 1.
From our triple of directly measured quantities, we derive three additional quantities: the gas surface density Σ gas , the star formation surface density Σ SFR , and the free-fall time t ff . The first of these is straightforward: Σ gas = M gas /A. To derive Σ SFR , we adopt M PS ≈ 0.5 M for the mean mass of protostars in our catalogue (Evans et al. 2009), and the duration of the protostellar phase during which newborn stars will be included in our catalogue is t PS ≈ 0.5 Myr (Dunham et al. 2014(Dunham et al. , 2015. Consequently, we compute the star formation rate within each contour as SFR = N PS M PS /t PS , and the star formation rate per unit area as Σ SFR = SFR/A. In order to estimate the free-fall time, we follow Krumholz et al. (2012a) in computing the density of the material within each surface density contour by assuming that the unseen dimension along the line of sight is comparable to the two dimensions observed in the plane of the sky, so that ρ = 3 √ πM gas /4A 3/2 ; we then compute the freefall time as t ff = 3π/32Gρ. This amounts to assuming that the region being studied is a sphere in three dimensions.
For the best-fit analysis, we use the Orthogonal Distance Regression (ODR) method in Pokhrel et al. (2020) as well as in all the best-fit analyses performed in this Letter. Hence, biases caused by different fitting techniques when comparing the results from the two studies are minimized. In the ODR method, uncertainties in both axes are used to find the regression line that is orthogonal to the residuals in finding optimized parameters. Thus, this method is preferred over the Ordinary Least Squares method for our analysis. For the details of using ODR in astronomical datasets, see Isobe et al. (1990) and Akritas & Bershady (1996).
We estimate uncertainties on our derived quantities as follows. First, we find typical uncertainties of ∼30% in the Herschel derived column density maps and up to a factor of two uncertainty in the derived gas mass (see Pokhrel et al. 2016 for the details of uncertainty estimation). We propagated the uncertainties in column density to estimate uncertainties in derived Σ gas . For the uncertainty in the number of protostars enclosed by each N (H 2 ) contour, we assume Poissonian errors so the error on N PS is √ N PS (Khullar et al. 2019), and propagate this to obtain the uncertainty in Σ SFR . Finally, Hu et al. (2021) shows that the assumption of a uniform, spherical region that we use to estimate t ff is likely responsible for adding a scatter of ∼0.2 dex. However, because this is a systematic rather than a random error, we do not attempt to propagate it below; we defer attempts to correct for this effect to Hu et al. (2021b, in preparation).

Variation of Σ SFR with Σ gas
We begin by investigating the relationship between Σ SFR and Σ gas . In Pokhrel et al. (2020), we used the local YSO density at the location of protostars, as given by an n th nearest neighbor density to measure Σ SFR and the gas column density at that location to determine Σ gas . We found that for each of our 12 molecular clouds, Σ SFR ∝ Σ 2 gas . Thus, the analysis in Pokhrel et al. (2020) is different from the one we perform here, in that Pokhrel et al. (2020) examine the gas surface density around each protostar (i.e., a star-centric analysis), whereas here we are investigating the properties of regions defined by the clouds column density; the latter has the advantage that it allows us to investigate the dependence of the intracloud KS relation on cloud volume density. In Figure 2a, we plot the relationship between Σ SFR and Σ gas as defined by our contours. Clearly, the relation is approximately linear in log-log, and we report the best-fit results of the data to a linear functional form in Table 1. We fit only the data that comes from column density contours < 3 × 10 22 cm −2 , since above this limit the contours and number of protostars enclosed become very small, and Poisson errors in Σ SFR become large. Considering all 12 clouds, the average best-fit slope is 2.00 ± 0.27 and the average best-fit y-intercept is −4.11 ± 0.80. Near the center of the observed data range at a gas surface density Σ gas = 10 2.5 M pc −2 , the standard deviation of the measured values of log Σ SFR across all clouds is 0.30. We show the line corresponding to our average best-fit parameters, with this scatter, in Figure 2a. Individual scatter that is intrinsic to an individual cloud is not considered when calculating the standard deviation as they may be caused by observational uncertainties, while cloud-to-cloud scatter is more robust. For the star-centric approach in Pokhrel et al. (2020), we used the best-fit equations for each cloud (c.f. Table 3 in Pokhrel et al. 2020) and find the standard deviation of measured log Σ SFR to be 0.33. Furthermore, the average best-fit slope in Pokhrel et al. (2020) is 2.02 ± 0.20 and the average best-fit y-intercept is −3.88 ± 0.59. The best-fit results in these two approaches are well within 1-σ standard deviation. Such stark similarities in results using two distinct methods is strong evidence that correlation is not being biased by the method.

4.2.
Variation of Σ SFR with Σ gas /t ff While the surface densities of gas and star formation are the quantities most directly accessible from observations, most theoretical models that predict the existence of a KS relation for single clouds predict a dependence on the gas free-fall time (Krumholz & McKee 2005;Padoan et al. 2012;Federrath 2013;Krumholz et al. 2019), which depends on the volume density. Incorporating the freefall time also gives a tighter correlation when measuring the cloud-to-cloud KS relation (Krumholz et al. 2012a;Heyer et al. 2016). Inclusion of the volume density and normalization of Σ gas by the free-fall time (t ff ) is the primary difference of our approach over Pokhrel et al. (2020). Figure 2. a) log Σ SFR vs. log Σ gas for contours defined on each of the 12 sample clouds (as indicated in the legend) b) Same as (a), but using Σ gas /t ff on the horizontal axis. In both panels, black dashed lines show the median best fit relation, using the parameters shown in Table 1; for (b), the black dashed line shows the fit constrained to have a slope of unity, though the best fit for an unconstrained slope is nearly indistinguishable. The darker shaded region shows the standard deviation of the data (see Table 1) around the average best fit line. and the lighter shaded region represents two times the standard deviation. Figure 2b shows the relationship between Σ SFR and Σ gas /t ff , and Table 1 shows the results of fitting a power law relationship between these quantities. It is immediately clear that the scatter is much smaller for this relationship than for the one between Σ SFR and Σ gas alone; quantitatively, the standard deviation of the Σ SFR is reduced from 0.30 to 0.21 (computed at Σ gas /t ff = 10 2.5 M pc −2 Myr −1 ) by inclusion of the free-fall time. Moreover, the relationship is now linear, with a median best-fit slope of 0.99. This finding, coupled with the theoretical predictions for a linear relationship, motivates us to carry out a fit where we fix the slope to unity and fit only the offset, so the functional form is where ff is the fraction of the gas mass converted to stars per free-fall time. The resulting fits are indistinguishable within the error bars from those where we allow the slope to vary (see Table 1). We show the fit using the median value of ff ≈ 0.026 in Figure 2b. The reformulation of the intracloud KS relation in terms of log Σ SFR and log Σ gas /t ff has the advantage of both a linear dependence between the quantities and a lower dispersion between clouds. The linear correlation between Σ SFR and Σ gas /t ff implies that we can reformulate this as a relationship between the volume density of star formation and volume density of gas. Since the volume density is the more fundamental physical quantity for determining fragmentation scales and collapse times, and since the scatter between clouds is comparatively low, we propose that this KS law is a more fundamental, universal formulation of the intracloud KS relation.

Variation of ff with Σ gas
To demonstrate that the correlations seen in Figure 2b are not just an artifact created by comparing two quantities that are inversely proportional to the area, we further examine the value of ff as a function of column density for each of our clouds in Figure 3. We construct this figure following the method of Khullar et al. (2019), whereby we vary the contour level as shown in Figure 1, and within each contour we measure ff = SFR/(M gas /t ff ), where the values of SFR, M gas , and t ff are the values within the contour. The Figure shows how ff varies with mean gas column density within the corresponding contour Σ gas .
Note that t ff ∝ A 3/4 , while Σ gas ∝ A −1 , so if the correlation shown in Figure 2b were primarily due to the fact that both axes depend similarly on area, then in Figure 3 we would expect to find gas . Figure 3 clearly shows no such correlation, which strongly indicates that the correlation shown in Figure 2b is real rather than spurious. We also find no evidence for any threshold at which star formation becomes efficient, i.e., where ff rises substantially. This is contrary to some earlier analyses using much more limited data (Lada et al. 2010;Heiderman et al. 2010;Könyves et al. 2015). Instead we find that in almost all clouds ff is nearly constant over ≈ 1 decade in column density from ≈ 100-1000 M pc −2 , and that at column densities 1000 M pc −2 the value of ff decreases rather than increases.
The decrease in ff at high column density is contrary to the naive expectation that star formation should become more rather than less efficient in denser gas. However, it seems likely that the drop in apparent ff is not indicative of a true decline in star formation efficiency but is rather a result of one of three possible effects. One is that the YSOs we use to estimate the star formation rate and thence ff average over a finite timescale of t YSO ≈ 0.5 Myr, and this can induce bias in estimates of ff at high density. We discuss this in more detail in Section 5.1. A second possible explanation is that protostellar lifetimes might not be independent of density as we have assumed. Protostellar luminosities are observed to be higher in dense regions of molecular clouds (Kryukova et al. 2014;Dunham et al. 2014;Kirk et al. 2016), and it is possible that this is a signature of more rapid accretion that could, in turn, lead to more rapid progression through the evolutionary phase selected from our source catalog. In this case, our method would lead us to somewhat underestimate ff in the densest regions we survey.
A third possible explanation is that the densest parts of star-forming regions are also sites of bright and complex emission in the infrared that can contribute to locally reduced YSO sensitivity (Megeath et al. 2016). Our SENSA catalog is more sensitive than previous ones in these regions, as we discuss in the next section, but we may still suffer from some incompleteness in the densest regions. Again, this would cause us to underestimate ff in those regions, a feature that is consistently observed across all clouds in this analysis.   Table 1) is shown by a black dashed line.
As discussed in Section 1, the primary motivation for this study is to circumvent the biases inherent in studying star formation using tracers based on massive stars, which integrate over relatively long periods of 5 − 10 Myr, and likely alter the star-forming environment over such timescales. It is therefore important to investigate to what extent our results may suffer from similar evolutionary biases. Feedback effects from our low-mass protostars are likely small, but counting protostars still amounts to measuring the star formation rate integrated over a finite time t YSO ≈ 0.5 Myr. If either gas or stellar quantities evolve on this timescale, this could cause an error in our estimates of ff . For example, if the gas were collapsing such that the density has increased over the past ∼ 0.5 Myr, then the present-day density that we measure is higher than the mean density at the time when the YSOs formed, in which case we are underestimating t ff and thus overestimating ff . Similarly, if YSOs born inside one of our contours were to move out of it during our ∼ 0.5 Myr integration interval, then we would underestimate the SFR and thus ff .
The regions where we expect evolutionary effects to be significant correspond to those for which the free-fall time, t ff , is comparable to the integration time, t YSO ≈ 0.5 Myr. This is because the freefall time is both the fastest timescale over which gas properties are likely to change (e.g., becoming denser due to collapse), and the fastest timescale over which we expect YSO motion to be significant. 1 Thus evolutionary biases are a potential concern wherever t ff t YSO . In Figure 3, we mark the Σ gas contour at which t ff ∼ 0.5 Myr with a star. We see that, for most clouds, the decline in ff at higher Σ gas begins close to the marked point, which is strongly suggestive that evolutionary effects may be the reason that we see the decline in ff at higher Σ gas , in addition to the two observational reasons (non-constant protostellar lifetimes and incompleteness) discussed in Section 4.3. However, the converse conclusion also applies: evolution should not be a concern for lower Σ gas regions where t ff 0.5 Myr. Even if we limit ourselves to the parts of the ff curves that lie to the left of the stars in Figure 3, we still find that ff both varies little from cloud-to-cloud, and is nearly constant over almost an order of magnitude dynamic range in Σ gas .

Comparison with previous cloud-scale studies
As discussed in the introduction, we are not the first authors to search for a single-cloud KS relation (see Gutermuth et al. 2011;Lada et al. 2013;Willis et al. 2015 and so on). Using the c2d and Gould Belt Spitzer Legacy Program, Evans et al. (2014) found a correlation between Σ SFR and Σ gas /t ff , but with more scatter and a steeper slope of ∼1.47 for an ensemble of different clouds. It is therefore of interest to understand why we find a much distinct intracloud KS relation in Pokhrel et al. (2020) and in this Letter. The primary explanation for the difference is the sensitivity and depth of our Herschel -derived column density maps, with the depth of the SENSA YSO catalog as a secondary factor. To demonstrate this, we focus on the Perseus cloud as an example, and repeat our analysis using a column density map derived from extinction together with the c2d protostellar catalog (both from Evans et al. 2009); these are representative of the data quality available in earlier studies. We plot the correlation between Σ SFR , Σ gas , and Σ gas /t ff derived from these data in Figure 4; the Figure also shows our results derived from Herschel plus SESNA for comparison. The most obvious difference is that the older data cover a much smaller dynamic range -0.5 decades in Σ gas , and 1.5 decades in Σ gas /t ff , compared to 1 decade in Σ gas and 2 decades in Σ gas /t ff for our data. The difference is primarily a result of the extinction maps saturating at high column density, which prevents them from measuring the high values of Σ gas that we can probe using far-infrared dust emission (Pokhrel et al. 2016). A secondary contributor is that the SESNA YSO catalog is more complete in high-density regions. Furthermore, we have included a larger number of both low-mass star-forming clouds and high-mass star-forming GMCs in our sample and fit each cloud separately. In contrast, Evans et al. (2014) combined measurements from multiple clouds into a single fit and small differences between clouds may have affected the slope. For these reasons, the linear KS relation between Σ SFR and Σ gas /t ff apparent in our data was not favored by their analysis. Figure 4. Comparison of star-gas surface density correlation plots between our data (far-IR Herschel N (H 2 ) gas and SESNA protostars; shown in blue) and a mid-IR extinction map and together with protostars from the c2d catalog (shown in red) for the Perseus molecular cloud.

Comparison with previous Galactic-scale and extragalactic studies
In addition to comparing to previous searches for single-cloud KS relations, it is helpful to put our study in the context of whole-Galaxy and extragalactic studies of the star formation relation. We cannot meaningfully compare to studies using the traditional KS relation expressed in terms of surface densities of gas and star formation, since the surface density measured for a single, spatiallyresolved, beam-filling molecular cloud, as in our data, is not the same as the surface density measured in a beam containing many clouds with a filling factor 1, the typical situation in extragalactic studies. However, we can meaningfully compare distributions of ff , since for a collection of equaldensity clouds partly filling an observed beam, the measured value of ff for the whole beam is simply the SFR-weighted mean ff of the individual clouds. As discussed in Section 1, distributions of ff measured in Galactic-scale and extragalactic studies depend systematically on the size scale. Observations that average over regions of a few hundred pc or larger generally yield fairly small dispersions of star formation efficiency; for example Leroy et al. (2017) and Utomo et al. (2018) find σ log eff = 0.2 − 0.3 dex for galaxies in the PHANGS sample. By contrast, measurements on ∼ 10 − 100 pc scales yield contradictory results, with some reporting much larger dispersions than on larger scales; for example Lee et al. (2016) find σ log ff ≈ 0.8 − 0.9 dex for individual molecular clouds in the Milky Way, while Ochsendorf et al. (2017) find obtain σ log ff ≈ 0.6 dex in the LMC. Others, for example Vutisalchavakul et al. (2016) and Barnes et al. (2017), report smaller dispersions that are closer to those found on larger scales. Our measured dispersion, σ log eff = 0.18 dex on ∼ 1 − 10 pc scales, is firmly in the small dispersion camp, and it is interesting to ask why.
One potential explanation might be that the clouds in our study are unrepresentative of those sampled in the extragalactic or Galaxy-scale studies, which are dominated by massive star-forming regions. However, we can quickly rule out this possibility. Our sample does include a number of massive star-forming regions (Orion A, Cep OB3, Mon R2, Cygnus-X) that would be readily visible to Galaxy-scale or extragalactic studies; indeed, Cygnus-X contains > 10 6 M of molecular gas and > 2 × 10 4 YSOs, which would place it the top quartile of Lee et al.'s Milky Way GMC catalog or Ochsendorf et al.'s LMC catalog by both mass and star formation rate. Moreover, both the GMC mass function (e.g., Williams & McKee 1997;Heyer et al. 2016) and the star cluster mass function (e.g., Whitmore et al. 2014) are relatively flat, dN/dM ∼ M −1.7 and M −2 , respectively. This implies that small clouds and star clusters make a non-negligible contribution to the integrated total gas mass or star formation rate measured in partly-filled beams; quantitatively, for a GMC mass function with slope −1.7 and a mass range of 10 3 − 10 6 M , roughly 25% of the mass is contained in clouds below 10 4.5 M , which are the majority of our sample. Thus we cannot attribute the difference in σ log ff between our study and earlier Galactic-scale or extragalactic cloud studies to differences in the clouds being sampled.
Instead, a more likely explanation is that the large dispersion reported in earlier cloud-scale studies is simply an artifact of the observational errors inherent in measuring star formation rates using massive star-formation tracers, which integrate over relatively long timescales and thus sample times after which the star-forming environment has been significantly transformed by feedback. Indeed, there are already hints toward such a conclusion in the literature. For example, Gutermuth et al. (2009Gutermuth et al. ( , 2011 show that, for individual clouds, infrared luminosity does not correlate to the number of YSOs to better than an order of magnitude. Heyer et al. (2016) measure ff in ATLASGAL clumps using intermediate-mass YSO counts rather than ionizing or infrared luminosity, and find σ log ff ≈ 0.4 dex, closer to both our results and the large-scale extragalactic results than to the cloud-scale measurements using ionizing or infrared luminosity. Moreover, results based on the latter two tracers appear to depend sensitively on exactly how one assigns SFRs to individual clouds: for example, the difference in dispersion for Milky Way GMCs reported by Vutisalchavakul et al. (2016) compared to Lee et al. (2016) is almost entirely due to such differences, and in Ochsendorf et al. (2017)'s study of the LMC, simply using counts of massive YSOs ( 8 M ) rather than Hα luminosity as a star formation rate indicator, while leaving all other aspects of the analysis unchanged, reduces σ log ff by 0.1 − 0.2 dex . Together with our results here, these studies support the hypothesis that the primary explanation for the large scatter reported in some previous cloudscale estimates of ff is a failure of ionizing and IR luminosity as a tracer, rather than a physical change in the star formation process in going from galactic to cloud scales; instead, the same mechanisms regulate star formation at size scales from ∼ 1 − 1000 pc. This hypothesis will be directly testable in the next few years using JWST, which will be able to detect YSOs at substantially larger distances than Spitzer. If our hypothesis based on this study is correct, then repeating earlier cloud-scale studies using JWST -detected YSOs rather than ionizing or IR luminoisty as star formation rate indicators should yield substantially lower dispersions in ff .

On the uniformity of ff
Our sample consists of clouds whose masses and SFRs span multiple orders of magnitude, yet we find that all clouds have roughly constant ff , both from cloud to cloud and within a single cloud. This strongly suggests that star formation is regulated by local processes that are present in both lowmass and high-mass star forming regions. One candidate is magnetized, supersonic turbulence stirred and aided by feedback from low-mass stars (e.g., Krumholz & McKee 2005;Krumholz et al. 2012a;Padoan et al. 2012;Federrath & Klessen 2012). Even in low-mass star-forming regions, outflows can drive and maintain turbulence at parsec scales (Bally 2016;Offner & Chaban 2017), and modern simulations including turbulence, magnetic fields, protostellar outflows, and thermal radiation from low-mass stars -all processes that would be present even in our low-mass clouds -yield ff values of a few percent, roughly consistent with our measurements (Federrath 2015;Cunningham et al. 2018;Li et al. 2018). This is a plausible explanation for our findings. In this view, it is also possible that the bend in ff we see at the highest surface densities is associated with the transition from supersonic to subsonic turbulence (e.g., Federrath et al. 2021), since observations suggest that the surface densities at which we see the bend correspond roughly to those where the role of thermal motion in supporting the clouds begins to increase (e.g., Pokhrel et al. 2018).

CONCLUSIONS
We use Herschel-derived H 2 column density maps and the SESNA YSO catalog to explore the intracloud KS relation in star-forming molecular clouds that are <1.5 kpc away. Our main conclusions are summarized below.
1. We find that Σ SFR ∝ Σ 2 gas in all the clouds in our sample. The result is consistent with that reported by Pokhrel et al. (2020), who use a different and complementary analysis technique.
2. Incorporating volume density reduces the scatter between different clouds and reveals a linear relation: Σ SFR = ff Σ gas /t ff , where the proportionality constant ff is the free-fall efficiency.
Our results demonstrate that star formation within individual molecular clouds follows a tight KS relation at parsec scales, characterized by a linear relationship between star formation rate and mass normalized by free-fall time. This relationship is essentially the same in all the molecular clouds we studied. This is significant because the clouds themselves span a huge range of properties: for example, the Perseus and Ophiuchus clouds contain no stars with significant ionizing luminosities or winds, Cep OB3 and Mon R2 are sites of ongoing massive star formation, and Cygnus-X is comparable to large complexes observed in other galaxies. The latter two are comparable to the star-forming regions that are probed in extragalactic observations (at least for very nearby galaxies), while the former would be below the detection threshold of extragalactic star-formation studies.
The small scatter in log ff we have measured rules out models in which star formation is regulated only at galactic scales, and not within individual clouds. For example, Murray & Chang (2015) propose that molecular clouds are collapsing and that, as a result, the star formation rate within them increases with time as SFR ∝ t 2 ; Lee et al. (2016) show that the observed dispersion in log ff predicted by this model is 0.54, a factor of ≈ 3 larger than we observe. By contrast, a model in which there is no collapse and thus ff does not increase yields a dispersion in log ff of 0.16, very close to what we observe. Thus our observations strongly favor the existence of a mechanism that keeps ff close to constant across all local molecular clouds. Moreover, this mechanism must not depend on the feedback provided by massive stars such as radiation and winds, since many of the clouds we have observed contain no massive stars.
We gratefully acknowledge funding support for this work from NASA ADAP awards NNX11AD14G  The entries at the bottom of the table give the median and mean of the fit values, and the spread in the data as measured by their standard deviation evaluated near the centre of the observed data range, at Σ gas = 10 2.5 M pc −2 and t ff = 1 Myr. The entries for ff correspond to the mean and standard deviation of ff for each cloud (see Figure 3).