# A Census of Protostellar Outflows in Nearby Molecular Clouds

Duo Xu  Department of Astronomy, The University of Texas at Austin, Austin, TX
78712, USA;  Department of Astronomy, University of Virginia, Charlottesville,
VA 22904-4235, USA<EMAIL_ADDRESS>Stella S. R. Offner  Department of
Astronomy, The University of Texas at Austin, Austin, TX 78712, USA;
<EMAIL_ADDRESS>Robert Gutermuth  Department of Astronomy,
University of Massachusetts, Amherst, MA 01003, USA;  Shuo Kong  Steward
Observatory, University of Arizona, Tucson, AZ 85719, USA;  Department of
Astronomy, Yale University, New Haven, CT 06511, USA;  Hector G. Arce
Department of Astronomy, Yale University, New Haven, CT 06511, USA;

###### Abstract

We adopt the deep learning method casi-3d (Convolutional Approach to Structure
Identification-3D) to systemically identify protostellar outflows in 12CO and
13CO observations of the nearby molecular clouds, Ophiuchus, Taurus, Perseus
and Orion. The total outflow masses are 267 $M_{\mathord{\odot}}$, 795
$M_{\mathord{\odot}}$, 1305 $M_{\mathord{\odot}}$ and 6332
$M_{\mathord{\odot}}$ for Ophiuchus, Taurus, Perseus and Orion, respectively.
We show the outflow mass in each cloud is linearly proportional to the total
number of young stellar objects. The estimated total 3D deprojected outflow
energies are $9\times 10^{45}$ ergs, $6\times 10^{46}$ ergs, $1.2\times
10^{47}$ ergs and $6\times 10^{47}$ ergs for Ophiuchus, Taurus, Perseus and
Orion, respectively. The energy associated with outflows is sufficient to
offset turbulent dissipation at the current epoch for all four clouds. All
clouds also exhibit a break point in the spatial power spectrum of the outflow
prediction map, which likely corresponds to the typical outflow mass and
energy injection scale.

ISM: outflows – ISM: clouds – methods: data analysis – stars: formation

##  1 Introduction

Protostars launch energetic collimated bipolar outflows during the star
formation process. This outflow gas entrains and accelerates the surrounding
material to higher velocities than that of the ambient cloud gas, injecting a
substantial amount of energy into the host molecular cloud (Frank et al.,
2014; Bally, 2016). Numerical simulations suggest protostellar outflows not
only reduce the core and stellar mass (e.g., Hansen et al., 2012; Offner &
Chaban, 2017) but also depress the star formation rate and star formation
efficiency in molecular clouds (e.g., Federrath et al., 2014; Federrath, 2015;
Cunningham et al., 2018). Meanwhile, protostellar outflows likely control the
shape of the stellar initial mass function (Cunningham et al., 2018; Guszejnov
et al., 2021). Theoretical and numerical studies indicate that protostellar
outflows are highly efficient at driving turbulent motions (Matzner, 2007;
Nakamura & Li, 2007; Moraghan et al., 2013). In small to intermediate size
clouds, energy injected by protostellar outflows can compensate for rapid
turbulent dissipation over at least several dynamical time scales (Cunningham
et al., 2006; Li & Nakamura, 2006; Matzner, 2007; Carroll et al., 2009).
Consequently, protostellar outflows may stall the local gravitational collapse
of molecular clouds, thereby extending the cloud lifetime (Wang et al., 2010).
The role of protostellar feedback on the fate of molecular clouds is still
under debate (Bally, 2016).

Unfortunately, the mechanism by which protostellar outflows convert their
kinetic energy into turbulent energy remains poorly understood observationally
(Frank et al., 2014). In practice, astronomers usually evaluate outflow impact
by comparing the total kinetic energy from protostellar outflows with the
turbulent energy of their host cloud. For example, Arce et al. (2010)
identified 60 outflow candidates in the Perseus molecular cloud and concluded
the total outflow energy is sufficient to replenish the dissipation of
turbulence. A similar result was also found in Taurus (Li et al., 2015), in
Ophiuchus (Nakamura et al., 2011), and in Orion (Feddersen et al., 2020).

Accurately determining the impact of protostellar outflows, especially
quantifying their effect on the cloud energy budget, requires a complete
census of outflows. However, it is challenging to identify most of the
outflows that are deeply embedded in dense clouds (Arce et al., 2010; Dunham
et al., 2014a). The morphology of these embedded outflows are not as distinct
as isolated ones in low velcity channels and the outflowing gas is difficult
to separate from the ambient cloud. Protostars are often clustered and their
outflows interact, making visual identification challenging. One solution is
to calculate the outflow gas whose velocity is significantly above the cloud
velocity dispersion and extrapolate the total mass (Arce et al., 2010; Li et
al., 2015; Feddersen et al., 2020). However, this leads to an order of
magnitude uncertainty (Arce et al., 2010).

Recent developments in machine learning approaches have enabled automated
detection, which can separate outflows from cloud emission and conduct
systematic identification of outflow features (Zhang et al., 2020; Xu et al.,
2020b). Zhang et al. (2020) employed Support Vector Machines (SVM) to identify
molecular outflows in a dark cloud complex from molecular line emission.
Although SVM performs robustly in classification tasks, it requires
preprocessing of the raw data cubes, i.e., manually extracting feature vectors
that represent the raw data. The choice of input features are determined
subjectively by visual inspection, and extracting these features from data
cubes is non-trivial. Moreover, the manually extracted features discard part
of the information embedded in the raw spectral cubes, which introduces
uncertainty and affects the performance of classification. Convolutional
neural networks (CNNs) are a powerful new approach being applied to identify
structures or objects in astronomical data, such as exoplanets (Shallue &
Vanderburg, 2018), stellar feedback bubbles (Van Oort et al., 2019; Xu et al.,
2020a) and protostellar outflows (Xu et al., 2020b). Given a labeled training
set of images or spectral cube data, CNNs can be applied to efficiently
identify features in large surveys. Van Oort et al. (2019) developed a
Convolutional Approach to Shell/Structure Identification, casi, to identify
stellar feedback bubbles in 2D density slices and 12CO integrated intensity
maps. Xu et al. (2020a) and Xu et al. (2020b) successfully extended casi to
casi-3d, which is able to identify stellar feedback bubbles and protostellar
outflows in position-position-velocity (PPV) molecular line spectra cubes. Xu
et al. (2020b) applied casi-3d to 12CO observations of the Perseus molecular
cloud and identified all 60 previously visually identified outflows.
Additionally, casi-3d found 20 new high-confidence outflows. Apart from
structure identification, casi-3d successfully predicts hidden information in
the data cube, such as the fraction of mass associated with feedback, which
provides a more accurate feedback mass estimation. Xu et al. (2020a) showed
that the actual mass associated with stellar feedback in Taurus is an order of
magnitude smaller than the previous visual estimates. These results illustrate
the capability of CNNs to identify complex structures and infer embedded
information.

In this paper, we adopt the deep learning method casi-3d to systemically
identify protostellar outflows in the nearby molecular clouds, Ophiuchus,
Taurus, Perseus and Orion. We describe casi-3d and the CO observations of
these nearby molecular clouds in Section 2. In Section 3, we present the
performance of our CNN models in identifying protostellar outflows in the
observational data, calculate the physical properties of outflows and discuss
their impact on the clouds. We discuss the broader impact of outflows and
compare them with stellar wind driven bubbles in Section 4 and summarize our
results and conclusions in Section 5.

##  2 Data and Method

###  2.1 COMPLETE Survey

The 12CO J=1-0 (115.271 GHz) and 13CO J=1-0 (110.201 GHz) lines were observed
simultaneously in surveys of Ophiuchus, Taurus and Perseus between 2002 and
2005 using the 13.7 m Five College Radio Astronomy Observatory (FCRAO)
Telescope (Ridge et al., 2006; Narayanan et al., 2008). The main beam of the
antenna pattern has a FWHM of 45″ for 12CO and 47″ for 13CO. The data are
obtained on the fly (OTF), but they are resampled onto a uniform 23″ grid
(Ridge et al., 2006).

The Ophiuchus data has root-mean-square (RMS) antenna temperatures of 0.98 K
and 0.33 K for 12CO and 13CO, respectively. We resampled the spectra with a
lower velocity resolution of 0.125 km/s to ensure a uniform velocity
resolution for all regions and to match the resolution of the training set. If
the velocity resolution is too low ($>0.3$ km s-1), the structure generated by
stellar feedback is not distinguishable across multiple velocity channels. On
the other hand, if the velocity resolution is too high ($<0.1$ km s-1), the
training data will hit the limitation of GPU memory. The noise levels for the
new Ophiuchus 12CO and 13CO spectra are reduced by a factor of square root of
2 to 0.69 K and 0.23 K, respectively. The final Ophiuchus data cube has a
velocity range between -0.8 and 7.5 km s-1 with 67 channels.

The Taurus data has a RMS antenna temperature of 0.28 K for 12CO and 0.125 K
for 13CO. There are 80 and 76 channels with 0.26 and 0.27 km s-1 spacing for
12CO and 13CO, respectively. The velocity range of the Taurus data spans -5.1
to 14.9 km s-1.

The Perseus 12CO and 13CO data have a RMS antenna temperature of 0.25 K and
0.2 K, respectively. We resampled the spectra with a lower velocity resolution
of 0.125 km/s to ensure a uniform velocity resolution for all regions and to
match the resolution of the training set. The noise levels for the new 12CO
and 13CO spectra is reduced by a factor of square root of 2 to 0.17 K and 0.14
K, respectively. The final Perseus data cube has a velocity range between -2.0
and 15.0 km s-1 with 137 channels.

###  2.2 NRO45 Orion Survey

The 12CO J=1-0 and 13CO J=1-0 observations of Orion A were carried out from
2007 to 2017 by the Nobeyama Radio Observatory 45 m telescope (NRO 45m), using
two different receivers (Shimajiri et al., 2015a, b; Ishii et al., 2019;
Nakamura et al., 2019). The two maps were calibrated to the same intensity
scale and combined on a common grid with a pixel scale of 7″.5, which
corresponds to an effective angular resolution of 22″. Ishii et al. (2019)
smoothed the combined map to a velocity resolution of 0.22 km/s, and further
converted the intensity to the main beam temperature scale that is benchmarked
with data from FCRAO and IRAM 30m. The final sensitivity for the 12CO map is
0.35 K, while the sensitivity for the 13CO map is 0.40 K. The velocity range
of the Orion data spans -2.9 to 19.3 km s-1. More detailed description about
the data can be found in Kong et al. (2018) and Ishii et al. (2019).

###  2.3 YSO Catalogue

To validate our outflow identifications, we compare our feedback maps with the
observed distributions of YSOs. We use the YSO catalog for Ophiuchus and
Perseus from SESNA (Spitzer Extended Solar Neighborhood Archive, Gutermuth et
al., in prep) used by Pokhrel et al. (2020). SESNA uses an updated
implementation of the data treatment, source catalog construction, and YSO
identification and classification processes on Spitzer surveys (Gutermuth et
al., 2009). SESNA classified YSOs into four groups: deeply embedded
protostars, Class I YSOs, Class II YSOs, and transition disks. For further
analysis, we combine the former two groups as “Younger YSOs” and merge the
latter two groups as “Older YSOs”.

We adopt the YSO catalog from Rebull et al. (2010) for Taurus, which contains
215 previously identified YSOs and 148 newly identified YSOs. Rebull et al.
(2010) adopted observations from IRAC and MIPS between 2005 and 2007. Rebull
et al. (2010) applied color-magnitude criteria to select the YSO candidates
and classified these YSO candidates into four different classes: Class I, II,
III and Flat. For further analysis, we combine YSOs in Class I and Flat into a
single category of “Younger YSOs” and merge YSOs in Class II and III into a
single category of “Older YSOs”.

We adopt the YSO catalog from Megeath et al. (2012) for Orion, which contains
a total of 3479 dusty YSOs from IRAC and MIPS observations. Megeath et al.
(2012) classified these YSOs into three categories, which includes 2991 young
stars with disks, 428 protostars, and 50 faint candidate protostars. For
further analysis, we designate protostars as “Younger YSOs” and young stars
with disks as “Older YSOs”.

YSOs in the group of “Younger YSOs” have an average age smaller than 1 Myr,
while YSOs in the group of “Older YSOs” have an average age around 1-3 Myrs
(Dunham et al., 2014b). However, it is worth noting that truly younger YSOs
and older YSOs may be misclassified based on their orientation with respect to
the line of sight (Offner et al., 2012), which means the actual age of an
individual YSO in the group of “Younger YSOs” is not necessarily smaller than
that of an object in the group of “Older YSOs”.

Finally, we note that not all the surveys completely cover the area of
observed molecular emission. We focus our analysis on the areas with YSO
coverage and indicate the boundaries of the surveys on the maps.

###  2.4 Method

####  2.4.1 casi-3d

We adopt the previously trained casi-3d model from Xu et al. (2020b) to
identify protostellar outflows in 12CO. casi-3d is an encoder-decoder based
convolutional neural network combining both residual networks (He et al.,
2016) and a “U-net” (Ronneberger et al., 2015). casi-3d is trained on
synthetic 12CO data, which models forming stars that launch protostellar
outflows. We adopt the same training set and ranges of magnetohydrodynamic
(MHD) model properties as described in Xu et al. (2020b, Table 1). We adopt
three different 12CO abundances and two different cloud kinetic temperatures
when conducting synthetic observations as described in Xu et al. (2020b, Table
2). casi-3d takes 12CO data cubes as input and predicts the position of
outflows on the voxel level. Xu et al. (2020b) trained two casi-3d models to
identify protostellar outflows. One model, model ME1, is trained to predict
the 12CO emission that is associated with outflows. The other model, model MF,
is trained to predict the fraction of the mass that comes from stellar
feedback. A more detailed description of how we generate training data for
these two models can be found in Xu et al. (2020b).

####  2.4.2 Data Preprocessing

Before we apply our casi-3d models to the observational data, we apply the
same preprocessing steps that we adopt for the training data. Due to the
relatively large dynamic range of the 12CO emission, we take the logarithm of
the emission and then normalize the values by subtracting the mean and
dividing by the standard deviation of the full map. Figure 1 illustrates the
cumulative distribution of normalized emission for the four clouds. Although
the 12CO emission varies significantly between clouds, after preprocessing,
the normalized emission of the four clouds is similar. This allows casi-3d to
perform stably across a variety of conditions.

Xu et al. (2020b) examined the performance of casi-3d on conditions that are
not included in the training set, such as different kinetic temperatures, 12CO
abundances, beam sizes and noise levels. casi-3d is able to successfully
identify outflows across a variety of physical and observational conditions.
See Appendix C in Xu et al. (2020b) for more details. These tests suggest that
the casi-3d identifications are relatively robust across the range of
conditions expected to be characteristic of the observed clouds.

Figure 1: Cumulative distribution function (CDF) of the normalized emission
for the four clouds.

After preprocessing the observational data, we adopt the same strategy as Xu
et al. (2020a) to carry out the full 12CO map prediction. We crop the full
12CO map into a stack of three-dimensional chunks. Each chunk has a size of
$64\times 64\times 32$ (position-position-velocity). To reduce the bias due to
the position of outflows, each chunk has at least 84% volume overlap with
adjacent chunks. To construct the full-map prediction, we combine the
predictions for the individual chunks and adopt the maximum predicted value
for each voxel.

####  2.4.3 Mass Calculations

We follow the same strategy as Arce et al. (2010) to calculate the outflow
mass by combining both 12CO and 13CO data. If there is distinct 13CO emission
at the corresponding position, we use 13CO to calculate the outflow mass. We
assume the 13CO emission line is optically thin and the 12CO emission line is
optical thick. We adopt an excitation temperature of 25 K or the 12CO peak
temperature, whichever is higher, to calculate the mass (Arce et al., 2010;
Narayanan et al., 2012; Li et al., 2015; Feddersen et al., 2020). If there is
no distinct 13CO emission at the corresponding position, we assume the 12CO
emission line is optically thin to derive the mass. Under the assumption of
local thermodynamic equillibrium (LTE), the mass estimation scales linearly
with the excitation temperature. From previous feedback mass estimates, the
choice of excitation temperature ranges from 10 K to 50 K. This could
potentially introduce a factor of two uncertainty in the mass estimation. We
take 62 as the abundance ratio between 12CO and 13CO and 10-4 as the abundance
ratio between 12CO and H2 (Arce et al., 2010; Feddersen et al., 2020). We
verify that the 13CO emission is generally optically thin for all four clouds.
For example, Orion has the largest cloud mass and highest column density among
the four clouds. The optical depth of 13CO in Orion is less than 1 for 99.4%
of the pixels (see also Kong et al., 2018). Only 0.6% of the pixels have
$\tau_{13}>$ 1, and these are mostly in the OMC-2/3 and L1641-N regions. The
maximum value of $\tau_{13}$ is 7.8.

##  3 Results

###  3.1 Outflows Identified in the Full Map

In this section, we present the prediction by models ME1 and MF on the four
star-forming clouds (Ophiuchus, Taurus, Perseus and Orion). Most molecular
clouds have a global velocity gradient, which indicates there is no one unique
central velocity of a cloud. To better visualize the blue-shifted and red-
shifted lobes of outflows, we remove the large scale gradient to make local
outflow motions clearer, i.e., shifting the velocity zero point to the flux-
weighted central velocity of the spectrum for each pixel. We calculate the
flux-weighted central velocity from the 1st moment maps. In order to reduce
the velocity fluctuations on small scales, we apply a Gaussian kernel with a
FWHM that corresponds to 0.1 pc to convolve the 1st moment maps. For each
pixel, we show the integrated prediction over the channels that have absolute
velocities greater than the central velocity by the specified thresholds.

####  3.1.1 Ophiuchus

Figure 2 shows ME1 and MF model predictions for Ophiuchus. Both models
identify a large amount of outflow activity in the region of the large,
central star cluster. In contrast, the models identify much less outflow
activity just north of the cluster, which has only a few YSOs.

Figure 3 shows the predicted outflow activity by models ME1 and MF toward the
young star cluster region Lynds 1688 in Ophiuchus. The L1688 region has over
120 young stars and a large number of interacting outflows. Model ME1, which
performs similarly to human visual identification, identifies almost
everything as feedback in this region. In contrast, Model MF provides an
approach to disentangle the outflow emission from the ambient cloud emission.
Quantitively, model ME1 predicts 81% of the mass in this region as outflows.
However, model MF predicts that only 12% of the total mass is associated with
outflows. Figure 3 shows that both models ME1 and MF are able to identify
coherent high velocity structures in the position-velocity diagram.

It is challenging to separate outflow emission from the cloud emission in 12CO
(1-0) “by eye,” since the cloud emission dominates the 12CO (1-0) emission.
Consequently, even an expert astronomer can only separate the high-velocity
component from the host cloud emission but cannot identify the outflow
morphology near the rest frame velocity of the cloud.

Figure 2: Intensity of 12CO (1-0) integrated over all velocity channels for
Ophiuchus, overlaid with the model ME1 prediction (upper panel in red and blue
contours) and with the model MF prediction (lower panel in red and blue
contours). Red contours indicate the integrated prediction over the channels
that have velocities greater than $V_{\rm cen}+0.4$ km s-1 for model ME1 and
$V_{\rm cen}+0.2$ km s-1 for model MF. Blue contours indicate the integrated
prediction over the channels that have velocities smaller than $V_{\rm
cen}-0.4$ km s-1 for model ME1 and $V_{\rm cen}-0.2$ km s-1 for model MF. “Y”
and “O” indicates the location of YSOs, as described in Section 2.3. The
contours start at 25th percentile of the sorted pixel values and end at 99.7th
percentile of the sorted pixel values of the data, with 6 levels evenly
spaced. It is worth noting that the absolute values of the contour levels for
models ME1 and MF are different. The contour levels for the model ME1
prediction start at 2.8 K$\cdot$km/s, for the blue-shifted lobes and 2.2
K$\cdot$km/s, for the red-shifted lobes, and end at 19.3 K$\cdot$km/s, and
15.3 K$\cdot$km/s, for the blue- and red-shifted lobes, respectively. The
contour levels for the model MF prediction start at 0.80 K$\cdot$km/s, (blue)
and 0.96 K$\cdot$km/s, (red), and end at 7.2 K$\cdot$km/s, (blue) and 8.6
K$\cdot$km/s, (red). The method of plotting contours is the same for Figure 4,
5, 7 and 8. Figure 3: Position-velocity diagram of 12CO emission toward the
Ophiuchus L1688 region. Left panel: integrated intensity of 12CO over the the
full velocity range (from -0.8 km/s to 7.5 km/s) overlaid with the model ME1
and MF predictions in white contours. Letters “Y” and “O” mark YSO positions,
as described in Section 2.3. The purple line illustrates the cut direction of
the position-velocity diagram. Middle and right panel: position-velocity
diagram of 12CO emission overlaid with the model ME1 and MF predictions in
white contours.

####  3.1.2 Taurus

Figure 4 and 5 show the ME1 and MF model predictions for Taurus. The
structures in the outflow prediction maps are more discrete than those in
Ophiuchus. One possible reason is that YSOs in Taurus are more sparsely
distributed than those in Ophiuchus. Although the total number of YSOs in
Taurus is slightly larger, they are less clustered. In Ophiuchus, the outflows
driven by YSOs are more likely to interact and overlap with each other. While
in Taurus, the outflows are more isolated.

Figure 6 shows the outflow activity predicted by models ME1 and MF toward a
previously identified outflow, TMO_06 (Li et al., 2015), in Taurus. Several
young YSOs are located around this outflow. Both models are able to identify
the blueshifted and redshifted lobes of this outflow. However, the model MF
prediction is more extended in the position-velocity diagram. In high velocity
channels, the 12CO emission produced by outflows is usually faint while the
fraction of outflow mass is high. Model ME1 only identifies the location of
outflows and gives stronger emission a higher weight. In this subregion, model
ME1 predicts that 43% of the mass in this region is associated with outflows.
However, model MF predicts that only 8% of the total mass is outflows gas.

Figure 4: Intensity of 12CO (1-0) integrated over all velocity channels for
Taurus, overlaid with the model ME1 prediction (red and blue contours).
Letters “Y” and “O” mark YSO positions, as described in Section 2.3. The grey
line encloses the Spitzer coverage, where YSOs are identified (Rebull et al.,
2010). The contours start at the 25th percentile of the sorted pixel values
(0.47 K$\cdot$km/s, for blue contours and 0.49 K$\cdot$km/s, for red contours)
and end at the 99.7th percentile of the sorted pixel values of the data (2.9
K$\cdot$km/s, for blue contours and 3.4 for red contours), with 6 levels
evenly spaced. Note that the absolute values of the contour levels for models
ME1 and MF are different. Figure 5: Intensity of 12CO (1-0) integrated over
all velocity channels for Taurus, overlaid with the model MF prediction (red
and blue contours). Letters “Y” and “O” mark YSO positions, as described in
Section 2.3. The grey line encloses the Spitzer coverage, where YSOs are
identified (Rebull et al., 2010). The contours start at the 25th percentile of
the sorted pixel values (0.27 K$\cdot$km/s, for blue contours and 0.27
K$\cdot$km/s, for red contours) and end at the 99.7th percentile of the sorted
pixel values of the data (1.7 K$\cdot$km/s, for blue contours and 1.5 for red
contours), with 6 levels evenly spaced. Figure 6: Position-velocity diagram of
12CO emission toward a previously identified outflow, TMO_06 (Li et al.,
2015), in Taurus. Left panel: integrated intensity of 12CO over the the full
velocity range (from -1.5 km/s to 13.4 km/s) overlaid with the model ME1 and
MF predictions in white contours. Letters “Y” and “O” mark YSO positions, as
described in Section 2.3. The purple line illustrates the cut direction of the
position-velocity diagram. Middle and right panel: position-velocity diagram
of 12CO emission overlaid with the model ME1 and MF predictions in white
contours.

####  3.1.3 Perseus

Xu et al. (2020b) followed up previous visually identified outflow targets and
validated the model performance. Here we extend that analysis by carrying out
a “blind search” where we analyze the full cloud. Figure 7 shows the full map
predictions for both models. Both model ME1 and MF predictions are more
concentrated towards star clusters. The model ME1 prediction is more spatially
extended than that by model MF. This might be caused by a strong emission
region with a low outflow mass fraction, where model ME1 counts the entire
voxel as feedback while model MF recognizes that much of this emission is
cloud contamination and excludes it.

Figure 7: Intensity of 12CO (1-0) integrated over all velocity channels for
Perseus, overlaid with the model ME1 prediction (upper panel in red and blue
contours) and with the model MF prediction (lower panel in red and blue
contours). Letters “Y” and “O” mark YSO positions, as described in Section
2.3. The grey line encloses the Spitzer coverage, where YSOs are identified
(Gutermuth et al., in prep). The contour levels for the model ME1 prediction
start at 2.5 K$\cdot$km/s (blue) and 2.1 K$\cdot$km/s (red), and end at 17
K$\cdot$km/s (blue) and 15 K$\cdot$km/s (red). The contour levels for the
model MF prediction start at 0.96 K$\cdot$km/s (blue) and 0.90 K$\cdot$km/s
(red), and end at 8.6 K$\cdot$km/s (blue) and 8.0 K$\cdot$km/s (red).

####  3.1.4 Orion

Figure 8 shows the model ME1 and MF predictions for Orion. Due to the large
number of old evolutionary stage YSOs in Orion, we only show the younger YSOs
to reduce confusion. We show the position of all Orion YSOs in Figure 22 in
Appendix C. The prediction by both models covers most of the emission in
Orion, which indicates outflows exist everywhere. This is consistent with the
extremely dense distribution of YSOs in Orion.

Figure 9 shows the predicted outflow activity by models ME1 and MF toward a
previously identified outflow, Outflow No. 7 (Tanabe et al., 2019), in Orion.
Both models similarly highlight the coherent high-velocity structures. There
are several young YSOs located around this outflow, which could be the driving
source. For reference, we illustrate the high-velocity components in blue and
red contours, which are considered to be outflows by Tanabe et al. (2019). The
mass of these high-velocity components is 268 $M_{\mathord{\odot}}$. Model MF
predicts the outflow mass to be 308 $M_{\mathord{\odot}}$, which is similar to
previous estimates. In this subregion, model ME1 predicts 61% of the mass in
this region is associated with outflows. However, model MF predicts that only
13% of the total mass is outflow gas.

Figure 8: Intensity of 12CO integrated over all velocity channels for Orion,
overlaid with the model ME1 prediction (left panel in red and blue contours)
and with the model MF prediction (right panel in red and blue contours).
Letters “Y” mark young YSO positions, as described in Section 2.3. The contour
levels for the model ME1 prediction start at 2.1 K$\cdot$km/s (blue) and 2.6
K$\cdot$km/s (red), and end at 14 K$\cdot$km/s (blue) and 18 K$\cdot$km/s
(red). The contours levels for the model MF prediction start at 0.48
K$\cdot$km/s (blue) and 0.55 K$\cdot$km/s (red), and end at 4.3 K$\cdot$km/s
(blue) and 4.9 K$\cdot$km/s (red). Figure 9: Position-velocity diagram of 12CO
emission toward a previously identified outflow, Outflow No. 7 (Tanabe et al.,
2019), in Orion. Left panel: integrated intensity of 12CO over the the full
velocity range (from 2.3 km/s to 17.1 km/s) overlaid with the model ME1 and MF
predictions in white contours. The blue and red contours indicate the high
velocity components, integrated over the raw emission with velocity ranges of
2.3-9.5 km s-1and 12.3-17.0 km s-1, respectively. Letters “Y” and “O” mark YSO
positions, as described in Section 2.3. The purple line illustrates the cut
direction of the position-velocity diagram. Middle and right panel: position-
velocity diagram of 12CO emission overlaid with the model ME1 and MF
predictions in white contours.

###  3.2 Physical Properties of Outflows

In this section, we study the physical properties of the outflows we
identified and analyze the correlation between the physical quantities of
outflows and the number of YSOs in the clouds. This analysis help us evaluate
the robustness of the outflow identification by casi-3d. If the outflows are
accurately identified by casi-3d, we should expect a linear correlation
between the outflow mass and the number of YSOs. This is exactly what we find
in this section and Section 3.3. Because our analysis aims to correlate the
physical properties of outflows with the number of YSOs, we only calculate
mass, energy and momentum of outflows in the source catalog covered area.

We follow the method in Section 2.4.3 to calculate the mass of outflows. We
adopt the model MF estimates for the fiducial values. Figure 10 shows the mass
estimates for the outflows in the four molecular clouds. When there are more
YSOs in the region, the predicted outflow mass is higher. In Ophiuchus, on
average, each YSO (including both young and old) contributes 0.43
$M_{\mathord{\odot}}$ outflow mass to the host cloud. While the values are 1.7
$M_{\mathord{\odot}}$, 2.8 $M_{\mathord{\odot}}$, and 3.4
$M_{\mathord{\odot}}$ for Taurus, Perseus and Orion, respectively. On average,
in each of the four clouds, the total mass associated with feedback is about
10% of the mass of the host cloud within the area studied.

We define the 1D (line-of-sight, LOS) momentum as the sum of the gas mass in
each channel multiplied by the channel velocity, where we have shifted the
mean cloud velocity to zero. To better quantify LOS momentum, we subtract the
central velocity along each sightline to reduce the effect of large velocity
gradients across the entire cloud, as described in Section 3.1. Subtracting
the velocity gradient provides a better estimate of LOS momentum, which
reduces the impact of large-scale motions on our estimates. For comparison, we
also calculate the LOS moment without subtracting the velocity gradient. The
LOS momenta of the outflows and the host clouds are calculated by equations:

$P(VGsub)=\sum_{i,j,k}M_{{\rm CO},i,j,k}(v_{i,j,k}-\overline{v_{i,j}}),$ (1)

$P(VGnonsub)=\sum_{i,j,k}M_{{\rm CO},i,j,k}(v_{i,j,k}-\overline{v_{\rm
global}}).$ (2)

Similarly, we define the LOS kinetic energy as

$E(VGsub)=\sum_{i,j,k}\frac{1}{2}M_{{\rm
CO},i,j,k}(v_{i,j,k}-\overline{v_{i,j}})^{2},$ (3)

$E(VGnonsub)=\sum_{i,j,k}\frac{1}{2}M_{{\rm
CO},i,j,k}(v_{i,j,k}-\overline{v_{\rm global}})^{2}.$ (4)

Figure 11 shows the LOS outflow momentum estimates for the four molecular
clouds. We list two sets of momentum estimates, one that subtracts the
velocity gradient and a second that does not subtract the velocity gradient.
We adopt the momentum estimate by model MF, which subtracts the velocity
gradients as the fiducial estimate. The trend between the momentum driven by
outflows and the number of YSOs is similar to that of the mass estimates for
four clouds. In Ophiuchus, the 1D momentum driven by outflows is around 10% of
the total 1D momentum of the host cloud. In Taurus, this momentum ratio is 8%.
In Perseus and Orion, the 1D momentum driven by outflows is around 16% of the
total 1D momentum, which indicates there are likely more high velocity
structures identified as outflows. This is consistent with the star formation
history of these clouds. Taurus has a relatively low star formation rate and
mainly forms low mass stars. On the other hand, Ophiuchus and Orion are home
to clusters forming higher-mass stars and have more active star formation.
Perseus also contains some intermediate mass stars, including a couple of
B-type stars (Arce et al., 2011).

Ophiuchus, Perseus and Orion have larger velocity gradients than Taurus, which
indicates these clouds contain more complex cloud kinematic structure. These
velocity gradients are likely caused by large scale processes, such as gas
accretion (Klessen & Hennebelle, 2010) or converging atomic gas flows caused
by spiral density waves (Wolfire et al., 1995; Heitsch et al., 2006). Given
the magnitude of these bulk motions, we expect that the momentum estimate that
includes the velocity gradient overestimates the contribution of outflows,
while the momentum that subtracts the velocity gradient provides a better
estimate.

Figure 12 shows the LOS outflow kinetic energy for the four molecular clouds.
We adopt the kinetic energy estimates by model MF that subtracts the velocity
gradient for the comparisons. In Ophiuchus, the LOS kinetic energy injected by
outflows is around 6% of the total LOS kinetic energy of the host cloud. In
Taurus, this energy ratio is 4%. In Perseus and Orion, the LOS kinetic energy
injected by outflows is around 14% of the total LOS kinetic energy. Since
kinetic energy is proportional to $v^{2}$, one caveat here is that a small
amount of mass located at high velocities likely dominates the kinetic energy
estimates. However, due to observational limits, some high velocity gas
emission vanishes into background noise, such that our kinetic energy
calculation is likely an underestimate. As discussed in Xu et al. (2020b),
this could potentially introduce a factor of two uncertainty.

Figure 10: Outflow mass (solid green colors) and the number of YSOs (hatched
blue and red) in the four molecular clouds. Numbers in the hatched bars
indicate the number of YSOs. Numbers in the unhatched bars indicate the gas
mass in units of $10^{3}$ $M_{\mathord{\odot}}$. The total mass indicates the
mass of the host cloud within the area studied.  Figure 11: Outflow line-of-
sight momentum (solid colors) and the number of YSOs (hatched colors ) in the
four molecular clouds. Numbers in the hatched bars indicate the number of
YSOs. Numbers in the unhatched bars indicate the gas momentum in units of
$10^{3}$ $M_{\mathord{\odot}}$ km s-1. The total momentum indicates the
momentum of the host cloud within the area studied. Figure 12: Outflow line-
of-sight kinetic energy (solid colors) and the number of YSOs (hatched colors)
in the four molecular clouds. Numbers in the hatched bars indicate the number
of YSOs. Numbers in the unhatched bars indicate the kinetic energy in units of
$10^{46}$ ergs. The total energy indicates the energy of the host cloud within
the area studied.

###  3.3 Correlation Between the Outflow Properties and the Number of YSOs

In this section, we study the correlation between the total mass, momentum and
energy associated with outflows and the number of YSOs. Due to the extended
and overlapping nature of the outflow emission it is not possible to
definitively associate particular sources with particular outflows. Instead,
we consider an approach that independently counts the number of YSOs and the
outflow impact in a region. We define a “window” to scan through the entire
cloud region. We adopt window sizes of 0.5 pc$\times$0.5 pc, 1 pc$\times$1 pc
and 2 pc$\times$2 pc to scan Ophiuchus, and 2 pc$\times$2 pc, 3 pc$\times$3 pc
and 5 pc$\times$5 pc to scan Taurus and Orion. Since Perseus has the poorest
physical resolution, and its full map is narrow, we only adopt window sizes of
2 pc$\times$2 pc and 3 pc$\times$3 pc to prevent the window exceeding the
region. We examine the effect of our choice of window sizes in Appendix A. We
set a scanning step size such that each box overlaps with its neighbors by at
least 80%. We calculate the mass, momentum and energy inside the window and
count the number of YSOs inside. We consider the outflow relationship between
the number of both young YSOs and total YSOs.

Figure 13 shows the correlation between the mass associated with feedback and
the number of YSOs in the four clouds. The outflow mass grows linearly with
the number of young YSOs for all four clouds. However, the clouds have
different offsets, indicating that the outflow mass per young YSO varies.
Orion has the most mass associated with outflows per young YSO. However, Orion
also has a very large population of older YSOs, which may also drive outflows.
The right panel of Figure 13 shows that the correlation between the outflow
mass and the number of total YSOs is tighter, and most of the separation
between different clouds disappears. There is a linear relation between the
outflow mass and the number of all YSOs, where on average, each YSO is
associated with one $M_{\mathord{\odot}}$ of outflow material. The outflow
mass directly launched by the protostar is expected to be 10-30% of the
accreted gas, i.e., 10-30% of the star mass (e.g., Shu et al., 1988; Pelletier
& Pudritz, 1992). However, here our outflow estimate also counts entrained
mass, which is much higher (e.g., Offner & Chaban, 2017). Our results suggest
one solar mass of outflow material per source. Assuming a one solar mass star
for example, theoretical models for outflow launching predict that the outflow
mass directly launched would be about 0.1-0.3 $M_{\mathord{\odot}}$(e.g., Shu
et al., 1994; Bontemps et al., 1996). This implies that 0.7-0.9
$M_{\mathord{\odot}}$ of the outflow gas we identify is entrained. This gives
a mass-loading factor of 2.3-9. This is comparable to the mass-loading factor
estimated from simulations of individual protostars and the average core-mass-
function to stellar IMF offset (Machida & Hosokawa, 2013; Offner & Arce, 2014;
Offner & Chaban, 2017). It is worth noting that our outflow estimate is larger
than previous observational estimates. The main reason is that casi-3d is able
to capture the outflow material that has relatively low velocities. The low-
velocity outflowing gas entrains a significant amount of mass that also
contributes to the momentum and kinetic energy of the host cloud, but it is
not possible to separate this gas from the ambient gas by eye. Previous work
mainly considered the high-velocity outflow mass only.

Younger sources power stronger outflows (Bally, 2016), so we might expect a
tighter relation between the outflow mass and the number of young YSOs.
However, we see better correlation between outflow mass and the total YSO
number. There are several possible explanations. First, low-mass sources, with
age up to 3 to 4 Myr (Bally, 2016) continue to launch outflows, albeit weaker
ones. We discuss an example in Section 4.1, where an outflow is likely driven
by old evolutionary stage YSOs. Second, although young YSOs have a higher
mass-loss rate, the outflow mass is a combination of both mass-loss rate and
the launching duration. The older YSOs likely have ejected a significant
amount of gas by the time they are observed. In addition, our method may
include outflow “relics.” The dissipation timescale of outflow features
produced by young stars is much longer than the lifetime of the driving
source. These fossil outflows might remain even as their driving sources
evolve to older evolutionary stages. Cunningham et al. (2006) found a similar
scenario in numerical simulations, where fossil outflows retain speeds above
the turbulent velocity for a timescale that is 10 times the duration of their
driving source. Similarly, Offner & Chaban (2017) found in simulations that
although the high accretion and active outflow phase lasts for $\sim$ 0.05
Myr, their impact on the velocity dispersion of the host clouds remains for
several 0.1 Myr.

In our outflow mass estimate, a significant amount of the total outflow mass
derives from emission in the cloud velocity channels ($|v|\lesssim 1$ km s-1).
This raises the concern that our results may be contaminated by dense, non-
outflow material. In Appendix E we evaluate the effect of dense gas
contamination on the derived outflow masses and YSO relations. We show that
our outflow mass estimate is not significantly contaminated by dense gas.

We show the correlation between the outflow momentum and the number of YSOs in
Figure 13. The trend is similar to that of outflow mass. On average, one YSO
injects 1 $\rm{M_{\sun}km/s}$. This indicates the mass weighted LOS velocity
is 1 km/s, which is even smaller than the cloud turbulent velocity derived
from 12CO. This is because the 13CO emitting region dominates the mass of the
cloud and the outflows, which usually has a velocity dispersion only half that
of 12CO. For example, after subtracting the velocity gradients, the 12CO
velocity dispersion in Orion is $\sim$1.5 km/s, while that of 13CO is only 0.8
km/s. The mass-weighted LOS velocity dispersion in Orion is 0.98 km/s. Model
MF predicted outflow gas has a LOS velocity dispersion of 1.02 km/s in Orion.
All these numbers are similar for the four clouds. Although the outflow 12CO
emission spans 2-3 km/s across the velocity channels, most of the mass is
associated with gas within 0.8 km/s of the cloud mean velocity, i.e., the 13CO
emitting regions. Near the rest-frame velocity channels, the gas mass
increases significantly but the fraction of gas associated with feedback
drops. Consequently, there is a competition between these two effects. Xu et
al. (2020b) examined the outflow mass located near the rest-frame velocity,
which was ignored in previous outflow surveys. As shown in Figure 16 and 17 in
Xu et al. (2020b), the rest-frame gas that is associated with the outflows
accounts for almost 75% of the total outflow mass. Thus, the mass-weighted LOS
velocity is relatively low, because high-velocity gas contributes only a small
portion of the total outflow mass. In addition, the resolution and low-signal-
to-noise of the data makes detecting high-velocity outflow emission difficult.
Consequently, the mass weighted LOS velocity of outflows is only 1 km/s for
all four clouds. We show a position-velocity diagram of an outflow in 13CO in
Appendix D.

Our mass-weighted LOS outflow velocity of 1 km/s is low compared to typical
outflow velocities estimated in previous works of 1-4 km/s (e.g. Arce et al.,
2010). This is caused by the mass-weighting. In most cases, observers discard
the rest-frame gas emission, whose outflow morphology is difficult to
identify. This leads to a higher LOS velocity estimate. However, casi-3d is
able to pick out the outflow emission from the confused ambient gas emission
near the rest-frame velocity. Although the fraction of mass associated with
feedback is low, perhaps a few percent, the total mass near rest-frame,
especially in 13CO emission channels, is significantly higher than that
located at high-velocity channels. After correcting for the contamination, the
outflow gas near the rest frame still dominates the total mass, which is
consistent with the conclusion in Xu et al. (2020b). On the other hand, high-
velocity gas is less dense and may be lost in the background noise. Thus, we
likely underestimate the total outflow momentum by at least 10% Xu et al.
(2020b).

The correlation between the outflow kinetic energy and the number of YSOs in
the four clouds is similar to that of mass, as shown in Figure 13.

Figure 13: Correlation between the mass (top row)/momentum (middle row)/energy
(bottom row) associated with feedback and the number of young YSOs (left
panel) or the number of all YSOs (right panel) in the four clouds. The filled
circles indicate the total outflow mass (top row)/momentum (middle row)/energy
(bottom row) and the total number of YSOs (young YSOs in the left panel, all
YSOs in the right panel) in the four clouds. The window sizes are 0.5
pc$\times$0.5 pc, 1 pc$\times$1 pc and 2 pc$\times$2 pc for Ophiuchus, and 2
pc$\times$2 pc, 3 pc$\times$3 pc and 5 pc$\times$5 pc for Taurus and Orion.
The window sizes for Perseus are 2 pc$\times$2 pc and 3 pc$\times$3 pc.

###  3.4 Quantifying the Impact of Feedback with Turbulent Statistics

Feedback, including stellar winds and outflows, injects kinetic energy into
the host cloud. The input energy influences the shape of spatial power
spectrum (SPS) of the integrated intensity map of 12CO. Xu et al. (2020a)
found that stellar wind-generated bubbles flatten the SPS of the 13CO
integrated intensity map, which indicates mass and energy are injected at
small scales. In this section, we investigate how protostellar outflows affect
the SPS of the 12CO integrated intensity maps.

The SPS is defined as the square of the 2D Fourier transform of an image:

$\begin{split}\mathcal{P}(k)&=\sum_{|\vec{k}|=k}|\mathcal{M}_{0}({\vec{k}})|^{2}\\\
&=|\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}M_{0}({\vec{x}})e^{-2\pi
j{\vec{k}\vec{x}}}d{\vec{x}}|^{2},\end{split}$ (5)

where $\mathcal{M}_{0}$ is the 0th moment (integrated intensity) of 12CO.

Figure 14 shows the SPS of Taurus where the emission is above 0.2 K (i.e.,
excluding noise) and the SPS of the model ME1 and MF predicted feedback
regions. Figure 14 shows that the SPS of the feedback predictions breaks into
two power laws. The break point corresponds to a physical scale of $\sim 0.5$
pc, which might indicate the typical outflow mass and energy injection scale.
For reference, we show the scale bar of this injection length in a subregion
of Taurus in Figure 14. This outflow injection scale is comparable to the
typical outflow size in Taurus.

We find similar broken power laws in Ophiuchus, Perseus and Orion. We list the
SPS fits for all four regions in Table 1. For a given cloud, we find the model
ME1 break point scale to be similar and usually slightly higher, than that of
model MF. Table 1 also lists the median of the distance between young YSOs and
their four nearest young companions.

There are several competing effects that likely influence the location of the
break point. For isolated outflows, the break point size is close to the
outflow physical scale. However, the typical outflow size is not necessarily
the same in all regions. It is influenced by the gas density and how readily
the outflow can expand, as well as the source mass and age. In very clustered
regions, outflows interact with each other, since the separation between
sources is smaller. This causes the outflow emission to blend together, which
increases the break point.

The physical scale of the break point is small for Orion and Ophiuchus, and
Taurus has a intermediate break point scale, while Perseus has the largest.
Among the four clouds, Taurus has the most distributed and least clustered
YSOs. Taurus also has the lowest average gas density. The mass/energy
injection scale in Taurus reflects the physical size of individual molecular
outflows. While in Orion and Ophiuchus, the average distance between stars is
small, which increases the probability that outflows from different sources
interact. The mean gas densities are also higher, and high external pressure
from the surrounding gas may also act to limit the propagation of outflows.
For example, Kirk et al. (2017) found all dense cores in Orion are pressure
confined. All of these effects help to explain the relatively small
mass/energy injection scale.

The physical scale of the break point in Perseus is significantly larger than
that of the other regions. This might be due to both its typical density,
which is lower than that of Orion and Ophiuchus, and the presence of multiple
clusters in Perseus. Perseus hosts several intermediate mass star-forming
clusters, such as NGC 1333 and IC348. Outflows in these clusters interact and
blend together, such that mass/energy injection occurs on a relatively larger
scale compared to individual stars. The average young YSO separation is
slightly larger in Perseus than that in Orion and Ophiuchus, in part because
the two main clusters contain few young YSOs. Consequently, the combination of
these effects may explain the large mass/energy injection scale in Perseus.

To investigate how YSO clusters affect the SPS break point, we conduct a SPS
analysis towards two star clusters, NGC1333 in Perseus and L1688 in Ophiuchus.
We adopt different window sizes to explore the location of the break point in
the cluster when viewed on different size scales. When the window size is
smaller than the cluster size, we find there is no break point in the SPS for
both clustered region and the slope is steeper. When the window size is two or
more times larger than the cluster size, the break point appears. However, the
break point scale is not neatly correlated with the size of the cluster nor
the window size. NGC1333 has a break point between 0.6 pc and 0.7 pc, while
L1688 has a break point scale between 0.3 pc and 0.4 pc. These scales are
similar to those of the full predictions, which suggests that the outflows in
the clustered regions are influencing the location of the overall break point.

We also explore the effect of different observation resolutions on the
results. The observations of Ophiuchus, Taurus and Orion have a similar
physical scale per pixel, $0.013-0.015$ pc/pixel. However, the physical scale
per pixel of Perseus is twice that of the others. We convolve the Ophiuchus
and Orion observations to a similar effective resolution as Perseus, and apply
the casi-3d models to these convolved data cubes and conduct the same SPS
analysis on the convolved data. However, we find only minor changes in the
break point. For example, the break point remains the same in the convolved
Ophiuchus model MF prediction and increases by 0.1 pc in the convolved
Ophiuchus model ME1 prediction. However, it decreases by 0.08 pc in the
convolved Orion ME1 prediction, while increasing by 0.08 pc in the Orion MF
prediction. We find no significant trend when convolving with a larger beam
and conclude the analysis is relatively insensitive to the resolution.
However, this test suggests that the break point has an uncertainty of $\pm
0.1$ pc, which is smaller than the difference between the break point of
Perseus and those of the other clouds.

Finally, to test whether (non-feedback) gas structures influence the break
point location, we compute the SPS for the dense regions of the clouds that
are traced by 13CO. We mask out the low intensity pixels to isolate only the
dense gas. We compare the results for two clouds, Orion and Taurus. When all
pixels with intensities below 5 K are removed, there is no break point in the
Orion SPS. For the densest regions, i.e., as defined using an intensity cutoff
of 10 K, we find a break point at 0.6 pc. This is two times larger than that
of the break point of the outflow maps. We find a similar result in Taurus.
When applying a high intensity cut-off, we find a break point around 1 pc,
which is a size-scale twice that of the break point identified by the Taurus
SPS outflow analysis. We conclude that our outflow predictions are not
correlated with dense cloud structures.

Figure 14: Upper panel: intensity of 12CO (1-0) integrated over all velocity channels for a subregion in Taurus, overlaid with the model MF prediction (red and blue contours). The physical scale of the break point in the SPS is also shown for reference. Lower panels: the spatial power spectrum (SPS) of the regions in Taurus where the emission is above 0.2 K (excluding noise) and the SPS of the model ME1 and MF predicted feedback regions. The top right number in each panel indicates the slope of the fit. The middle bottom number in each panel indicates the physical scale of the broken point.  Table 1: Fitting Results of the Spatial Power Spectrum Cloud | Full Mapa | ME1b | MFb | Young YSO
---|---|---|---|---
Slope | Slope | Break (pc) | Slope | Break (pc) | Seperationc (pc)
Ophiuchus | -3.17$\pm$0.02 | -2.47$\pm$0.05 | 0.36 | -1.82$\pm$0.05 | 0.38 | 0.24
-3.71$\pm$0.07 | -3.78$\pm$0.08
Taurus | -2.92$\pm$0.02 | -1.59$\pm$0.02 | 0.56 | -1.49$\pm$0.02 | 0.49 | 0.85
-3.26$\pm$0.04 | -3.47$\pm$0.04
Perseus | -3.07$\pm$0.02 | -2.18$\pm$0.04 | 1.03 | -2.21$\pm$0.03 | 0.65 | 0.31
-3.40$\pm$0.06 | -3.47$\pm$0.05
Orion | -2.95$\pm$0.02 | -2.29$\pm$0.03 | 0.36 | -2.06$\pm$0.03 | 0.27 | 0.22
-3.75$\pm$0.05 | -4.83$\pm$0.05
Notes:
a Single power-law fit results for the spatial power spectra for the emission
regions (excluding noise).
b Broken power-law fit results for the spatial power spectra applied to the
ME1 and MF feedback maps.
c Median of the separation between YSOs and their four nearest companions.

##  4 Discussion

###  4.1 Outflows without Driving Sources and False Detections

Confirming outflow identifications requires ancillary data, e.g., YSO
catalogs. With the help of YSO locations, we are able to increase confidence
in the casi-3d predictions, which are based on 12CO data only. In this
section, we discuss two cases predicted by casi-3d: a region in Taurus with no
detected YSOs, and a region in Orion with no young YSOs.

Figure 15 presents the channel by channel prediction by model MF on a region
with no YSOs in Taurus. This region is considered to be a trans-Alfvenic flow
regulated by magnetic fields (Heyer & Brunt, 2012; Heyer et al., 2016), which
has low surface brightness 12CO emission. There are no known YSOs observed by
GAIA or Herschel (Roccatagliata et al., 2020). Consequently, we believe there
is no outflow activity. However, model MF predicts outflow activity in this
region. Inspection of the channel map shows that a coherent jet-like structure
exists across several velocity channels. These structures actually visually
resemble outflows, which might be reason for the model failure.

Figure 16 shows the model MF prediction along a position-velocity cut of the
fake outflow in Figure 15. In the position-velocity diagram, we identify two
faint coherent high velocity blobs between 4 and 5 km/s, Their morphology is
similar to some of the faint outflows (e.g., TMO22 in Li et al., 2015). Since
these 12CO structures are indistinguishable from true outflows, casi-3d is not
able to recognize these false detections based only on 12CO morphology. In
terms of morphology, casi-3d performs robustly in identifying coherent high
velocity features that are similar to outflows. But many mechanisms may cause
coherent high velocity features, including but not limited to cloud formation,
cloud-cloud collision, gas phase transition near the cloud boundary or gas
flows regulated by MHD waves (Heyer & Brunt, 2012; Motte et al., 2014;
Nakamura et al., 2014; Heyer et al., 2016).

Figure 17 shows an example of the performance of models ME1 and MF toward a
region where there are no young YSOs in Orion. In the position-velocity
diagram, we identify a clear outflow-like feature above 10 km/s, which is also
identified by the two models. We notice that there are several old
evolutionary stage YSOs nearby, which might act as driving sources. One reason
might be that the boundary between young and old is not so well defined, which
is discussed in Section 3.3. Or these outflow structures might exist for a
much longer timescale than that of their driving sources. We recognize that
there are two velocity components in the molecular cloud at this position. A
narrow gas bridge located around 8 km/s connects these two gas components,
which is likely an outflow. We cannot confirm the origin of this high velocity
component. It might be caused by the fossil outflows from YSOs or by two
converging gas flows. All these mechanisms cause similar high-velocity
features, which cannot be easily distinguished either visually or using
casi-3d.

In order to validate that our casi-3d models are sensitive to gas velocity
structures and morphology but not only the dense regions, we examine the
casi-3d performance on several previously identified filaments in Taurus.
Panopoulou et al. (2014) identified 10 filametary structures in 13CO emission
in Taurus, and illustrated the PPV diagrams of the filaments in Figure 15 in
their paper. The filaments have very coherent motions and exhibit small
velocity dispersions. Their filaments 3, 4, 5, 6, and 7 do not have
significant high-velocity components. When comparing with our outflow
predictions (Figure 4 and 5 in this work) and bubble predictions (Figure 17
and 18 in Xu et al., 2020a), we find that both models predict little feedback
along these filaments. However, when we look at Filament 2 (also known as
L1495/B213), we see a clear high-velocity structure with a “U” shape. This
morphology is consistent with the theoretical bubble morphology shown in
Figure 5 in Arce et al. (2011). Unsurprisingly, both casi-3d models that
identify bubbles and outflows predict the presence of feedback at the location
of filament 2. casi-3d may return a false detection if the morphology of a
structure is similar to that of outflows or bubbles. In this case, visual
inspection suggests that filament 2 also contains a bubble, so it is likely
that the feedback identified by our models at this location is real. We
emphasize that machine learning models are not perfect tools. They must be
applied with caution and checked thoroughly.

Figure 15: 12CO channel map of a region with no YSOs in Taurus overlaid with
the prediction by model MF in white contours.  Figure 16: Position-velocity
diagram of 12CO emission toward a region with no YSOs in Taurus. Left panel:
integrated intensity of 12CO over the the full velocity range (from -1.5 km/s
to 13.4 km/s) overlaid with the model ME1 and MF predictions in white
contours. The purple line illustrates the cut direction of the position-
velocity diagram. Middle and right panel: position-velocity diagram of 12CO
emission overlaid with the model ME1 and MF predictions in white contours.
Figure 17: Position-velocity diagram of 12CO emission toward a region without
young YSOs in Orion. Left panel: integrated intensity of 12CO over the full
velocity range (from 2.3 km/s to 17.1 km/s) overlaid with the model ME1 and MF
predictions in white contours. Letters “O” mark old YSO positions. The purple
line illustrates the cut direction of the position-velocity diagram. Middle
and right panel: position-velocity diagram of 12CO emission overlaid with the
model ME1 and MF predictions in white contours.

###  4.2 Outflows versus Bubbles

In this section we investigate the overlap between a model identifying outflow
feedback and one identifying stellar wind feedback. We present the prediction
by two sets of models: those trained to identify protostellar outflows and
those trained to identify stellar wind driven bubbles.

The morphology of bubbles is more symmetric and arc-like compared to that of
outflows. However, when the line of sight is parallel to an outflow launching
axis, i.e., looking through the cavity of the outflow, the morphology of the
outflow can resemble that of a bubble. Figure 18 illustrates an example of a
synthetic outflow whose launching axis is parallel to the line of sight.
Consequently, models trained to identify outflows are likely to also identify
some bubble structures.

Figure 19 shows the prediction by models trained to identify stellar wind
driven bubbles towards Ophiuchus L1688 Region. Comparison with Figure 3 shows
that the prediction by models trained to identify protostellar outflows covers
most of the regions that are identified as bubbles. Quantitively, 68% of the
voxels are identified as feedback by both MF models, 21% of the voxels are
identified as only outflow feedback, while 11% of the voxels are identified as
only belonging to bubbles. Statistically, we might neglect 13% of the feedback
gas associated with stellar winds when adopting the mass estimates predicted
by models trained to identify protostellar outflows. On average, the outflow
mass predicted by models can serve as a decent estimate for feedback mass.

To better investigate the difference between outflows and bubbles, we employ
an unsupervised learning method t-SNE (t-distributed Stochastic Neighbor
Embedding) to study the similarities between the spectra of outflows and that
of bubbles. We discuss t-SNE and the results in detail in Appendix B. Our main
finding is that there is a variation in spectra shape for outflows and
bubbles. And even the spectra that are dominated by outflow gas show various
spectral shapes. Meanwhile, there are some “transitional spectra” with a shape
similar to both. We note that the t-SNE representation encodes only 1D
information about the spectra and does not represent any correlations in the
spatial dimension. Consequently, this analysis provides intuition for
understanding the features that the models select when they identify feedback,
and it illustrates the challenge of distinguishing outflows and bubbles based
only on the LOS spectrum.

Figure 18: A synthetic outflow whose launching axis is parallel to the line of
sight.  Figure 19: Position-velocity diagram of 12CO emission towards
Ophiuchus L1688 Region, similar to Figure 3 but predicted by models trained to
identify stellar wind driven bubbles.

###  4.3 Impact of Outflows on Molecular Clouds

In this section, we place our energy estimations (Section 3.2) in the context
of prior work and discuss the broader implications of the global impact of
outflows. Most prior work derives the total outflow energy by assuming a
moderate inclination angle and multiplying the 1D energy estimate by a factor
of 3. For example, Dunham et al. (2014a) adopted an inclination angle of
57${}^{\circ}.$3 which yields a energy correction factor of 3.4. Here we make
the same assumption and report the 3D total, where we multiply our 1D estimate
by a factor of 3.

To investigate the impact of outflows on turbulence, a typical approach is to
compare the turbulence dissipation rate and the outflow energy injection rate
(e.g., Arce et al., 2010; Li et al., 2015; Feddersen et al., 2020). However,
the outflow energy injection timescale is highly uncertain. As discussed in
Section 3.3, we find the outflow mass is linearly correlated with the total
number of YSOs instead of only the younger sources, which indicates that the
driving timescale of outflows is even less clear. Instead, we compute the
total outflow kinetic energy – a quantity estimated as part of previous
analyses in the determination of the outflow injection rate – and compare our
new values to those in the prior published studies. This allows us to conclude
whether the impact of outflows is relatively similar to, more or less than
that in prior analyses.

Nakamura et al. (2011) conducted observations of 12CO (3-2) and 12CO (1-0)
towards an active cluster-forming clump in Ophiuchus and identified six
molecular outflows in both data. Nakamura et al. (2011) derived the outflow
mass from 12CO (3-2) and 12CO (1-0) emission located in the high-velocity line
wings. They derived a total outflow kinetic energy of $6\times 10^{44}$ ergs.
After comparing the total outflow energy injection rate with the dissipation
rate of the supersonic turbulence, Nakamura et al. (2011) concluded that
outflows inject significantly more energy than needed to offset the
dissipation of turbulence. In our work, we find that the total outflow kinetic
energy is $9\times 10^{45}$ ergs, which is an order of magnitude larger than
the value in Nakamura et al. (2011). A larger area coverage in our work and
using the optically thin tracer 13CO in mass estimates might explain this
difference. 12CO (3-2) and 12CO (1-0) are likely optically thick, which could
cause their outflow mass, momentum and energy to be underestimated. More
importantly, our method does not discard the outflow material ejected
perpendicular to the line of sight that is located near the cloud central
velocity. Since our energy calculation leads to a higher energy estimate, our
result confirms the conclusion that outflow kinetic energy is sufficient to
compensate for the dissipation of turbulence in Ophiuchus.

Li et al. (2015) conducted a feedback census in Taurus and identified 55
outflows and 37 bubbles. They derived a total outflow energy of $3.9\times
10^{45}$ ergs, which is 1% of the Taurus cloud turbulent energy. However, Li
et al. (2015) found that bubbles dominate the feedback energy. When also
combining the kinetic energy from bubbles, they found the total kinetic energy
from stellar feedback in Taurus is 30% of the cloud turbulent energy. Since
the energy injection rate from stellar feedback is comparable to the
dissipation rate of the cloud turbulence, they concluded that stellar
feedback, mainly in the form of bubbles, is sufficient to maintain turbulence
in the current epoch. However, Xu et al. (2020a) showed that the Li et al.
(2015) bubble energy estimate was substantially overestimated due to LOS
contamination, which is caused by gas emission that is not associated with
feedback being included in the outflow/bubble estimate. After correction, the
kinetic energy from bubbles decreased by a factor of four. In our current
work, we find the total outflow kinetic energy is $5.4\times 10^{46}$ ergs,
which is an order of magnitude larger than the value in Li et al. (2015). This
difference is likely caused by a similar reason to that in Ophiuchus, i.e.,
our method includes the outflow material located around the cloud central
velocity that is excluded in traditional approaches. Our calculation indicates
that the kinetic energy from outflows is 14% of the cloud turbulent energy. If
we combine the contribution from bubbles (Xu et al., 2020a), the total kinetic
energy from stellar feedback is 22% of the cloud turbulent energy. Li et al.
(2015) concluded that the outflow energy injection rate is marginally
comparable to the dissipation rate of turbulence in Taurus, although outflow
energy is only 1% of the Taurus cloud turbulent energy. With our updated
feedback energy, which is 22% of the cloud turbulent energy, we conclude that
the feedback energy injection rate can compensate for the turbulent
dissipation rate (Li et al., 2015; Xu et al., 2020a).

Arce et al. (2010) identified 60 outflow candidates in Perseus and derived a
total outflow energy of $2\times 10^{46}$ ergs, which is 12.5% of the total
cloud turbulent energy. However, it is worth noting that the method of
calculating cloud turbulent energy used in Arce et al. (2010) is different
from our method. Arce et al. (2010) adopted an average line width of 2 km/s to
calculate the turbulent energy. While in this work, we calculate the turbulent
energy channel by channel followed by Equation 3. The cloud turbulent energy
derived by Arce et al. (2010) is $1.6\times 10^{4}7$ ergs, while our approach
indicates an estimate of $8\times 10^{4}7$ ergs, which is a factor of 5 larger
than the estimate by Arce et al. (2010). To make a fair comparison, we adopt
our turbulent energy estimate to discuss the impact of outflows in Perseus.
Under these circumstances, the kinetic energy of the previously identified
outflows are only 2.4% of the total cloud turbulent energy. Xu et al. (2020b)
applied casi-3d to the same outflow catalog and predicted a total outflow
energy of $7.8\times 10^{45}$ ergs, which is one fourth of the value in Arce
et al. (2010). This is caused by the lower mass-weighted LOS velocity as
discussed in Section 3.3. In our present work, we derive a total outflow
energy of $1.2\times 10^{47}$ ergs, which is 6 times larger than that in Arce
et al. (2010). There are several possible reasons for this higher value.
First, there are newly identified isolated outflows as discussed in Xu et al.
(2020b). It should be noted that Xu et al. (2020b) just followed up the 60
previous outflows and identified a few new outflows in the vicinity of these,
whereas our current work make a prediction using the full CO map, which
includes the entire cloud and a couple of clusters. This likely increase the
outflow kinetic energy by a factor of a few. Second, casi-3d is able to
identify all feedback structures in clustered regions, which are not fully
identified in previous work. This also may enlarge the outflow energy by a
factor of a few. Moreover, outflows identified by Arce et al. (2010) are all
distinct outflows with obvious coherent high-velocity structures that are
likely powered by strong driving sources. However, in our work, casi-3d also
identifies less distinct outflows whose velocity is small compared to
previously identified ones. These weaker, less distinct outflows might be
driven by more evolved or lower mass sources. This indicates that casi-3d
provides a more inclusive outflow sample. Arce et al. (2010) concluded that
the 60 outflow candidates they identified play an important role in
maintaining turbulence at the current epoch in Perseus, although the kinetic
energy of outflows is 12.5% of the total cloud turbulent energy in their
study, or 2.4% of the total cloud turbulent energy in our turbulence
calculation. However, our updated outflow kinetic energy estimate is 15% of
the cloud turbulent energy. This indicates that outflows still play an
important role in maintaining turbulence in Perseus.

Feddersen et al. (2020) identified 45 outflows near Herschel Orion Protostar
Survey (HOPS) protostars (Furlan et al., 2016) in Orion but skipped the OMC-1
region. They derived a total outflow energy of $0.7-1.7\times 10^{46}$ ergs.
They found that the total outflow energy injection rate is comparable to the
dissipation rate of turbulence. In our work, we derive a total outflow energy
of $2.1\times 10^{47}$ ergs, which is an order of magnitude larger than
previous estimates. The main reason for the difference is that our survey
covers a larger area and does not only focus on gas near the position of YSOs.
In particular, we include the most active star forming region OMC-1. casi-3d
predicts the kinetic energy from outflows in OMC-1 is $7.7\times 10^{46}$
ergs, which is 37% of the entire outflow energy in Orion. Thus, we confirm the
prior conclusion that outflows are sufficient to maintain turbulence in Orion.

In conclusion, our outflow kinetic energy estimates for all four clouds are
larger than previous estimates. Consequently, we confirm the previous
conclusion that outflows are an important agent to maintain turbulence at the
current epoch.

###  4.4 Uncertainties of Outflow Estimates

In Section 3.2, we calculate the masses and dynamical properties of the
outflows for the four clouds. In this section, we attempt to quantify the
uncertainty of these outflow estimates.

First, we assume a moderate excitation temperature of 25 K or the 12CO peak
temperature, whichever is higher, to calculate the outflow mass. Most studies
adopt values of excitation temperature in the range of 10-50 K (Dunham et al.,
2014a; Feddersen et al., 2020). Since the mass scales linearly with the
excitation temperature, the choice of excitation temperature introduces a
factor of two uncertainty.

Meanwhile, Xu et al. (2020b) found that model MF predicts the outflow mass
within a scatter of 0.41. The most extreme offset case over-/under- estimates
the mass by a factor of two. We did not find any offset that would cause a
systematic over or under-estimation. Meanwhile, we carry out one additional
test to constrain the degree of possible contamination from cloud emission at
rest-frame gas velocities. Appendix F shows an analysis of a more complex
synthetic observation, which contains a number of interacting outflows. We
conclude our casi-3d model is generally able to exclude contamination from the
ambient gas near the rest frame velocity and may even underestimate the
outflow contribution from low velocities in complex star-forming areas.

Furthermore, high-velocity, low-density gas emission is easily missed in high-
velocity channels due to the low signal-to-noise. Xu et al. (2020b) used the
simulations to estimate that 10% of the outflow gas is lost. The correction
factors for LOS momentum and LOS energy are 1.3 and 1.8 (Xu et al., 2020b),
respectively.

In addition, the inclination angle of the outflow plays an important role in
converting LOS momentum/energy to 3D momentum/energy. In this work, we adopt
correction factors of $\sqrt{3}$ and 3 to convert LOS momentum and LOS energy
to 3D estimates. These correspond to an inclination angle of 55∘. If the
inclination angle is between 20∘ and 70∘, this leads to a factor of 2
uncertainty for 3D momentum and a factor of 4 uncertainty for 3D energy.

Finally, chemical conditions, i.e., 12CO and 13CO abundances, have an
influence on the mass estimate. In our mass calculation, we adopt constant
abundance ratios, [12CO]/[H2]=$10^{-4}$ and [12CO]/[13CO]=62. However, these
abundance ratios vary across diffident clouds. We cannot evaluate this
uncertainty without sophisticated chemical modeling for each cloud.

Consequently, by combining these uncertainties and using error propagation, we
find a factor of 3 uncertainty for mass, a factor of 3.4 uncertainty for 3D
momentum, and a factor of 5 uncertainty for 3D energy.

##  5 Conclusions

We apply the deep learning method casi-3d to four nearby molecular clouds,
Ophiuchus, Taurus, Perseus and Orion, to systemically identify protostellar
outflows and study the impact of outflows on their host clouds. Our main
findings are the following:

  1. 1.

The total outflow masses are 267 $M_{\mathord{\odot}}$, 795
$M_{\mathord{\odot}}$, 1305 $M_{\mathord{\odot}}$ and 6332
$M_{\mathord{\odot}}$ for Ophiuchus, Taurus, Perseus and Orion, respectively.
On average, the mass associated with feedback is around 10% of the host cloud
mass for all four clouds.

  2. 2.

The total 3D outflow energies are $9\times 10^{45}$ ergs, $6\times 10^{46}$
ergs, $1.2\times 10^{47}$ ergs and $6\times 10^{47}$ ergs for Ophiuchus,
Taurus, Perseus and Orion, respectively.

  3. 3.

The outflow mass is linearly proportional to the total number of YSOs for all
four clouds. On average, each YSO is associated with one solar mass of outflow
material.

  4. 4.

The outflow momentum is linearly proportional to the total number of YSOs for
all four clouds. On average, each outflow has a mass-weighted LOS velocity of
1 km/s. The relatively low value is because more mass is located near the
cloud central velocity.

  5. 5.

We compute the spatial power spectrum of the outflow prediction map. We find
all four clouds exhibit a break point, which ranges from 0.27 pc (Orion) to
0.65 (Perseus) pc. The break point likely indicates the typical outflow mass
and energy injection scale.

  6. 6.

Models trained to identify outflows are likely to also identify bubble
structures, which may be produced by main sequence massive stars.

  7. 7.

We compare the energy associated with outflows to the rate of turbulent
dissipation and conclude that feedback is sufficient to maintain turbulent
dissipation at the current epoch for all four clouds.

D.X., S.S.R.O., and R.A.G. acknowledge support by NSF grant AST-1812747. D.X.
acknowledges support from David Alan Benfield Memorial Scholarship in
Astronomy and from the Virginia Initiative on Cosmic Origins (VICO). S.S.R.O.
acknowledges support from NSF Career grant AST-1748571. R.A.G. acknowledges
support from NASA ADAP grant NNX17AF24G. HGA acknowledges support from NSF
grant AST-1714710. The Texas Advanced Computing Center (TACC) at the
University of Texas at Austin provided HPC resources that have contributed to
the research results reported within this paper.

##  Appendix A Discussion on the Effect of Box Sizes on the YSO-Mass Relation

In this section, we examine the effect of our choice of “window” sizes when
studying the correlation between the total mass associated with outflows and
the number of YSOs. We define windows with physical sizes of 0.5 pc$\times$0.5
pc, 1 pc$\times$1 pc, 2 pc$\times$2 pc, 3 pc$\times$3 pc and 5 pc$\times$5 pc
to scan through the entire cloud region. We set a step size such that each box
overlaps with its neighbors by at least 80%. Then we calculate the mass inside
the window and count the number of YSOs enclosed. We consider the outflow
relationship between the number of both young YSOs and total YSOs.

Figure 20 shows the correlation between the model MF outflow prediction and
the number of young YSOs for the five window sizes applied to the four clouds.
For Ophiuchus, Perseus and Orion, when the window size is above certain
threshold, the correlations between outflow mass and the number of YSOs become
similar, which indicates the choice of window size does not significantly
change the results. Taurus is more sensitive to the window size than the other
regions. This is likely caused by the large separations between YSOs. However,
it still shows a clear trend of smaller separation between different
correlations when the window size gets larger.

The trends are similar when considering total YSOs. Consequently, we conclude
that the study of the correlation between the outflow properties and the
number of YSOs is robust.

Figure 20: Correlation between the mass associated with feedback and the
number of young YSOs. The different colors indicate different window sizes for
the estimation. The filled circles indicate the total outflow mass and the
total number of young YSOs in the four clouds.

##  Appendix B Spectra of Outflows and Bubbles

In this section, we employ an unsupervised learning algorithm t-SNE
(t-distributed Stochastic Neighbor Embedding) to compare the spectra
identified as outflow feedback to those identified as belonging to bubbles.
t-SNE is a tool to visualize high-dimensional data by converting the
similarities between data points into a low-dimensional manifold. It has been
used for classification and outlier detection in a variety of astronomy data
(Jofré et al., 2017; Reis et al., 2018; Lochner & Bassett, 2020; Fluke &
Jacobs, 2020). For example, Reis et al. (2018) successfully applied t-SNE to
cluster stars based on their spectra and identify outliers based on the t-SNE
map.

We take the L1688 star cluster region in Ophiuchus as an example. Figures 3
and 19 show the prediction of this star cluster region by models trained to
identify only outflows and only bubbles, respectively. In some cases the
models clearly identify only one or the other type of feedback, however, there
is also significant overlap. In order to classify the spectra, we calculate
outflow mass and bubble mass of each spectra. If the outflow mass or bubble
mass for a given line of sight is above 5$\sigma$ compared with the background
noise, we consider that spectrum to contain feedback. If the outflow mass is
larger than the bubble mass, we label it an outflow spectrum. If the bubble
mass is larger than the outflow mass, we label it a bubble spectrum. If there
is emission but neither outflow mass nor bubble mass exceeds the 5$\sigma$
threshold, or the spectrum has less than three channels that contain feedback,
we label this spectrum as non-feedback. We neglect any spectra that have no
emission above 3$\sigma$ sigma.

We use the spectra of the L1688 region at each line of sight as input to
t-SNE. t-SNE projects the higher dimensional data into a lower dimensional
(2D) space and organizes it based on similarities (e.g., the line shape and
intensity of a spectrum) in the high-dimension space. Figure 21 shows the
t-SNE clusters. We color the data points based on the casi-3d predictions,
where spectra predicted to be primarily associated with outflows are in red,
those primarily associated with bubbles are in blue and non-feedback in black.
It is clear that the spectral shape for bubbles and outflows are different.
The spectra of bubbles show two Gaussian components with two peaks of
different intensities. While the spectra of outflows have various shapes,
e.g., broader line wings, two Gaussian components but with almost equal peaks.
It’s worth noting that broad line wings are often used to visually identify
outflows(Li et al., 2015; Tanabe et al., 2019). t-SNE separates some bubble
spectra from the outflow spectra. Meanwhile, we notice that there are some
blended/transiting regions, where outflow spectra and bubble spectra are
clustered together. In these cases, we find that bubble and outflow
identifications have similar spectral shapes. This demonstrates that we cannot
distinguish outflows and bubbles based only on single spectra. There is also a
mixed region where non-feedback spectra are blended with outflow spectra. The
spectral shape of this region indeed show less outflow features, where there
is a faint broad line-wing on the blue shifted side.

This analysis explores the similarity between outflow and bubble feedback
using an unsupervised machine learning approach. In some cases, the casi-3d
predictions are clearly distinct. In other cases, where both models identify
feedback emission, it is clear that the spectra are truly similar. This may be
because both types of feedback can produce such velocity features or because
such regions genuinely contain both types of feedback. In the current analysis
it is not possible to distinguish between these two possibilities. However,
this t-SNE clustering gives confidence that feedback is shaping these spectra
and that casi-3d is correctly identifying feedback signatures.

Figure 21: t-SNE map of the 12CO spectra in the L1688 star cluster region in
Ophiuchus. Red dots indicate outflow spectra. Blue dots indicate bubble
spectra. Black dots indicate non-feedback spectra. The sub-plots show the
average spectra of each subregion.

##  Appendix C YSOs in Orion

Orion has an extensive population of older YSOs, so in Figure 8 we plot only
the young YSOs for readability. Here, Figure 22 shows the location of both
young and old YSOs in Orion. While we expect younger YSOs to contribute more
feedback per source, the large number of older YSOs indicates that they likely
contribute feedback throughout the cloud, and this is consistent with the
prediction in Section 3.1.4.

Figure 22: Intensity of 12CO integrated over all velocity channels for Orion.
Letters “Y” and “O” mark YSO positions, as described in Section 2.3.

##  Appendix D Outflows in 13CO

In some regions 12CO is optically thick and we combine both 12CO and 13CO to
derive the kinematic estimates rather than using 12CO only (e.g., Arce et al.,
2010). In most cases, the velocity channels where there is 13CO emission
dominate the mass in each spectrum. Here we compare the identification from
the 12CO map to the 13CO emission. Figure 23 shows the 13CO emission of the
outflow in Figure 9. We can clearly see the velocity range of 13CO is
substantially smaller than that of 12CO. However, the mass in the 13CO
emitting region dominates the total mass of the cloud. This also explains why
the mass weighted LOS velocity is smaller than the velocity calculated based
on 12CO only.

Figure 23: Same as Figure 9 but where the background colorscale indicates the
13CO integrated intensity. White contours show the prediction by models ME1
and MF on the 12CO data.

##  Appendix E Correlation Between the Gas Mass and the Number of YSOs

In this appendix, we examine the relationship between the number of YSOs, the
total gas mass, and the dense gas mass. Prior work has shown that the number
of YSOs is well-correlated with the gas column density, which is
$\Sigma_{*}\propto\Sigma_{\rm gas}^{2}$, or $M_{\rm gas}\propto N_{*}^{0.5}$
(e.g., Gutermuth et al., 2011; Pokhrel et al., 2020, 2021), where the gas mass
was estimated from dust emission. Pokhrel et al. (2020, 2021) found that the
surface density of protostars (i.e., younger YSOs) are tightly correlated with
gas column following the relation above. Gutermuth et al. (2011) found a
similar but less tight relation when considering all YSOs. Consequently, we
expect YSOs and dense gas to be correlated. Figure 24 shows a nearly linear
relation between the total mass and the number of all YSOs. We find the slope
is 0.93, which is slightly smaller than the slope of the outflow mass-YSO
relation, 1.03.

To check whether the outflow mass estimate is correlated with the dense gas
due to contamination, we examine the relation excluding the rest-frame
velocity gas, which is located within the turbulent velocity range ($\pm 1$
km/s). Figure 25 shows the correlation between the total high velocity gas
($|v-v_{cen}|>1$ km/s) mass and the number of all YSOs, and the correlation
between the high velocity outflow mass ($|v-v_{cen}|>1$ km/s) and the number
of all YSOs. The slopes are similar to that of the total gas/outflow gas
relation including the rest-frame velocity gas. Therefore, we conclude the
correlation between gas mass and YSOs exists because outflows are linearly
proportional to the number of YSOs launching them. YSOs form from dense gas
and thus by nature are more likely to be located in denser regions. By
focusing only on the high-velocity material, which traces the outflows most
directly, we show that the trend arises directly from the expected correlation
between YSOs and outflows rather than indirectly from contamination by the
dense gas. This result provides further evidence that our method performs well
in dense regions where contamination from the cloud would otherwise be a
serious problem.

Figure 24: Left: correlation between the total mass and the number of all
YSOs. Right: correlation between the outflow mass and the number of all YSOs.

Figure 25: Left: correlation between the total high velocity gas
($|v-v_{cen}|>1$ km/s) mass and the number of all YSOs. Right: correlation
between the high velocity outflow mass ($|v-v_{cen}|>1$ km/s) and the number
of all YSOs.

##  Appendix F casi-3d Prediction on a Complex Synthetic Observation

In this appendix, we examine a prediction by the two casi-3d models on a
synthetic observation of a simulation snapshot with a number of interacting
outflows. This synthetic observation and the simulation it derives from is not
included in our previous training or testing, so it constitutes a more
challenging test of the model performance. This simulation is run using a
different initial turbulent seed and with twice the inital gas density as the
other simulations in the training set. The mean magnetic field strength is 0.8
$\mu$G. The evolutionary time of the snapshot is 0.7 $t_{\rm ff}$. There are
13 stars launching outflows in the simulation box, whereas our training set is
constructed from snapshots with 1 to 5 stars. Consequently, this output
provides an independent check on the casi-3dmodel performance.

Figure 26 shows the integrated CO emission, tracer fractions and predictions.
We find the casi-3d model does not actually over-predict the amount of
emission in the cloud velocity channels but instead slightly under-predicts
the location and mass of the outflows. More quantitatively, model MF
identifies 76% of the outflow emission at high-velocity channels (above 2km
s-1). However, model MF only identifies 40% outflow emissions near the rest-
frame velocity (within 2 km s-1). This implies casi-3d might underestimate the
outflow mass by a factor of 2. This is consistent with the model uncertainty
discussed in Section 4.4. This analysis gives further confidence that casi-3d
correctly excludes contamination from the ambient gas near the rest frame
velocity, even in relatively clustered and complex regions.

Figure 26: Synthetic 12CO observations and casi-3d prediction in different
velocity channels.

## References

  * Arce et al. (2011) Arce, H. G., Borkin, M. A., Goodman, A. A., Pineda, J. E., & Beaumont, C. N. 2011, ApJ, 742, 105, doi: 10.1088/0004-637X/742/2/105
  * Arce et al. (2010) Arce, H. G., Borkin, M. A., Goodman, A. A., Pineda, J. E., & Halle, M. W. 2010, ApJ, 715, 1170, doi: 10.1088/0004-637X/715/2/1170
  * Bally (2016) Bally, J. 2016, ARA&A, 54, 491, doi: 10.1146/annurev-astro-081915-023341
  * Bontemps et al. (1996) Bontemps, S., Andre, P., Terebey, S., & Cabrit, S. 1996, A&A, 311, 858
  * Carroll et al. (2009) Carroll, J. J., Frank, A., Blackman, E. G., Cunningham, A. J., & Quillen, A. C. 2009, ApJ, 695, 1376, doi: 10.1088/0004-637X/695/2/1376
  * Cunningham et al. (2006) Cunningham, A. J., Frank, A., Quillen, A. C., & Blackman, E. G. 2006, ApJ, 653, 416, doi: 10.1086/508762
  * Cunningham et al. (2018) Cunningham, A. J., Krumholz, M. R., McKee, C. F., & Klein, R. I. 2018, MNRAS, 476, 771, doi: 10.1093/mnras/sty154
  * Dunham et al. (2014a) Dunham, M. M., Arce, H. G., Mardones, D., et al. 2014a, ApJ, 783, 29, doi: 10.1088/0004-637X/783/1/29
  * Dunham et al. (2014b) Dunham, M. M., Stutz, A. M., Allen, L. E., et al. 2014b, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 195, doi: 10.2458/azu_uapress_9780816531240-ch009
  * Feddersen et al. (2020) Feddersen, J. R., Arce, H. G., Kong, S., et al. 2020, ApJ, 896, 11, doi: 10.3847/1538-4357/ab86a9
  * Federrath (2015) Federrath, C. 2015, MNRAS, 450, 4035, doi: 10.1093/mnras/stv941
  * Federrath et al. (2014) Federrath, C., Schrön, M., Banerjee, R., & Klessen, R. S. 2014, ApJ, 790, 128, doi: 10.1088/0004-637X/790/2/128
  * Fluke & Jacobs (2020) Fluke, C. J., & Jacobs, C. 2020, Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery, 10, e1349
  * Frank et al. (2014) Frank, A., Ray, T. P., Cabrit, S., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 451, doi: 10.2458/azu_uapress_9780816531240-ch020
  * Furlan et al. (2016) Furlan, E., Fischer, W. J., Ali, B., et al. 2016, ApJS, 224, 5, doi: 10.3847/0067-0049/224/1/5
  * Guszejnov et al. (2021) Guszejnov, D., Grudić, M. Y., Hopkins, P. F., Offner, S. S. R., & Faucher-Giguère, C.-A. 2021, MNRAS, 502, 3646, doi: 10.1093/mnras/stab278
  * Gutermuth et al. (2009) Gutermuth, R. A., Megeath, S. T., Myers, P. C., et al. 2009, ApJS, 184, 18, doi: 10.1088/0067-0049/184/1/18
  * Gutermuth et al. (2011) Gutermuth, R. A., Pipher, J. L., Megeath, S. T., et al. 2011, ApJ, 739, 84, doi: 10.1088/0004-637X/739/2/84
  * Gutermuth et al. (in prep) Gutermuth et al. in prep
  * Hansen et al. (2012) Hansen, C. E., Klein, R. I., McKee, C. F., & Fisher, R. T. 2012, ApJ, 747, 22, doi: 10.1088/0004-637X/747/1/22
  * He et al. (2016) He, K., Zhang, X., Ren, S., & Sun, J. 2016, in Proceedings of the IEEE conference on computer vision and pattern recognition, 770–778
  * Heitsch et al. (2006) Heitsch, F., Slyz, A. D., Devriendt, J. E. G., Hartmann, L. W., & Burkert, A. 2006, ApJ, 648, 1052, doi: 10.1086/505931
  * Heyer et al. (2016) Heyer, M., Goldsmith, P. F., Yıldız, U. A., et al. 2016, MNRAS, 461, 3918, doi: 10.1093/mnras/stw1567
  * Heyer & Brunt (2012) Heyer, M. H., & Brunt, C. M. 2012, MNRAS, 420, 1562, doi: 10.1111/j.1365-2966.2011.20142.x
  * Ishii et al. (2019) Ishii, S., Nakamura, F., Shimajiri, Y., et al. 2019, PASJ, 71, S9, doi: 10.1093/pasj/psz065
  * Jofré et al. (2017) Jofré, P., Traven, G., Hawkins, K., et al. 2017, MNRAS, 472, 2517, doi: 10.1093/mnras/stx1877
  * Kirk et al. (2017) Kirk, H., Friesen, R. K., Pineda, J. E., et al. 2017, ApJ, 846, 144, doi: 10.3847/1538-4357/aa8631
  * Klessen & Hennebelle (2010) Klessen, R. S., & Hennebelle, P. 2010, A&A, 520, A17, doi: 10.1051/0004-6361/200913780
  * Kong et al. (2018) Kong, S., Arce, H. G., Feddersen, J. R., et al. 2018, ApJS, 236, 25, doi: 10.3847/1538-4365/aabafc
  * Li et al. (2015) Li, H., Li, D., Qian, L., et al. 2015, ApJS, 219, 20, doi: 10.1088/0067-0049/219/2/20
  * Li & Nakamura (2006) Li, Z.-Y., & Nakamura, F. 2006, ApJ, 640, L187, doi: 10.1086/503419
  * Lochner & Bassett (2020) Lochner, M., & Bassett, B. A. 2020, arXiv e-prints, arXiv:2010.11202.  https://arxiv.org/abs/2010.11202
  * Machida & Hosokawa (2013) Machida, M. N., & Hosokawa, T. 2013, MNRAS, 431, 1719, doi: 10.1093/mnras/stt291
  * Matzner (2007) Matzner, C. D. 2007, ApJ, 659, 1394, doi: 10.1086/512361
  * Megeath et al. (2012) Megeath, S. T., Gutermuth, R., Muzerolle, J., et al. 2012, AJ, 144, 192, doi: 10.1088/0004-6256/144/6/192
  * Moraghan et al. (2013) Moraghan, A., Kim, J., & Yoon, S. J. 2013, MNRAS, 432, L80, doi: 10.1093/mnrasl/slt044
  * Motte et al. (2014) Motte, F., Nguyên Luong, Q., Schneider, N., et al. 2014, A&A, 571, A32, doi: 10.1051/0004-6361/201323001
  * Nakamura & Li (2007) Nakamura, F., & Li, Z.-Y. 2007, ApJ, 662, 395, doi: 10.1086/517515
  * Nakamura et al. (2011) Nakamura, F., Kamada, Y., Kamazaki, T., et al. 2011, ApJ, 726, 46, doi: 10.1088/0004-637X/726/1/46
  * Nakamura et al. (2014) Nakamura, F., Sugitani, K., Tanaka, T., et al. 2014, ApJ, 791, L23, doi: 10.1088/2041-8205/791/2/L23
  * Nakamura et al. (2019) Nakamura, F., Ishii, S., Dobashi, K., et al. 2019, PASJ, 71, S3, doi: 10.1093/pasj/psz057
  * Narayanan et al. (2008) Narayanan, G., Heyer, M. H., Brunt, C., et al. 2008, ApJS, 177, 341, doi: 10.1086/587786
  * Narayanan et al. (2012) Narayanan, G., Snell, R., & Bemis, A. 2012, MNRAS, 425, 2641, doi: 10.1111/j.1365-2966.2012.21579.x
  * Offner & Arce (2014) Offner, S. S. R., & Arce, H. G. 2014, ApJ, 784, 61, doi: 10.1088/0004-637X/784/1/61
  * Offner & Chaban (2017) Offner, S. S. R., & Chaban, J. 2017, ApJ, 847, 104, doi: 10.3847/1538-4357/aa8996
  * Offner et al. (2012) Offner, S. S. R., Robitaille, T. P., Hansen, C. E., McKee, C. F., & Klein, R. I. 2012, ApJ, 753, 98, doi: 10.1088/0004-637X/753/2/98
  * Panopoulou et al. (2014) Panopoulou, G. V., Tassis, K., Goldsmith, P. F., & Heyer, M. H. 2014, MNRAS, 444, 2507, doi: 10.1093/mnras/stu1601
  * Pelletier & Pudritz (1992) Pelletier, G., & Pudritz, R. E. 1992, ApJ, 394, 117, doi: 10.1086/171565
  * Pokhrel et al. (2020) Pokhrel, R., Gutermuth, R. A., Betti, S. K., et al. 2020, ApJ, 896, 60, doi: 10.3847/1538-4357/ab92a2
  * Pokhrel et al. (2021) Pokhrel, R., Gutermuth, R. A., Krumholz, M. R., et al. 2021, ApJ, 912, L19, doi: 10.3847/2041-8213/abf564
  * Rebull et al. (2010) Rebull, L. M., Padgett, D. L., McCabe, C. E., et al. 2010, ApJS, 186, 259, doi: 10.1088/0067-0049/186/2/259
  * Reis et al. (2018) Reis, I., Poznanski, D., Baron, D., Zasowski, G., & Shahaf, S. 2018, MNRAS, 476, 2117, doi: 10.1093/mnras/sty348
  * Ridge et al. (2006) Ridge, N. A., Di Francesco, J., Kirk, H., et al. 2006, AJ, 131, 2921, doi: 10.1086/503704
  * Roccatagliata et al. (2020) Roccatagliata, V., Franciosini, E., Sacco, G. G., Randich, S., & Sicilia-Aguilar, A. 2020, A&A, 638, A85, doi: 10.1051/0004-6361/201936401
  * Ronneberger et al. (2015) Ronneberger, O., Fischer, P., & Brox, T. 2015, in International Conference on Medical image computing and computer-assisted intervention, Springer, 234–241
  * Shallue & Vanderburg (2018) Shallue, C. J., & Vanderburg, A. 2018, AJ, 155, 94, doi: 10.3847/1538-3881/aa9e09
  * Shimajiri et al. (2015a) Shimajiri, Y., Kitamura, Y., Nakamura, F., et al. 2015a, ApJS, 217, 7, doi: 10.1088/0067-0049/217/1/7
  * Shimajiri et al. (2015b) Shimajiri, Y., Sakai, T., Kitamura, Y., et al. 2015b, ApJS, 221, 31, doi: 10.1088/0067-0049/221/2/31
  * Shu et al. (1994) Shu, F., Najita, J., Ostriker, E., et al. 1994, ApJ, 429, 781, doi: 10.1086/174363
  * Shu et al. (1988) Shu, F. H., Lizano, S., Ruden, S. P., & Najita, J. 1988, ApJ, 328, L19, doi: 10.1086/185152
  * Tanabe et al. (2019) Tanabe, Y., Nakamura, F., Tsukagoshi, T., et al. 2019, PASJ, 71, S8, doi: 10.1093/pasj/psz100
  * Van Oort et al. (2019) Van Oort, C. M., Xu, D., Offner, S. S. R., & Gutermuth, R. A. 2019, ApJ, 880, 83, doi: 10.3847/1538-4357/ab275e
  * Wang et al. (2010) Wang, P., Li, Z.-Y., Abel, T., & Nakamura, F. 2010, ApJ, 709, 27, doi: 10.1088/0004-637X/709/1/27
  * Wolfire et al. (1995) Wolfire, M. G., Hollenbach, D., McKee, C. F., Tielens, A. G. G. M., & Bakes, E. L. O. 1995, ApJ, 443, 152, doi: 10.1086/175510
  * Xu et al. (2020a) Xu, D., Offner, S. S. R., Gutermuth, R., & Van Oort, C. 2020a, ApJ, 890, 64, doi: 10.3847/1538-4357/ab6607
  * Xu et al. (2020b) Xu, D., Offner, S. S. R., Gutermuth, R., & Van Oort, C. 2020b, ApJ, 905, 172, doi: 10.3847/1538-4357/abc7bf
  * Zhang et al. (2020) Zhang, S., Yang, J., Xu, Y., et al. 2020, ApJS, 248, 15, doi: 10.3847/1538-4365/ab879a