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

We use the ${}^{12}\mathrm{CO}$(1−0), ${}^{13}\mathrm{CO}$(1−0), and C¹⁸O(1−0) spectral-line maps of Orion A from the Combined Array for Millimeter/Submillimeter Astronomy (CARMA) Nobeyama Radio Observatory (NRO) Orion Survey (Kong et al. 2018). These maps were obtained by combining interferometric images from the CARMA with single-dish maps from the NRO 45 m 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).⁵ The ${}^{12}\mathrm{CO}$ and C¹⁸O maps have a beam full-width at half-maximum (FWHM) of 10'' × 8'' and the ${}^{13}\mathrm{CO}$ map has a beam FWHM of 8'' × 6''. The ${}^{12}\mathrm{CO}$ velocity resolution is 0.25 $\mathrm{km}\,{{\rm{s}}}^{-1}$, while ${}^{13}\mathrm{CO}$ and C¹⁸O have a velocity resolution of 0.22 $\mathrm{km}\,{{\rm{s}}}^{-1}$. For more details on the CARMA–NRO Orion data, see Kong et al. (2018).

To compare the statistics of turbulence with feedback impact, we would ideally like to have a control—a cloud that 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 analog 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 H ii 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 L1641N (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 deg² (5 pc²). We also add a subregion defined by Ungerechts et al. (1997) that 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 ${}^{12}\mathrm{CO}$ peak temperature, and Table 1 defines the extent of each subregion.

Zoom In Zoom Out Reset image size

Table 1.  Subregion Results

Name	R.A. (J2000)	Decl. (J2000)	${\dot{P}}_{\mathrm{shell}}$ a	${\dot{P}}_{\mathrm{out}}$ b	${\mathrm{log}}_{10}{n}_{\mathrm{YSO}}$ c	12CO αSCF	12CO αSPS
Northd	5h37m35s to 5h33m06s	−5°19'26'' to −4°48'19''	17 [6, 41]	6.2 ± 0.5	3.2	−0.096 ± 0.004	−3.07 ± 0.04
Centrald	5h37m35s to 5h33m06s	−5°34'03'' to −5°19'21''	33 [14, 64]	0.13 ± 0.04	3.4	−0.139 ± 0.002	−3.1 ± 0.06
Southd	5h37m36s to 5h33m07s	−6°17'03'' to −5°34'03''	41 [19, 77]	0.7 ± 0.2	2.9	−0.087 ± 0.006	−3.43 ± 0.03
L1641Nd	5h39m01s to 5h34m22s	−7°11'44'' to −6°17'03''	24 [10, 52]	2.2 ± 0.3	2.7	−0.059 ± 0.002	−3.0 ± 0.03
NGC1977e	5h35m47s to 5h34m28s	−4°57'08'' to −4°37'23''	33 [14, 64]	0	3.3	−0.121 ± 0.006	−3.6 ± 0.14
OMC23e	5h36m03s to 5h34m43s	−5°16'37'' to −4°56'37''	28 [4, 85]	19 ± 1	3.4	−0.138 ± 0.003	−3.42 ± 0.07
OMC1f	5h35m24s to 5h35m06s	−5°29'35'' to −5°17'35''	93 [33, 186]	0	4.2	−0.158 ± 0.001	−3.83 ± 0.29
OMC4e	5h35m41s to 5h34m20s	−5°46'56'' to −5°33'00''	162 [78, 279]	0.19 ± 0.05	3.4	−0.111 ± 0.002	−4.46 ± 0.19
OMC5e	5h36m00s to 5h34m41s	−6°15'07'' to −5°50'06''	18 [6, 41]	2.9 ± 0.6	2.9	−0.079 ± 0.002	−4.04 ± 0.1
HH34e	5h36m00s to 5h34m40s	−6°40'07'' to −6°17'00''	23 [13, 40]	0.18 ± 0.04	2.7	−0.082 ± 0.002	−4.18 ± 0.11
L1641Ne	5h37m20s to 5h36m00s	−6°42'15'' to −6°19'02''	36 [15, 71]	5.0 ± 0.9	2.8	−0.077 ± 0.002	−4.04 ± 0.1
V380e	5h37m05s to 5h36m00s	−6°55'00'' to −6°37'00''	22 [4, 75]	7 ± 1	2.7	−0.070 ± 0.002	−4.22 ± 0.12

Notes.

^aMomentum injection rate surface density of shells. Units are 10⁻³ M_⊙ km s⁻¹ yr⁻¹ deg⁻². Values in brackets are summed lower and upper limits from Table 3 of Feddersen et al. (2018). ^bMomentum injection rate surface density of outflows. Units are 10⁻³ M_⊙ km s⁻¹ yr⁻¹ deg⁻². Uncertainties are calculated by adding in quadrature the individual outflow uncertainties in Table 7 of Tanabe et al. (2019). ^cSurface density of YSOs. Units are deg⁻². Because of catalog incompleteness, this is likely an underestimate in the Central and OMC 1 regions. ^dSubregions from Feddersen et al. (2018). ^eSubregions from Davis et al. (2009). ^fSubregions from Ungerechts et al. (1997).

Download table as:  ASCIITypeset image

Feddersen et al. (2018) identified 42 expanding shells in Orion A using the CO maps from the NRO 45 m 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 that 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}}_{\mathrm{shell}}$). The total shell momentum injection rate in each subregion is the sum of ${\dot{P}}_{\mathrm{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}}_{\mathrm{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.

Young stars that are actively accreting material launch collimated bipolar outflows that 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 45 m data has been carried out by Tanabe et al. (2019). 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) 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 OMC 1 are incomplete.

To quantify the impact of outflows on the cloud, we use the momentum injection rates (${\dot{P}}_{\mathrm{out}}$) tabulated in Table 6 of Tanabe et al. (2019). 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). 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). 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) first integrate the ${}^{12}\mathrm{CO}$ emission. ${}^{12}\mathrm{CO}$ 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 ${}^{12}\mathrm{CO}$, other tracers must be observed. For example, Zhang et al. (2016) use ${}^{13}\mathrm{CO}$ and C¹⁸O combined with ${}^{12}\mathrm{CO}$ to measure outflow masses more accurately. Tanabe et al. (2019) only detect a few outflows in both ${}^{12}\mathrm{CO}$ and ${}^{13}\mathrm{CO}$. 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) 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) as well as the variations in the ratio between the integrated intensities of ${}^{12}\mathrm{CO}$ and ${}^{13}\mathrm{CO}$ 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.

As an independent estimate of feedback we consider the young stellar objects (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_YSO to be a proxy for the relative strength of feedback in different regions.

To measure n_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_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 that have ejected outflows or powered shells no longer detectable as coherent structures in the cloud, then n_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 OMC 1 subregions is likely higher relative to the surface density in other subregions. From Figures 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⁻¹ yr⁻¹ deg⁻²) and the YSO surface density (deg⁻²).

PCA is a statistical technique used to reduce the dimensionality of a data set. 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) ≡ 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

$$\begin{eqnarray}&&{S}_{{jk}}=S({v}_{j},{v}_{k})=\displaystyle \frac{1}{n}\displaystyle \sum _{i=1}^{n}{T}_{{ij}}{T}_{{ik}}.\end{eqnarray} \tag{ 1 }$$

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 corresponds 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 $\mathrm{km}\,{{\rm{s}}}^{-1}$ in the covariance matrix of the simulation with winds that 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, 2002b, 2013).

The 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 Δ r (between two pixels) is defined to be

$$\begin{eqnarray}&&\mathrm{SCF}({\boldsymbol{r}},{\rm{\Delta }}{\boldsymbol{r}})=1-\sqrt{\displaystyle \frac{{\displaystyle \sum }_{v}{\left[T({\boldsymbol{r}},v)-T({\boldsymbol{r}}+{\rm{\Delta }}{\boldsymbol{r}},v)\right]}^{2}}{{\displaystyle \sum }_{v}T{\left({\boldsymbol{r}},v\right)}^{2}+{\displaystyle \sum }_{v}T{\left({\boldsymbol{r}}+{\rm{\Delta }}{\boldsymbol{r}},v\right)}^{2}}},\end{eqnarray} \tag{ 2 }$$

where  r is the position of a pixel, Δ r is the lag vector, v is velocity, and T is the temperature (or intensity). The SCF is then averaged over all pixels  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).

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 that 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) that 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 ${}^{13}\mathrm{CO}$ 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 that 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.

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 $\mathrm{km}\,{{\rm{s}}}^{-1}$ that are not present in the simulations without winds. In Figure 2, we show the covariance matrices of the Orion A ${}^{12}\mathrm{CO}$ subregions (the ${}^{13}\mathrm{CO}$ and C¹⁸O 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 $\mathrm{km}\,{{\rm{s}}}^{-1}$ as in the B16 wind simulations. Essentially, these covariance peaks mean there is emission in these subregions that is spatially correlated but separated by a few $\mathrm{km}\,{{\rm{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.

Zoom In Zoom Out Reset image size

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 $\mathrm{km}\,{{\rm{s}}}^{-1}$) dominates emission to the southeast of L1641N while the red component (7–12 $\mathrm{km}\,{{\rm{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 $\mathrm{km}\,{{\rm{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.

Figure 3 shows the ${}^{12}\mathrm{CO}$ 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–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–17 pixels (approximately 0.02–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 ${}^{13}\mathrm{CO}$ maps of Ophiuchus and Perseus. We describe the SCF at scales larger than 0.12 pc in Appendix C.

Zoom In Zoom Out Reset image size

The SCF slopes in the Orion A subregions range from −0.15 to −0.06. Gaches et al. (2015) calculated the ${}^{12}\mathrm{CO}$ 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 ${}^{13}\mathrm{CO}$ 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 data set, 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 comparisons between different data sets 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.

Figure 4 shows the ${}^{12}\mathrm{CO}$ 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).

Zoom In Zoom Out Reset image size

The power-law fits to the SPS in the Orion A subregions have slopes ranging between about −3 and −4 in ${}^{12}\mathrm{CO}$ and ${}^{13}\mathrm{CO}$ with somewhat shallower slopes between about −2.2 and −3.4 in C¹⁸O (see Appendix A for ${}^{13}\mathrm{CO}$ and C¹⁸O 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 ${}^{12}\mathrm{CO}$ and ${}^{13}\mathrm{CO}$ SPS in molecular clouds have shown slopes close to this optically thick limit of −3 (e.g., Stutzki et al. 1998; Padoan et al. 2006; Sun et al. 2006; Pingel et al. 2018), 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' × 50', or 4.4 pc × 4.4 pc at a distance of 300 pc (Ortiz-León et al. 2018; Zucker 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.

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 that include chemistry.

We quantify both SCF and SPS by the slope of their power-law fits that are shown for ${}^{12}\mathrm{CO}$ in Figures 3 and 4. Figures 5–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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 ${}^{12}\mathrm{CO}$ 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 ${}^{12}\mathrm{CO}$ SCF in their simulations. If the true feedback impact is positively correlated with n_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_YSO and ${}^{12}\mathrm{CO}$ 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 ${}^{12}\mathrm{CO}$ emission (e.g., opacity or excitation temperature) may still explain the trend we see between n_YSO and ${}^{12}\mathrm{CO}$ SCF.

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

$$\begin{eqnarray}\begin{array}{rcl}{\alpha }_{\mathrm{SCF}} & = & (-0.060\pm 0.006)\,\mathrm{log}({n}_{\mathrm{YSO}}[{\deg }^{-2}])\\ & & +\ (0.090\pm 0.020),\end{array}\end{eqnarray} \tag{ 3 }$$

where α_SCF is the ${}^{12}\mathrm{CO}$ SCF slope. This fit has a correlation coefficient of r² = 0.84. If we exclude the OMC 1 subregion, which has an order of magnitude higher n_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_YSO, which is also consistent with B16. In ${}^{13}\mathrm{CO}$ and C¹⁸O, neither SCF nor SPS appear to be correlated with n_YSO. It is unclear why the ${}^{12}\mathrm{CO}$ SCF slope is correlated with n_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 Figure 5 through 7. In particular, the anticorrelation between ${}^{12}\mathrm{CO}$ SCF slope and n_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.

Zoom In Zoom Out Reset image size

The power spectra in Figure 4 are calculated using the ${}^{12}\mathrm{CO}$ 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 ${}^{13}\mathrm{CO}$ 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 ${}^{13}\mathrm{CO}$ 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–5.6 $\mathrm{km}\,{{\rm{s}}}^{-1}$, the central cloud velocity range as 6.6–7.8 $\mathrm{km}\,{{\rm{s}}}^{-1}$, and the redshifted line wing as 10.8–12.1 $\mathrm{km}\,{{\rm{s}}}^{-1}$. We find no sign of peaks in the ${}^{12}\mathrm{CO}$ or ${}^{13}\mathrm{CO}$ line wing power spectra of any subregion in Orion A.

Zoom In Zoom Out Reset image size

While the ${}^{12}\mathrm{CO}$ and ${}^{13}\mathrm{CO}$ power spectra are featureless power laws, we do find peaks in the ${}^{18}\mathrm{CO}$ 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 ${}^{13}\mathrm{CO}$ 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''–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) 12 m telescope. The spacing between BIMA pointings is 45'' and the FWHM of the ARO 12 m beam is 55''. These angular scales are similar to the scale of the peak in the power spectrum, suggesting that an artifact introduced in the ${}^{13}\mathrm{CO}$ 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.

We show the ${}^{13}\mathrm{CO}$ and C¹⁸O covariance matrices in Figure 10. The off-axis peaks seen in the ${}^{12}\mathrm{CO}$ covariance matrices are also present in ${}^{13}\mathrm{CO}$ and C¹⁸O. 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 preprocessed cube does not show covariance peaks, implying that excitation and optical depth effects enhance this feature of their winds. The mapping between these preprocessed cubes and observed optically thin CO lines is unclear. We suggest future models incorporate multiple observable transitions to better compare modeled and observed covariance.

Zoom In Zoom Out Reset image size

We show the ${}^{13}\mathrm{CO}$ and C¹⁸O SCF in Figure 11 and the ${}^{13}\mathrm{CO}$ and C¹⁸O SPS in Figure 12. The ${}^{13}\mathrm{CO}$ SCF spectra have shapes similar to the ${}^{12}\mathrm{CO}$ spectra. They follow power laws up to lags of about 20 pixels (40'' or 0.08 pc), where they steepen toward larger scales.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The C¹⁸O SCF spectra look very different. They show a sharp decrease between the shortest lags of 3 and 5 pixels (6''–10'' or 0.01–0.02 pc). Many of the C¹⁸O 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 C¹⁸O 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.

As in ${}^{12}\mathrm{CO}$, the ${}^{13}\mathrm{CO}$ SPS (Figure 12) closely follows a power law down to scales of a few beamwidths. The slopes of the ${}^{13}\mathrm{CO}$ SPS fits range from −3 to −4. Overall, the ${}^{13}\mathrm{CO}$ slopes are slightly shallower than ${}^{12}\mathrm{CO}$, 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 ${}^{13}\mathrm{CO}$ is not very optically thick in Orion A, with only 0.6% of pixels having ${\tau }_{{}^{13}\mathrm{CO}}\gt 1$.

The C¹⁸O SPS have significantly shallower slopes than ${}^{12}\mathrm{CO}$ or ${}^{13}\mathrm{CO}$. 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 C¹⁸O, we do not interpret these power spectra.