The CARMA-NRO Orion Survey: Statistical Signatures of Feedback in the Orion A Molecular Cloud

TL;DR

This study examines how stellar feedback influences the structure and kinematics of the Orion A molecular cloud by analyzing gas statistics and feedback indicators, finding limited direct correlation but some signatures suggestive of feedback effects.

Contribution

It introduces a comprehensive analysis of feedback signatures in Orion A using multiple statistical methods and compares observations with molecular cloud simulations.

Findings

No correlation between shells/outflows and gas statistics.

Anti-correlation between young star density and spectral correlation function slope.

Potential feedback signatures in principal component covariance peaks.

Abstract

We investigate the relationship between turbulence and feedback in the Orion A molecular cloud using maps of $^{12}$CO(1-0), $^{13}$CO(1-0) and C$^{18}$O(1-0) from the CARMA-NRO Orion survey. We compare gas statistics with the impact of feedback in different parts of the cloud to test whether feedback changes the structure and kinematics of molecular gas. We use principal component analysis, the spectral correlation function, and the spatial power spectrum to characterize the cloud. We quantify the impact of feedback with momentum injection rates of protostellar outflows and wind-blown shells as well as the surface density of young stars. We find no correlation between shells or outflows and any of the gas statistics. However, we find a significant anti-correlation between young star surface density and the slope of the $^{12}$CO spectral correlation function, suggesting that feedback may influence this statistic. While calculating the principal components, we find peaks in the covariance matrix of our molecular line maps offset by 1-3 km s$^{-1}$ toward several regions of the cloud which may be produced by feedback. We compare these results to predictions from molecular cloud simulations.

Tables

Table 1. Table 1: Subregion Results

Name	Right Ascension (J2000)	Declination (J2000)	${\dot{P}}_{\mathrm{shell}}$ddMomentum injection rate surface density of shells. Units are ${10}^{-3}$ M⊙ km s-1 yr-1 deg-2. Values in brackets are summed lower and upper limits from Table 3 of Feddersen et al. (2018).	${\dot{P}}_{\mathrm{out}}$eeMomentum injection rate surface density of outflows. Units are ${10}^{-3}$ M⊙ km s-1 yr-1 deg-2. Uncertainties are calculated by adding in quadrature the individual outflow uncertainties in Table 7 of Tanabe et al. (2019, submitted).	${\log }_{10}{n}_{\mathrm{YSO}}$ffSurface density of YSOs. Units are deg-2. Because of catalog incompleteness, this is likely an underestimate in the Central and OMC1 regions.	12CO ${\alpha }_{\mathrm{SCF}}$	12CO ${\alpha }_{\mathrm{SPS}}$
NorthaaSubregions from Feddersen et al. (2018)	$5^h{37}^m{35}^s$ to $5^h{33}^m{06}^s$	$-5°19\prime 26″$ to $-4°48\prime 19″$	17 [6, 41]	$6.2\pm 0.5$	3.2	$-0.096\pm 0.004$	$-3.07\pm 0.04$
CentralaaSubregions from Feddersen et al. (2018)	$5^h{37}^m{35}^s$ to $5^h{33}^m{06}^s$	$-5°34\prime 03″$ to $-5°19\prime 21″$	33 [14, 64]	$0.13\pm 0.04$	3.4	$-0.139\pm 0.002$	$-3.1\pm 0.06$
SouthaaSubregions from Feddersen et al. (2018)	$5^h{37}^m{36}^s$ to $5^h{33}^m{07}^s$	$-6°17\prime 03″$ to $-5°34\prime 03″$	41 [19, 77]	$0.7\pm 0.2$	2.9	$-0.087\pm 0.006$	$-3.43\pm 0.03$
L1641NaaSubregions from Feddersen et al. (2018)	$5^h{39}^m{01}^s$ to $5^h{34}^m{22}^s$	$-7°11\prime 44″$ to $-6°17\prime 03″$	24 [10, 52]	$2.2\pm 0.3$	2.7	$-0.059\pm 0.002$	$-3.0\pm 0.03$
NGC1977bbSubregions from Davis et al. (2009)	$5^h{35}^m{47}^s$ to $5^h{34}^m{28}^s$	$-4°57\prime 08″$ to $-4°37\prime 23″$	33 [14, 64]	$0$	3.3	$-0.121\pm 0.006$	$-3.6\pm 0.14$
OMC23bbSubregions from Davis et al. (2009)	$5^h{36}^m{03}^s$ to $5^h{34}^m{43}^s$	$-5°16\prime 37″$ to $-4°56\prime 37″$	28 [4, 85]	$19\pm 1$	3.4	$-0.138\pm 0.003$	$-3.42\pm 0.07$
OMC1ccSubregions from Ungerechts et al. (1997)	$5^h{35}^m{24}^s$ to $5^h{35}^m{06}^s$	$-5°29\prime 35″$ to $-5°17\prime 35″$	93 [33, 186]	$0$	4.2	$-0.158\pm 0.001$	$-3.83\pm 0.29$
OMC4bbSubregions from Davis et al. (2009)	$5^h{35}^m{41}^s$ to $5^h{34}^m{20}^s$	$-5°46\prime 56″$ to $-5°33\prime 00″$	162 [78, 279]	$0.19\pm 0.05$	3.4	$-0.111\pm 0.002$	$-4.46\pm 0.19$
OMC5bbSubregions from Davis et al. (2009)	$5^h{36}^m{00}^s$ to $5^h{34}^m{41}^s$	$-6°15\prime 07″$ to $-5°50\prime 06″$	18 [6, 41]	$2.9\pm 0.6$	2.9	$-0.079\pm 0.002$	$-4.04\pm 0.1$
HH34bbSubregions from Davis et al. (2009)	$5^h{36}^m{00}^s$ to $5^h{34}^m{40}^s$	$-6°40\prime 07″$ to $-6°17\prime 00″$	23 [13, 40]	$0.18\pm 0.04$	2.7	$-0.082\pm 0.002$	$-4.18\pm 0.11$
L1641NbbSubregions from Davis et al. (2009)	$5^h{37}^m{20}^s$ to $5^h{36}^m{00}^s$	$-6°42\prime 15″$ to $-6°19\prime 02″$	36 [15, 71]	$5.0\pm 0.9$	2.8	$-0.077\pm 0.002$	$-4.04\pm 0.1$
V380bbSubregions from Davis et al. (2009)	$5^h{37}^m{05}^s$ to $5^h{36}^m{00}^s$	$-6°55\prime 00″$ to $-6°37\prime 00″$	22 [4, 75]	$7\pm 1$	2.7	$-0.070\pm 0.002$	$-4.22\pm 0.12$

Equations

$$S_{jk}=S(v_j,v_k)=\frac{1}{n}\sum _{i=1}^n{T}_{ij}{T}_{ik}.$$

$$\mathrm{SCF}(\mathbf{r},\Delta \mathbf{r})=1-\sqrt{\frac{{\sum }_v[T(\mathbf{r},v)-T(\mathbf{r}+\Delta \mathbf{r},v){]}^2}{{\sum }_{v}T(\mathbf{r},v{)}^2+{\sum }_{v}T(\mathbf{r}+\Delta \mathbf{r},v{)}^2}},$$

$$\begin{array}{rl}{\alpha }_{\mathrm{SCF}}=& (-0.060\pm 0.006)\log (n_{\mathrm{YSO}}[{\mathrm{deg}}^{-2}])\\ +& (0.090\pm 0.020)\end{array}$$

Full text

ISM: clouds — ISM: individual objects (Orion A) — stars: formation — stars: pre-main sequence — turbulence

††facilities: CARMA, No:45m††software: TurbuStat (Koch et al., 2017)

\turnoffedit

1 Introduction

Stars form deep inside giant molecular clouds (GMCs) (McKee & Ostriker, 2007; Heyer & Dame, 2015). Young stars output mechanical and thermal energy, or feedback, into their birth clouds (Krumholz et al., 2014). This feedback may decrease the efficiency of star formation by counteracting gravitational collapse (Federrath, 2015) and help drive and maintain turbulence (Nakamura & Li, 2007; Offner & Liu, 2018).

GMCs are turbulent, with supersonic linewidths that increase with physical size (Zuckerman & Evans, 1974; Larson, 1981; Mac Low & Klessen, 2004). However, turbulence decays rapidly (Mac Low et al., 1998; Stone et al., 1998; Padoan & Nordlund, 1999), and must be maintained by some mechanism. Feedback from young stars may help maintain the turbulence of molecular clouds (Li & Nakamura, 2006; Matzner, 2007; Carroll et al., 2009).

Offner & Liu (2018) proposed a mechanism for translating local feedback into large scale turbulent driving. They used magnetohydrodynamic simulations to show that feedback effects may be propagated through a cloud by magnetic fields. Upon injecting winds into the simulation, they showed that the velocity dispersion outside of the wind-blown shells was increased and the velocity power spectrum flattened. These effects were caused by magnetosonic waves coming from the compressed wind-blown shells. These magnetosonic waves could explain how feedback drives turbulence at larger scales.

Feedback in molecular clouds has mostly been studied by cataloging individual features such as protostellar outflows (e.g., Arce et al., 2010; Plunkett et al., 2015), photon-dominated regions (PDRs) with far-UV heating and photoablation (e.g., Bally et al., 2018), and stellar wind-blown shells (e.g., Arce et al., 2011; Nakamura et al., 2012; Li et al., 2015; Feddersen et al., 2018). The physical characteristics of these features can then be measured to estimate the impact of feedback on a molecular cloud. However, visually cataloging feedback features is time-consuming and prone to significant bias and difficulty of separating features from the rest of the cloud.

Recently, several studies have considered the impact of feedback by measuring statistics of molecular gas structure and motions. Nakamura & Li (2007) and Carroll et al. (2009) showed that the presence of outflows modifies the velocity power spectrum of simulated molecular clouds, producing peaks at the scale where outflows inject energy into the cloud. Padoan et al. (2009) investigated the power spectrum of the molecular cloud NGC 1333, finding no evidence for a departure from a power-law near the outflow energy injection scales predicted by the above-mentioned simulations.

Swift & Welch (2008) computed the power spectrum of the red and blue line wings of 13CO in the L1551 molecular cloud, which hosts several outflows. They found a feature in these power spectra at a scale of about 0.05 pc, indicating a preferential scale of energy injection into the cloud. Sun et al. (2006) compared the CO power spectrum in different regions of the Perseus molecular cloud. They found the most actively star-forming region NGC 1333 had a steeper slope than the quiescent dark cloud L1455. Swift & Welch (2008) also measured a flat linewidth-size relationship in the low-mass molecular cloud L1551, suggesting that turbulent motions originate to a large degree at small-spatial scales. This is contrary to the turbulent cascade that is usually assumed for molecular clouds, where driving happens at large scales and dissipation occurs at the smallest scales via gas viscosity (Mac Low & Klessen, 2004).

Principal component analysis (Heyer & Schloerb, 1997) has also been used to investigate the effect of feedback on molecular clouds. Brunt et al. (2009) showed that the ratio between different-order principal components of a simulated cloud is sensitive to the driving scale of turbulence. Adding outflows to their simulations did not change this ratio, implying that feedback in the form of outflows was not driving turbulence in their simulated cloud. However, Carroll et al. (2010) showed that this analysis is biased towards the largest scales. The presence of significant smaller-scale turbulent driving from feedback may then be hidden in the principal component analysis.

The spectral correlation function quantifies the similarity of pairs of spectra as a function of their separation. Ballesteros-Paredes et al. (2002) applied the spectral correlation function to 21 cm HI spectra of the North Celestial Pole Loop region. Instead of averaging the spectral correlation function (see Equation 2 below) over all pixels, they constructed a map of the local SCF of each pixel. Their map of the SCF highlighted the edge of the expanding supernova remnant HI shell. They suggested using the SCF as a tool for finding shells.

Boyden et al. (2016) (hereafter B16) applied the TurbuStat111http://turbustat.readthedocs.io/en/latest/ (Koch et al., 2017) suite of statistical measures to molecular cloud simulations from Offner & Arce (2015). These simulations tested the effect of stellar winds on the structure of a simulated molecular cloud and successfully reproduced the expanding shells found in Perseus by Arce et al. (2011). B16 found statistical measures sensitive to the mass-loss rate of the injected stellar winds. In their simulation with winds, the covariance matrix (see Section 4.1) showed peaks of spatially correlated emission separated by 1-3 km s*-1* and the spectral correlation function (see Section 4.2) had a steeper slope compared to simulations without winds. They found the power spectrum was insensitive to stellar winds.

In this paper, we test whether statistical measures of CO in the Orion A GMC trace feedback in the cloud, as predicted by B16. In Section 2, we describe the CARMA-NRO Orion observations and split the cloud into subregions. In Section 3, we define three ways of quantifying feedback in the cloud. In Section 4, we introduce the statistics we use to summarize the CO data. In Section 5, we present the results of these statistics and compare with previous studies. In Section 6, we discuss the relationship between the statistics and feedback in the cloud. In Section 7, we summarize the conclusions of the paper and suggest future directions for the statistical study of feedback.

2 Data

2.1 CARMA-NRO Orion CO Maps

We use the 12CO(1-0), 13CO(1-0), and C18O(1-0) spectral-line maps of Orion A from the CARMA-NRO Orion Survey (Kong et al., 2018). These maps were obtained by combining interferometric images from the Combined Array for Millimeter/Submillimeter Astronomy (CARMA) with single-dish maps from the Nobeyama Radio Observatory (NRO) 45m radio telescope. This method preserves the angular resolution of CARMA while also recovering large-scale structure. The combined maps probe physical scales of 0.01 - 10 pc at a distance of 414 pc (Menten et al., 2007).222The part of Orion A covered by our survey has been located about 5% closer (at about 380-400 pc) with GAIA Data Release 2 parallaxes (Kounkel et al., 2018; Großschedl et al., 2018; Kuhn et al., 2019). Adopting a distance of 390 pc would decrease the projected distances and the derived physical quantities (momentum injection rates) in this paper by a factor of 0.94. Our conclusions are independent of distance. The 12CO and C18O maps have a beam full-width at half-maximum (FWHM) of $10\arcsec\times 8\arcsec$ and the 13CO map has a FWHM of $8\arcsec\times 6\arcsec$. The 12CO velocity resolution is 0.25 km s*-1* while 13CO and C18O have a velocity resolution of 0.22 km s*-1*. For more details on the CARMA-NRO Orion data, see Kong et al. (2018).

2.2 Subregions

To compare the statistics of turbulence with feedback impact, we would ideally like to have a control - a cloud which is identical in every way to Orion A except with no stellar feedback. In this ideal case, any statistical differences between the clouds could be attributed solely to the action of feedback. This ideal scenario was simulated by B16 but no such true control cloud exists for Orion A, though Lada et al. (2009) proposed that the California Molecular Cloud is an Orion A analogue with an order of magnitude lower star formation rate. In this study, we compare different regions of Orion A with different amounts of feedback, assume that this is the only relevant difference, and look for trends in the statistical methods that resemble those found in the simulated clouds of B16. A similar approach is used by Sun et al. (2006) in their comparison of power spectra in different regions of the Perseus molecular cloud.

We divide the Orion maps into several subregions, guided by previous studies. Feddersen et al. (2018) split Orion A into four subregions - North, Central, South, and L1641N - to compare the impact of expanding shells with protostellar outflows and cloud turbulence. The North subregion includes the NGC 1977 PDR (Peterson & Megeath, 2008), OMC-2/3, and the M43 HII region. The Central subregion covers a wide variety of environments, including the Orion Bar PDR (Goicoechea et al., 2016), OMC-1 (the densest part of the cloud), the explosive Orion BN/KL outflow (Bally et al., 2017), and more diffuse gas to the east and west. The South subregion covers OMC 4/5 and the pillar-shaped PDRs to the east dubbed the dark lane south filament (DLSF) by Shimajiri et al. (2011). The L1641N subregion covers the low-mass cluster L1641-N (Nakamura et al., 2012) and the reflection nebula NGC 1999 (Stanke et al., 2010).

In addition to these regions, we consider several smaller subregions focused on individual parts of the cloud based on the molecular hydrogen outflow survey of Davis et al. (2009). These smaller subregions focus on specific clusters or groups of young stars in the Orion A cloud and all have similar projected areas on the sky of about 0.1 deg2 (5 pc2). We also add a subregion defined by Ungerechts et al. (1997) which is restricted to the densest core of Orion A - OMC 1, BN-KL, and the Orion Bar. We define these subregions following conventional definitions in Orion A to avoid cherry-picking regions that happen to show the correlations between statistics and feedback that we are looking for. Figure 1 shows the subregions on a map of 12CO peak temperature and Table 1 defines the extent of each subregion.

3 Quantifying Feedback

In order to relate gas statistics to the impact of feedback in the cloud, we attempt to measure this impact. B16 quantified feedback using the mass-loss rate of stellar winds injected into their simulations. Offner & Arce (2015) designed these simulations to reproduce the expanding shells observed by Arce et al. (2011) in the Perseus molecular cloud. Thus, to make the closest comparison to B16, we first consider the expanding shells in Orion A identified by Feddersen et al. (2018).

3.1 Expanding Shells

Feddersen et al. (2018) identified 42 expanding shells in Orion A using the CO maps from the NRO 45m telescope. The authors visually identified expanding structures in the CO channel maps and matched many of these structures with low- and intermediate-mass young stars. Similar shells have been found in the Perseus (Arce et al., 2011) and Taurus (Li et al., 2015) molecular clouds. While the origin of such shells is unclear, one explanation is spherical stellar winds from young stars which entrain the cloud material into expanding shells. We show the shells of Feddersen et al. (2018) in Figure 1.

The impact of shells on the cloud is summarized by Table 3 in Feddersen et al. (2018). In this study, we quantify a shell’s impact on the cloud by its momentum injection rate ($\dot{P}_{\rm shell}$). The total shell momentum injection rate in each subregion is the sum of $\dot{P}_{\rm shell}$ for each shell centered inside that subregion. The true shell impact is uncertain, as it is difficult to extract shell emission cleanly from the rest of the cloud. Furthermore, most shells are incomplete. They are only detected over a fraction of the expected volume in the spectral cube. Lower and upper limits on $\dot{P}_{\rm shell}$ are estimated in Table 3 of Feddersen et al. (2018) and can span a factor of several above or below the median value. We ignore these lower and upper limits and instead focus on the relative shell impacts between subregions, but the uncertainties mean any results based on the shell momenta are inconclusive. Because of these large uncertainties and the potential for false-positive bias in identifying expanding shells, we use independent methods to quantify feedback.

3.2 Protostellar Outflows

Young stars that are actively accreting material launch collimated bipolar outflows which impact the surrounding cloud material and may help drive or maintain turbulence (Arce et al., 2007; Frank et al., 2014). In Orion A, outflows have been identified by many authors (e.g. Morgan et al., 1991; Williams et al., 2003; Stanke & Williams, 2007; Takahashi et al., 2008). A comprehensive census of outflows found in our NRO 45m data has been carried out by Tanabe et al. (2019, submitted). They have detected about 50% more outflows than previously known. Notably, they identified 11 outflows in the poorly studied OMC 4/5 region where none were previously known. Tanabe et al. (2019, submitted) excluded the area around OMC-1 from their search because the YSOs in this region are crowded and the cloud velocity width is too broad to disentangle outflow emission. Therefore, the outflow measurements for our regions around OMC1 are incomplete.

To quantify the impact of outflows on the cloud, we use the momentum injection rates ($\dot{P}_{\rm out}$) tabulated in Table 6 of Tanabe et al. (2019, submitted). Many of the outflows in this table have multiple lobes listed separately. In these cases, we sum the individual lobes of each outflow. We assign outflows to the subregion containing the driving source in Table 3 of Tanabe et al. (2019, submitted). The only case where the driving source of an outflow lies in a different subregion from part of the outflow emission is Outflow 19 in Tanabe et al. (2019, submitted). A small portion of the emission from this outflow falls south of the OMC23 subregion, where its driving source is located. In every other case, the outflow emission and its driving source are in the same subregion.

To measure the outflow mass, Tanabe et al. (2019, submitted) first integrate the 12CO emission. 12CO is often optically thick, which means line intensity is no longer directly proportional to column density and mass. To account for the optical depth of 12CO, other tracers must be observed. For example, Zhang et al. (2016) use 13CO and C18O combined with 12CO to measure outflow masses more accurately. Tanabe et al. (2019, submitted) only detect a few outflows in both 12CO and 13CO. They calculate the average optical depth of these few outflows and apply this correction factor to every outflow in their catalog, ignoring any variation in optical depth.

After they find the mass of an outflow, Tanabe et al. (2019, submitted) they calculate momentum by multiplying this mass by the line-of-sight velocity of the outflow. But if the outflow axis is inclined relative to the line-of-sight, this will underestimate the true momentum.

Variable optical depth and inclination angle both introduce uncertainty in the outflow momentum injection rates. The optical depth of outflows in Orion A is likely not uniform. This fact is evidenced by the different outflow optical depths measured by Tanabe et al. (2019, submitted) as well as the variations in the ratio between the integrated intensities of 12CO and 13CO throughout the cloud (see Figure 25 in Kong et al., 2018). A variable outflow optical depth will affect our comparisons of outflow impact between subregions.

3.3 Young Stars

As an independent estimate of feedback we consider the young stars (YSOs) in each subregion. YSOs are more directly a measure of star formation rate than feedback impact. However, YSOs are ultimately responsible for both outflows and shells. Therefore we consider the surface density of YSOs $n_{\rm YSO}$ to be a proxy for the relative strength of feedback in different regions.

To measure $n_{\rm YSO}$, we use the Spitzer Orion catalog of YSOs (Megeath et al., 2012). They classified stars as protostars or pre-main sequence stars with disks on the basis of their mid-IR colors. In this classification, 86% of the Spitzer Orion YSOs are pre-main sequence stars and the remaining 14% are protostars. Shells are likely driven by pre-main sequence stars (Arce et al., 2011; Feddersen et al., 2018) while outflows are more likely driven by protostars. Therefore, we calculate $n_{\rm YSO}$ using all YSOs in the catalog.

Because the YSO catalog is a mixture of more evolved pre-main sequence stars and younger protostars, this measure traces a wider range of timescales than either shells or outflows alone. As noted in Section 5.2.3 of Arce et al. (2010), the cumulative impact of feedback in a cloud over the course of star formation may be greater than what is traced by the currently active outflows and shells. Thus, if there are YSOs which have ejected outflows or powered shells no longer detectable as coherent structures in the cloud, then $n_{\rm YSO}$ may be a better tracer of the potential link between feedback and turbulence than shells or outflows alone.

To compare subregions of different sizes, we calculate the surface density of each feedback measure, dividing by the projected area of each subregion. The Spitzer Orion catalog suffers from incompleteness toward regions with bright IR nebulosity (Megeath et al., 2016). Thus the YSO surface density in the Central and OMC1 subregions is likely higher relative to the surface density in other subregions. From Figure 2 and 3 in Megeath et al. (2016), the typical completeness fraction in the ONC is approximately 0.5. Therefore, the true central YSO surface density may be up to about twice the value reported here. Our three feedback measures are the momentum injection rate surface density of shells/outflows (M⊙ km s*-1* yr*-1* deg*-2*) and the YSO surface density (deg*-2*).

4 Statistical Methods

We use the TurbuStat333https://turbustat.readthedocs.io/en/latest package, described in detail by Koch et al. (2017), to compute the principal component analysis, spectral correlation function, and spatial power spectrum of the 12CO, 13CO, and C18O data in Orion A. B16 also found several other statistics that were sensitive to wind mass loss rate (see their Table 3). However, Boyden et al. (2018) incorporated astrochemical models into these simulations and found that many of the statistics were also sensitive to the chemical complexity, radiation field, and molecular tracer used. The statistics which are least sensitive to these effects while remaining sensitive to feedback strength were principal component analysis and the spectral correlation function (see Section 5.3 in Boyden et al. 2018). While the spatial power spectrum was not sensitive to feedback in B16, its form has been studied for possible signatures of feedback (e.g. Swift & Welch, 2008) so we include it here. TurbuStat also provides a distance metric for each statistic, which allows comparison of different cubes. We briefly describe these statistical methods below.

4.1 Principal Component Analysis

Principal component analysis (PCA) is a statistical technique used to reduce the dimensionality of a dataset. Ungerechts et al. (1997) and Heyer & Schloerb (1997) first applied PCA to molecular line maps to study the chemistry and turbulence in molecular clouds.

The method implemented by TurbuStat to compute the PCA of a spectral cube comes from Heyer & Schloerb (1997). A spectral cube with $n$ pixels can be expressed as a set of $n$ spectra with a number $p$ of velocity channels. This can be represented as a matrix $T(r_{i},v_{j})\equiv T_{ij}$ with $n$ rows and $p$ columns, where $r_{i}$ is the position of pixel $i$ and $v_{j}$ is the velocity of channel $j$.

First, each element of the covariance matrix $S$ is calculated to be

[TABLE]

This covariance matrix is then diagonalized to find its eigenvalues and eigenvectors. Projecting the spectral cube onto each eigenvector gives a set of eigenimages called principal components. Each eigenvalue correspond to the amount of variance recovered by each principal component.

B16 found the covariance matrix to be sensitive to the injected stellar winds. Their Figure 3 shows peaks at 1-3 km s*-1* in the covariance matrix of the simulation with winds which are absent in the simulation without winds. We test whether such features appear in the Orion covariance matrices.

To compare their simulated spectral cubes, B16 used the Turbustat PCA distance metric. Each set of eigenvalues for a particular cube is sorted in descending order; then the eigenvalues are normalized by dividing them by their sum. The PCA distance metric between two cubes is then the Euclidean distance between the two normalized sets of eigenvalues, or the square root of the sum of the square differences. B16 found a strong correlation between the PCA distance metric and winds. Boyden et al. (2018) incorporated gas chemistry into these simulations and found that PCA also varied between models with and without chemistry. However, the covariance peaks remained a unique signature of winds. For more details on PCA, see Brunt & Heyer (2002a, b, 2013).

4.2 Spectral Correlation Function

The spectral correlation function (SCF) was first introduced by Rosolowsky et al. (1999) and refined by Padoan et al. (2001). The SCF measures the similarity of spectra as a function of their spatial separation, or lag. The SCF at a specific lag vector $\Delta\mathbf{r}$ (between two pixels) is defined to be

[TABLE]

where $\mathbf{r}$ is the position of a pixel, $\Delta\mathbf{r}$ is the lag vector, $v$ is velocity, and $T$ is the temperature (or intensity). The SCF is then averaged over all pixels $\mathbf{r}$. Repeating this for various lag vectors, a 2D spectral correlation surface can be constructed. An azimuthal average of this surface, or equivalently an average over all rotated lag vectors of the same length, is a 1D spectrum of the SCF as a function of spatial separation. We refer to this 1D spectrum as the SCF in this paper. Padoan et al. (2001) showed that the SCF is well characterized by a power-law in both simulated and observed molecular clouds, over a wide range of physical scales. They also showed that the slope of the SCF is independent of velocity resolution and signal-to-noise.

The SCF distance metric is defined in TurbuStat as the sum of the square differences between the SCF surfaces, weighted by the inverse square distance from the center of the surface. B16 found SCF to be sensitive to the strength of feedback. In their simulations with winds injected, the SCF has a significantly steeper slope than in the wind-free simulations. More recently, Boyden et al. (2018) found that including chemistry in these simulations flattened the SCF slope as opposed to the steepening seen with winds. For more details on SCF, see Rosolowsky et al. (1999) and Padoan et al. (2001).

4.3 Spatial Power Spectrum

The power spectrum is the Fourier transform of the two-point autocorrelation (or square) of the integrated intensity. We azimuthally average this two-dimensional power spectrum to arrive at a one-dimensional power spectrum which we hereafter refer to as the SPS. See Stutzki et al. (1998) for the general $n$-dimensional derivation and Pingel et al. (2018) for a detailed description of the SPS implementation used here.

Because our maps have emission at the edges, the Fourier transform is affected by strong ringing (e.g. Brault & White, 1971; Muller et al., 2004) which can be seen at small spatial scales in the power spectrum. To correct for this ringing, we taper the integrated intensity with a Tukey window where the outer 20% of the map is gradually reduced to zero. This tapering also reduces the noise at the edges of the observed maps. Additionally, the observed beam introduces artificial correlation into the map (e.g. Dickey et al., 2001). To correct for this, we divide the power spectrum of the integrated intensity by the power spectrum of the ellipsoidal Gaussian beam. This introduces a divergence at very small scales but allows us to extend the dynamic range over which the underlying power spectrum is recovered. These corrections are implemented in TurbuStat and described in tutorials included in the package documentation.

B16 found that the SPS was not sensitive to feedback in their simulations of winds, while Boyden et al. (2018) found that temperature variation flattened the slope. However, some studies have suggested that feedback induces a break in the power spectrum. Swift & Welch (2008) reported a “bump” in the 13CO SPS of L1551 at a scale of about 0.05 pc. They attributed this peak to the energy injection scale of the outflows in the cloud. However, they used the power spectrum of the line wings which may be more sensitive to outflow emission than the full integrated intensity of the cloud. We discuss the line wing power spectra in Section 6.2.

5 Results

5.1 Principal Component Analysis

B16 showed that the covariance matrix (Equation 1), an intermediate product of PCA on spectral cubes, was sensitive to the presence of feedback. In their simulations with winds included, the covariance matrix shows several peaks at $1-3$ km s*-1* which are not present in the simulations without winds. In Figure 2, we show the covariance matrices of the Orion A 12CO subregions (the 13CO and C18O covariance matrices are described in Appendix A). In several regions, most prominently NGC 1977, L1641N, OMC 4, and V380, we find covariance peaks offset from the diagonal axis by $1-3$ km s*-1* as in the B16 wind simulations. Essentially, these covariance peaks mean there is emission in these subregions which is spatially correlated but separated by a few km s*-1* in velocity. This spatially correlated (but velocity-offset) emission can clearly be seen in the channel maps of these regions (Kong et al., 2018). The source of these features could be feedback, as it is in the B16 simulations, or some other effect.

Nakamura et al. (2012) proposed that the L1641N cluster is located at the intersection of two colliding clouds. This idea is based on the overlapping velocity components in this region. The blue velocity component (4-6 km s*-1*) dominates emission to the southeast of L1641N while the red component (7-12 km s*-1*) extends north of the cluster. Nakamura et al. (2012) noticed these velocity components in channel maps, but these are the same velocity components that appear as peaks in our covariance matrices of the L1641N and V380 regions. It is unclear whether the covariance peaks are the result of this cloud collision scenario or the expansion of wind-blown shells like those simulated by B16 or some combination of the two effects.

In the NGC 1977 region, the molecular cloud is excited by FUV radiation from the H II region to the north (Kutner et al., 1985; Makinen et al., 1985). This FUV may be responsible for photoablation (Ryutov et al., 2003) of the northern edge of the molecular cloud, accelerating it to a few km s*-1* and generating the covariance peaks seen in this region. Well known examples of photoablative flows can be found in the Orion Bar (Goicoechea et al., 2016) and the Horsehead Nebula (Bally et al., 2018). Detailed comparison of the covariance matrix in both simulated and observed molecular clouds is needed to fully understand the mechanisms behind these peaks.

B16 used the PCA distance metric to further show that PCA was sensitive to the winds in their simulations (see their Figure 16). We find no correlation between any of our feedback measures and the PCA distance metric between subregions. This could be because our feedback measures are not good proxies for the winds simulated by B16, or because another mechanism unrelated to feedback is driving the velocity structure traced by the covariance peaks.

5.2 Spectral Correlation Function

Figure 3 shows the 12CO SCF and power-law fits of each subregion. The fit slopes are tabulated in Table 1. We compute the two-dimensional SCF surface at lags between 0 and 30 pixels (0 to 0.12 pc) in intervals of 3 pixels (0.01 pc). Using a lag interval of 1 pixel does not significantly change the resulting SCF, so we save computational time by only calculating the SCF surface every 3 pixels. We then average the SCF surface in equally-spaced annuli to arrive at the one-dimensional SCF spectrum shown in Figure 3. We calculate the SCF over the same range of spatial scales as shown in Figure 6 of B16 for the most direct comparison. Each subregion’s SCF follows a power-law closely up to a lag of approximately 20 pixels. We compute a weighted least-squares fit to each SCF between lags of 5 to 17 pixels (approximately 0.02 to 0.07 pc). Unlike in the B16 simulations, the Orion SCF steepens at larger lags. Gaches et al. (2015) also found a steepening SCF in 13CO maps of Ophiuchus and Perseus. We describe the SCF at scales larger than 0.12 pc in Appendix C.

The SCF slopes in the Orion A subregions range from -0.15 to -0.06. Gaches et al. (2015) calculated the 12CO SCF of the Perseus and Ophiuchus molecular clouds and found slopes around -2, steeper than in any of our Orion A subregions. However, limited by the angular resolution of their data they fit the SCF at larger scales: 0.1 - 1 pc. Padoan et al. (2003) found 13CO SCF slopes between -0.1 and -0.5 in various molecular clouds. Most of their SCF spectra are fit at larger scales ($>0.1$ pc) than those presented here. But the smallest maps in their dataset, L1512 and L134a, have similar spatial resolution to our Orion A maps and also have the shallowest SCF slopes at -0.18 and -0.13, respectively. The simulations in B16 also have steeper SCF slopes compared to our data. The SCF is well known to vary with spatial resolution (Gaches et al., 2015) making comparison between different datasets difficult. B16 found the SCF slope was sensitive to their simulated wind mass loss rate. In Section 6.1 we compare the SCF slope between subregions and discuss the relationship between SCF slope and feedback impact.

5.3 Spatial Power Spectrum

Figure 4 shows the 12CO power spectra and power-law fits for each subregion. The fit slopes are tabulated in Table 1. The beam correction described in Section 4.3 causes the sharp upturn at scales smaller than about twice the beamwidth. We restrict the power-law fits to spatial scales greater than about five times the beamwidth. We also exclude the largest scale point in each power spectrum from the fits as some of the power spectra show slight deviations at this largest scale. In the fitted regime, the power spectra closely follow a power-law with no evidence of the peaks seen in Swift & Welch (2008).

The power-law fits to the SPS in the Orion A subregions have slopes ranging between about -3 and -4 in 12CO and 13CO with somewhat shallower slopes between about -2.2 and -3.4 in C18O (see Appendix A for 13CO and C18O SPS). The SPS slope of an optically thick medium is predicted to saturate to -3 for a wide range of physical conditions (such as sound speed and magnetic field strength) (Lazarian & Pogosyan, 2004; Burkhart et al., 2013). Previous studies of the 12CO and 13CO SPS in molecular clouds have shown slopes close to this optically thick limit of -3 (Stutzki et al., 1998; Padoan et al., 2006; Sun et al., 2006; Pingel et al., 2018, e.g.), shallower than our SPS slopes. However, these studies measure the SPS of entire clouds instead of smaller regions within clouds, averaging over larger areas than our subregions. Sun et al. (2006) reported the power spectrum slope in several regions within the Perseus Molecular Cloud spanning $50\arcmin\times 50\arcmin$, or $4.4~{}\mathrm{pc}\times 4.4~{}\mathrm{pc}$ at a distance of 300 pc (Zucker et al., 2018; Ortiz-León et al., 2018). Most of these regions have SPS slopes steeper than -3, more similar to our subregions than to cloud-wide power spectra. Also, our largest subregions have slopes closest to the theoretically predicted value of -3. In Section 6.1, we compare the SPS slope to feedback impact in each subregion.

6 Discussion

6.1 Spectral Slopes and Feedback Impact

In the simulations of B16, SCF responds strongly to the strength of feedback but is independent of evolutionary time and magnetic field strength. They found that SCF slope steepened with increased wind mass-loss rate. Boyden et al. (2018) showed that SCF is sensitive to gas chemistry, but including chemistry flattens the SCF slope - the opposite impact of feedback. Thus, SCF may still probe the relative strength of feedback in different regions, especially if chemistry is similar between regions.

SPS, on the other hand, has only a weak dependence on feedback in B16 but varies strongly with evolutionary time and magnetic field strength, making SPS a poor diagnostic for feedback. Further, Boyden et al. (2018) showed that SPS is also sensitive to temperature variations in their simulations which include chemistry.

We quantify both SCF and SPS by the slope of their power-law fits which are shown for 12CO in Figure 3 and Figure 4. Figures 5, 6, and 7 plot the SCF and SPS slopes in each subregion against the shell momentum injection rate, outflow momentum injection rate, and YSO surface densities, respectively. We look for systematic trends between our feedback impact measures and either statistic.

In our most direct comparison to B16, we find no correlation between the SCF or SPS slopes and shell momentum injection rate surface density in the Orion A subregions (Figure 5). In particular, we do not find the predicted SCF steepening with increased wind feedback. Using only those shells with the highest confidence score (i.e., a score of 5) from Table 2 in Feddersen et al. (2018), there is still no significant trend between shell momentum injection rate and SCF or SPS slope. We also find no correlation between the spectral slopes and the outflow momentum injection rate surface density (Figure 6). Because both of these feedback measures are quite uncertain (see Section 3), this analysis does not rule out an underlying correlation between the momentum injection rates and the spectral slopes. Using the more objective measure of YSO surface density, we do find a significant correlation with SCF slope.

The 12CO SCF steepens in subregions with higher YSO surface density (Figure 7). While B16 does not model YSO populations, they show that increased feedback strength steepens the 12CO SCF in their simulations. If the true feedback impact is positively correlated with $n_{\rm{YSO}}$, then our result is consistent with the SCF prediction by B16. We stress that feedback is not the only possible driver of the relationship between $n_{\rm{YSO}}$ and 12CO SCF. We tested the effects of column density on these relationships using the Herschel maps from Stutz & Kainulainen (2015). The median column density of a subregion does not affect the trends shown in Figures 5 through 7. However, some other underlying variable that correlates strongly with stellar density and influences 12CO emission (e.g., opacity or excitation temperature) may still explain the trend we see between $n_{\rm{YSO}}$ and 12CO SCF.

We fit the 12CO SCF slopes with a weighted least-squares regression, shown in the upper left panel of Figure 7. The best-fit line is

[TABLE]

where $\alpha_{\mathrm{SCF}}$ is the 12CO SCF slope. This fit has a correlation coefficient of $r^{2}=0.84$. If we exclude the OMC-1 subregion, which has an order of magnitude higher $n_{\mathrm{YSO}}$ than any other subregion, the slope of the best-fit line steepens slightly but the correlation coefficient does not change. Future studies of the SCF in molecular clouds should test this correlation.

We do not find any correlation between SPS slope and $n_{\rm{YSO}}$, which is also consistent with B16. In 13CO and C18O, neither SCF nor SPS appear to be correlated with $n_{\rm{YSO}}$. It is unclear why the 12CO SCF slope is correlated with $n_{\rm{YSO}}$ while the other CO lines are not.

When measuring the SCF slopes, we fit the SCF up to about 0.12 pc scales. We chose this fitting range to be most consistent with B16. The SCF steepens toward larger scales. However, fitting the SCF at larger scales does not change the relative trends between SCF slope and feedback impact reported in Figures 5 through 7. In particular, the anticorrelation between 12CO SCF slope and $n_{\rm{YSO}}$ remains. We show the SCF at larger scales in Appendix C.

Figure 8 shows the three feedback measures plotted against each other. There is no significant correlation between any of the three feedback measures among the Orion A subregions. This lack of correlation is unsurprising. As discussed in Section 3, shells, outflows, and YSOs each trace different populations of forming stars. For this reason, along with the significant uncertainties in the feedback measures also described in Section 3, it is premature to conclude that the feedback mechanisms are not spatially correlated.

6.2 Line Wing Power Spectrum

The power spectra in Figure 4 are calculated using the 12CO integrated intensity maps. Because outflows are often most prominent in the line wings (channels blueward and redward of the main cloud velocity range), their influence may be largest in power spectra restricted to these velocity ranges.

Swift & Welch (2008) studied the line wing power spectrum of the low-mass star-forming cloud L1551. They found a peak at an angular scale of 1′ (0.05 pc) in the power spectrum of 13CO integrated over the line wings. They attributed this feature to a preferential scale at which protostellar outflows deposit energy into the cloud.

Padoan et al. (2009) showed that the integrated intensity power spectrum in the NGC 1333 molecular cloud is nearly a perfect power law with no features or flattening despite the many outflows present (Plunkett et al., 2013). Padoan et al. (2009) interpret this to mean that all turbulent driving takes place at large scales, presumably by external driving forces, and that outflows do not drive enough momentum into the cloud to affect the turbulent cascade. However, Padoan et al. (2009) focused on the power spectrum of integrated intensity maps, while Swift & Welch (2008) specifically examined the power spectra in the line wings.

To investigate the line wing power spectra in Orion A, we visually define the velocity ranges of the line wings where the main cloud emission disappears. Then, we integrate the emission in these velocity ranges and compare their power spectra. Figure 9 shows the 13CO line wing power spectrum toward the L1641N subregion defined in Feddersen et al. (2018), which contains several outflows. In L1641N, we define the blueshifted line wing as 4.3 to 5.6 km s*-1*, the central cloud velocity range as 6.6 to 7.8 km s*-1*, and the redshifted line wing as 10.8 to 12.1 km s*-1*. We find no sign of peaks in the 12CO or 13CO line wing power spectra of any subregion in Orion A.

While the 12CO and 13CO power spectra are featureless power-laws, we do find peaks in the 18CO power spectrum at 30-60″. In Appendix B, we show that these peaks are present in the power spectrum of emission-free channels and are therefore artifacts of the data combination process.

We suggest that the peak in the 13CO line wing power spectra found by Swift & Welch (2008) may be caused by a similar numerical artifact in their data. This peak occurs over angular scales of approximately 40 to 100″. Swift & Welch (2008) combined an interferometric mosaic from the Berkeley-Illinois-Maryland Association (BIMA) array with a map from the Arizona Radio Observatory (ARO) 12m telescope. The spacing between BIMA pointings is 45″ and the FWHM of the ARO 12m beam is 55″. These angular scales are similar to scale of the peak in the power spectrum, suggesting that an artifact introduced in the 13CO data reduction or combination may be responsible for this feature. To test if the peak is real or a data artifact, the power spectrum of emission-free channels should be compared. If the peak appears in this noise power spectrum, it is not intrinsic to the emission and does not indicate a characteristic scale imposed by outflows. We suggest that future studies of molecular cloud power spectra first consider the noise power spectrum before interpreting any deviations from a smooth power-law.

7 Conclusions

This study is one of the first attempts to quantify the connection between feedback and the statistics of turbulent motion in a molecular cloud. Previous studies of feedback in molecular clouds have focused on cataloging and measuring individual features like outflows and shells.

We find spatially correlated emission at relative velocities of 1-3 km s*-1* in the covariance matrix toward the NGC 1977, OMC4, L1641N, and V380 regions. These features resemble those found by B16 in their simulations of stellar winds. It is unclear whether winds are responsible for these features in Orion, or whether they arise from some other mechanism, such as the cloud-cloud collision proposed for L1641-N by Nakamura et al. (2012). We suggest a detailed comparison of the covariance matrix in simulations that incorporate both feedback and a colliding-cloud model to clarify what mechanism produces these features in Orion A.

Contrary to the predictions of the B16 simulations, we do not find a relationship between the slope of the spectral correlation function in Orion A and the momentum injection rate of shells or outflows. However, the uncertainties inherent in both shell and, to a lesser extent, outflow momenta mean we cannot rule out an underlying relationship. A better accounting of the impact of these feedback processes is necessary to fully understand their relationship to the statistics of turbulence.

We find, for the first time, a significant trend between the spectral correlation function and the surface density of young stars. Regions with higher YSO surface density have steeper spectral correlation functions in 12CO. If higher YSO surface density correlates with greater feedback impact, the 12CO SCF slope may be a useful indicator of the importance of feedback in molecular clouds. Feedback is not the only possible underlying variable in this relationship. Any parameter that correlates strongly with stellar density (such as optical depth or excitation temperature) may contribute to the trend we see with the 12CO SCF.

We find no significant trend between the power spectrum in Orion A and the momentum injection rate of shells or outflows, in agreement with B16. We find no evidence for features or breaks in the integrated intensity power spectra or line wing power spectra. Thus, we find no evidence in Orion A for a preferential scale imposed by feedback. Future studies of molecular cloud power spectra should first examine the power spectrum of noise before interpreting any deviations from a smooth power-law.

The statistical study of feedback in molecular clouds is still in its infancy. Future simulations of molecular clouds should compare statistical measures of gas structure and kinematics with the strength of feedback, beyond merely its presence or absence. To aid this, simulations should explore parameter space more fully (see e.g. Yeremi et al., 2014) to disentangle the impacts of feedback from other factors like chemistry and magnetic fields. Future studies should also differentiate between the statistical effects of protostellar outflows, spherical stellar winds, and other types of feedback.

We thank Yoshihiro Tanabe for providing the outflow catalogue. We thank Jens Kauffmann, Jaime Pineda, and the entire CARMA-NRO Orion collaboration for helpful discussion. We thank the referee for thoughtful suggestions that greatly improved the paper. CARMA operations were supported by the California Institute of Technology, the University of California-Berkeley, the University of Illinois at Urbana-Champaign, the University of Maryland College Park, and the University of Chicago. The Nobeyama 45 m telescope is operated by the Nobeyama Radio Observatory, a branch of the National Astronomical Observatory of Japan. This project was supported by the National Science Foundation, award AST-1140063.

Appendix A Statistics of 13CO and C18O

A.1 Covariance Matrix

We show the 13CO and C18O covariance matrices in Figure 10. The off-axis peaks seen in the 12CO covariance matrices are also present in 13CO and C18O. Because these features are prominent in both optically thick and thin tracers, they likely result from real velocity structure in the cloud rather than simply from optical depth effects. B16 computed the covariance matrix of their simulated cube before modeling radiative transfer (Figure 21 in B16). This pre-processed cube does not show covariance peaks, implying that excitation and optical depth effects enhance this feature of their winds. The mapping between these pre-processed cubes and observed optically thin CO lines is unclear. We suggest future models incorporate multiple observable transitions to better compare modeled and observed covariance.

A.2 SCF

We show the 13CO and C18O SCF in Figure 11 and the 13CO and C18O SPS in Figure 12. The 13CO SCF spectra have shapes similar to the 12CO spectra. They follow power-laws up to lags of about 20 pixels ($40\arcsec$ or $0.08$ pc), where they steepen toward larger scales.

The C18O SCF spectra look very different. They show a sharp decrease between the shortest lags of 3 and 5 pixels (6-10$\arcsec$ or 0.01-0.02 pc). Many of the C18O SCF spectra also turn upward at larger scales. An upward sloping SCF means that distantly separated spectra are more similar than close pairs of spectra, an unphysical result. Because the C18O has low signal-to-noise, the upward slopes are more likely to be an artifact of the data reduction process rather than a real signature of the cloud spectra.

A.3 SPS

As in 12CO, the 13CO SPS (Figure 12) closely follows a power-law down to scales of a few beamwidths. The slopes of the 13CO SPS fits range from -3 to -4. Overall, the 13CO slopes are slightly shallower than 12CO, closer to the predicted value of -3 for an optically thick medium and the observational results compiled by Burkhart et al. (2013), although Section 5.2 in Kong et al. (2018) indicates that 13CO is not very optically thick in Orion A, with only 0.6% of pixels having $\tau_{{}^{13}\mathrm{CO}}>1$.

The C18O SPS have significantly shallower slopes than 12CO or 13CO. However, they show a clear peak near 0.1 pc. This peak remains in the power spectrum of emission-free channels (see Appendix B). Because of the low signal-to-noise of C18O, we do not interpret these power spectra.

Appendix B Noise Peak in C18O SPS

The C18O SPS has a peak near 0.1 pc which is not present in either 12CO or 13CO. This feature can also be seen in the integrated C18O delta-variance spectra in Figure 3 of Kong et al. (2018), where they speculate it arises from the low signal-to-noise of C18O coupled with problems in the map cleaning process. Swift & Welch (2008) found a similar feature in the 13CO line-wing power spectrum of the L1551 cloud and attributed this SPS feature to feedback (see Section 6.2). To rule out the influence of feedback in our C18O data, we examine the noise power spectrum. If the feature at 0.1 pc is present in the noise, it is not related to feedback.

To compute the noise power spectrum of C18O, we first sum the emission-free channels between 12 and 16 km s*-1*. Then we compute the spatial power spectrum as described in Section 4.3, including the beam correction and tapering. We show the power spectrum of noise in the L1641N region C18O map in Figure 13. The noise shows the 0.1 pc peak. Because the C18O map has much lower signal-to-noise than 12CO and 13CO, the shape of its noise power spectrum dominates the power spectrum of integrated emission. While the 12CO and 13CO noise power spectra also deviate from power-law, their higher signal-to-noise renders the noise power spectrum insignificant. Before interpreting features in the SPS of any new dataset, we suggest calculating the noise power spectrum.

Appendix C SCF at Large Scales

The SCF shown in Figure 3 only extends to lags of about 45 pixels, or 0.17 pc, to match the range in physical scales covered by the SCF in B16. However, other observational studies fit the SCF spectra at larger scales (e.g. Padoan et al., 2003; Gaches et al., 2015). In Figure 14, we show the 12CO SCF spectra of the Orion A subregions at scales up to 300 pixels (1.2 pc). The SCF steepens at larger scales, in closer agreement with the SCF measured by Padoan et al. (2003) and Gaches et al. (2015). Each subregion SCF flattens and turns upward at different lags. This shape at large scales is a numerical effect caused by the finite size of the subregions and makes it difficult to compare the SCF at large scales between maps of different sizes. If we fit each subregion SCF with a power-law between 30 and 50 pixels (0.12 to 0.2 pc), which avoids the SCF upturn for all subregions, the relationship between SCF slope and feedback impact remains qualitatively the same as shown in Figures 5 through 7.

Appendix D Full-Map Statistics

We show the covariance matrices, SCF, and SPS computed over the entire area covered by CARMA-NRO Orion maps in Figures 15, 16, and 17, respectively.