ALMA-IMF. IX. Catalog and Physical Properties of 315 SiO Outflow Candidates in 15 Massive Protoclusters

Starting from the fully processed cubes for spw1, we use the SpectralCube package ²⁴ to create subcubes for each field extending from V_LSR − 95 km s⁻¹ to V_LSR + 95 km s⁻¹, using a rest frequency of ν_SiO = 217.10498 GHz. This is the maximum velocity coverage common to all 15 fields in the sample based on the tuning for each target. We also created V_LSR ± 95 km s⁻¹ cutouts from the spw1 primary beam files produced by our imaging pipeline (see Ginsburg et al. 2022; Cunningham et al. 2023).

Our noise levels typically vary by 4%−8% between channels for a given field (columns 8 and 9 of Table 2), but can vary by as much as 11%–20% (G338.93, G351.77) or as little as 2.8% (W43-MM1). We take our noise in a single channel to be σ = 1.4826 × MAD. ²⁵ The noise is measured within an emission-free polygonal region near the center of each field of view; the same emission-free region is used for all channels in a given cube. The largest noise variations in our data appear to be caused by imaging artifacts.

In order to use the most accurate noise level for each outflow candidate, we used SpectralCube to create "noise cubes" whose values vary by both channel and pixel. We take the noise in each channel (σ_channel) and, at each pixel location, divide σ_channel by the value of the primary beam correction at that pixel. This has little effect on σ near the center of the image (where the primary beam response is ∼1) but will increase σ toward the edges of each mosaic. In this way, we created noise cubes in which the noise level varies with frequency according to σ measured for each channel, and varies spatially according to the effects of the primary beam.

Using these 3D noise cubes, we masked each image cube spectrally and spatially at the 3σ level. The resulting maps still contain some spurious emission, as expected for a 3σ cutoff given our typical ∼10⁹ pixels in a cube. In order to remove this emission, we used the scipy.ndimage package to perform 3D binary erosion (one iteration) and binary dilation (two iterations) on each mask. ²⁶ This procedure is equivalent to requiring that emission be present above the 3σ level in both spectral (≥3 consecutive channels) and spatial (≥3 pixels across) dimensions. The erosion/dilation step successfully removed nearly all of the spurious emission from our data cubes with minimal loss of true signal. Using these masked cubes, we created integrated-intensity (moment 0), intensity-weighted velocity (moment 1), intensity-weighted variance (moment 2), line width ($\sqrt{\mathrm{moment}\,2}$), and maximum-intensity maps for each field.

We performed an initial search for protostellar outflows using these maps and the unmasked line cubes. All examinations, including the initial inspection discussed above, were performed in the Cube Analysis and Rendering Tool for Astronomy (CARTA; Comrie et al. 2021). We identified candidates first by eye based on linear morphology in any map plus V − V_LSR > 5 km s⁻¹ in the moment 1 maps or line width > 10 km s⁻¹ in the line width maps. We then examined the line cubes directly to ensure that no regions with emission of V − V_LSR > 5 km s⁻¹ were missed in the moment-map examination. We do not require a continuum driving source to be positively identified in order to list a candidate, and we do not report specific driving sources for any outflows in this paper. However, the presence of a continuum source coincident with an outflow may increase our confidence that a particular candidate is indeed an outflow, or that a red- and blueshifted lobe share a driving source (see Section 3.3, below).

After initial identification, we drew a custom polygonal region around each candidate using its 3σ contours in the masked moment maps. We then modified (expanded) this region as appropriate based on by-eye examination of each channel in the line cube, to make sure no relevant emission was missed (e.g., faint or high-velocity). Each candidate's emission was then integrated spatially over this custom polygonal aperture in each channel to create an aperture-integrated spectrum, which was then examined by eye. We also generated a position–velocity (PV) diagram along each candidate's longest axis using the radio-astro-tools package PV Extractor. ²⁷ Finally, each candidate's overall structure was examined directly in the image cube channel by channel. Candidates were confirmed or rejected, and polygonal apertures modified, based on this final, spatial and kinematic examination.

For each outflow, we produce a summary figure showing its integrated intensity, velocity-weighted intensity, line width, and maximum-value maps, along with the aperture-integrated spectrum, PV diagram, and 1.3 mm continuum emission map. Examples are shown in Figures 1 and 2. The candidate shown in Figure 1 is a typical, symmetric bipolar outflow, while that shown in Figure 2 is bipolar but asymmetric, with a significantly larger, brighter, and higher-velocity redshifted lobe than blueshifted lobe. Detailed examination of individual candidates or fields is beyond the scope of this work.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (315 panels)

Tarball of figure set images

Zoom In Zoom Out Reset image size

Summary figures for all candidates can be viewed on Zenodo at doi:10.5281/zenodo.8350595. An initial list of outflow candidates was made available internally within the collaboration, with voluntary feedback requested. Following review and evaluation by 14 members, this initial list of candidates was refined.

In total, we detect 315 outflow candidates across 15 fields. Of the 315 candidates, 39 are classified as bipolar, 147 as blue monopolar, and 129 as red monopolar. ²⁸ We find a median of 17 candidates per field; the minimum number of candidates is three (G328.25), and the maximum is 47 (W51-IRS2). Our full catalog is presented in a machine-readable format; a representative example of the catalog is shown in Table 3. The full catalog can also be found in ESCV format on Zenodo at doi:10.5281/zenodo.8350595. For each outflow candidate, we report its approximate center position, its color, the total velocity range over which high-velocity emission is detected, the velocity at which the aperture-integrated spectrum peaks, the peak intensity of the aperture-integrated spectrum, the aperture- and velocity-integrated flux density, and our classification for that candidate.

The center positions are calculated as the center of the bounding box encompassing the polygonal aperture used for each outflow candidate. That is, the center R.A. for each candidate is the average of the minimum and maximum R.A. values in that candidate's polygonal aperture, and the center decl. is the average of the minimum and maximum decl. values. Candidates with colors listed as "red+blue" are bipolar. Candidates listed with just a single color ("red" or "blue") are either monopolar or, if potentially bipolar, a counterpart cannot be definitively identified from among multiple possibilities. The latter is especially common in regions with significant outflow activity and/or a high local number density of cores. Candidates are only listed as "red+blue" (bipolar) if the same 1.3 mm continuum driving source can be associated with both lobes with high confidence.

We identify velocity ranges by eye for each candidate based on its aperture-integrated spectrum, its PV diagram, and the unmasked line cube. In general, velocity ranges for each outflow exclude V_{LSR,candidate} ± 5 km s⁻¹ based on line shape, where V_{LSR,candidate} is the standard of rest velocity assessed locally for each candidate. We assess V_{LSR,candidate} locally because the clump V_LSR can vary across the field (see Table 1). In some rare cases, we exclude more or less of the velocity range around V_{LSR,candidate}, based on line shape in the integrated spectrum and channel-by-channel by-eye examination of candidate morphology. We derive physical properties for each candidate in Section 3.4. Although we restrict our search for outflows to velocities of V_LSR ± 95 km s⁻¹, the aperture-integrated spectra of five candidates (W43-MM2 Candidate #24; W51-E Candidates #9, #19, and #20; W51-IRS2 Candidate #16) show emission at even higher velocities. We use the full spw1 cubes for our analyses of these candidates only. For W43-MM1 Candidate #24, W51-E Candidate #20, and W51-IRS Candidate #16, the full outflow spectrum is cut off even in the spw1 cube, so the reported velocity ranges for those candidates should be considered lower limits.

We use the following three classifications: (1) likely, (2) possible, and (3) complex or cluster. "Likely" candidates are those we consider significantly likely to be protostellar outflows, based on their brightness, morphology, aperture-integrated spectrum, and structure in the PV diagram. Most of the candidates that have bright continuum sources in or near the polygonal aperture fall into this category. "Possible" candidates are those we consider likely or probable outflow candidates, but either their brightness, morphology, spectral structure, or PV structure is not quite definitive enough to place them in category #1. "Complex or cluster" candidates are those that clearly exhibit high-velocity emission but either (a) do not display the typical morphology of a protostellar outflow or (b) appear to be blended emission from multiple outflows and individual driving sources cannot be identified. In total, 129 candidates are classified as "likely," 180 are classified as "possible," and six are classified as "complex or cluster."

Outflow activity and outflow-core associations will be explored separately using a multispecies approach for certain individual targets in future publications (e.g., Armante et al. 2023 for G012.80/W33, and M. Valeille-Manet et al. 2024, in preparation, for all massive cores in the 15 ALMA-IMF regions).

We stress that each outflow candidate remains an outflow candidate. We also note that many of our candidates do not have clear 1.3 mm continuum peaks within the outflow path or immediately adjacent to one end or the other of the longest axis, i.e., do not have candidate driving sources. In these cases, we suggest that either the SiO emission is tracing only the leading edge of the flow, or that this is an observational limitation driven by our sensitivity in each field. Both possibilities are consistent with the picture of protostellar outflows in which SiO preferentially traces the (smaller) active shocks and warm/hot gas in outflows and lower-J CO transitions trace the (larger) coolest, outer layers of an outflow's cavity walls (Bally 2016).

The total number of bipolar candidates in our sample should be considered a lower limit. The low overall fraction of bipolar candidates is a result of the interactions between our identification method, classification criteria, and the highly clustered nature of the dust-core populations in most fields. These limitations (single-species data, clustered dust cores) also affect our ability to positively identify continuum driving sources for a number of our candidates. There are many fields (especially but not only G012.80, G333.60, W43-MM1, W51-E, and W51-IRS2) in which numerous red- and blueshifted outflows appear to be emanating from a single location in which two or more (usually >4) closely packed dust cores are also located. In these cases, even slight misalignments in position angle between outflow candidates introduce significant uncertainty as to which candidates are associated with which cores, and thus with each other. We elected to leave most candidates in these confused regions classified as red or blue monopolar, rather than bipolar, because we are confident in their candidacy as outflows but less confident as to their specific red/blue pairings.

Some targets have overlapping fields of view (W43-MM2 and W43-MM3; W51-E and W51-IRS2). Consequently, four outflow candidates are detected in more than one field: W43-MM3 candidate #1 and candidate #2 are also detected at the eastern edge of the W43-MM2 field of view, and W51-IRS2 candidate #46 and candidate #47 are also detected in the northwestern quadrant of the W51-E field of view. We analyze these candidates using the W43-MM3 and W51-IRS2 data cubes, respectively, as these cubes had higher signal-to-noise ratios at the locations of these candidates. These candidates appear in the catalog under W43-MM3 and W51-IRS2 only, i.e., the entries are not duplicated under the alternate fields.

There were several findings in our SiO data set that are scientifically interesting, but beyond the scope of this catalog paper. These findings are briefly described in Appendix A, and will be explored in future publications.

3.3.1. Flux Filtering on Large Spatial Scales

3.3.2. Comparison to Nony et al. (2023)

3.3.3. Crossing Outflows

3.3.4. Outflow Inclination Angle

For each outflow candidate, we derive its median SiO column density and its total mass, momentum, kinetic energy, mass rate, momentum rate, and energy rate. These properties are presented for each candidate in a machine-readable table; a representative example is shown in Table C1.

We first convert each cube from janskys per beam to kelvin, and then extract the spatially integrated spectrum of each candidate using the custom polygonal apertures described in Section 3.2. These aperture-integrated spectra are the basis of our derivation of the physical properties of each candidate. We use a channel-based calculation method (see, e.g., Maud et al. 2015) in which each physical property is calculated separately for each channel and then summed, rather than using velocity-integrated fluxes. This channel-based method has been shown to reduce the overestimation of outflow momenta and kinetic energies that can occur when aperture-integrated intensities are multiplied by an outflow's maximum velocity only (see Maud et al. 2015, and references therein).

3.4.1. Derivation of Column Density

After converting our cube to units of kelvin and extracting the aperture-integrated spectrum for a candidate, we mask out any channels in which the aperture-integrated brightness temperature is <3σ, where σ is the aperture-integrated noise level. This masking step helps to prevent high-velocity, low-signal features from disproportionately impacting both the derived momenta and energies and their associated uncertainties at later stages.

To derive column density in each channel individually, we adopt a discrete form of the general equation for molecular column density in the optically thin approximation (see Appendix B):

$$\begin{eqnarray}&&N(i)=\displaystyle \frac{8\pi {k}_{{\rm{B}}}{\nu }^{2}}{{{hc}}^{3}{A}_{\mathrm{ul}}}\displaystyle \frac{{Q}_{\mathrm{rot}}}{{g}_{J}{g}_{K}{g}_{I}}\exp \left(\displaystyle \frac{{E}_{u}}{{k}_{{\rm{B}}}{T}_{\mathrm{ex}}}\right){T}_{{\rm{B}},i}{\rm{\Delta }}v\end{eqnarray} \tag{ 1 }$$

where k_B is the Boltzmann constant, ν is the rest frequency of the SiO J = 5–4 transition, h is the Planck constant, c is the speed of light, A_ul is the Einstein coefficient of spontaneous emission from the upper state to the lower state, Q_rot is the partition function of the SiO molecule, g_J , g_K , and g_I together represent the total degeneracy of the rotational state, E_u is the energy of the upper state of the transition, T_ex is the excitation temperature, T_B,i is the aperture-integrated brightness temperature in channel i, and Δv is the channel width. N(i) is the aperture-integrated column density in channel i. We assume that the SiO emission is optically thin, and assume an excitation temperature of T = 50${}_{-20}^{+30}$ K for all candidates. A detailed discussion of our reasons for these choices can be found in Appendix B.

The nonlinear relationship between T_ex and N_SiO leads to asymmetric uncertainties in our column densities, and in all properties that are subsequently derived from them (mass, momentum, energy, and associated rates). We propagate this asymmetric uncertainty by calculating two Gaussian uncertainties for each N(i): one using ${\sigma }_{{T}_{\mathrm{ex},\mathrm{lower}}}$ = −20 K, and one using ${\sigma }_{{T}_{\mathrm{ex},\mathrm{upper}}}$ = +30 K. The total error distribution for each derived N(i) becomes the combination of two half-Gaussians: below the calculated N(i) value, it is a Gaussian with σ = σ_lower and centered at N(i), and above the calculated N(i), it is a Gaussian with σ = σ_upper centered at N(i). When summing the data either spatially or spectrally, the lower- and upper-bound uncertainties are summed in quadrature separately, i.e.,

$$\begin{eqnarray}&&{\sigma }_{N\,\mathrm{lower},\mathrm{sum}}={\left(\displaystyle \sum _{i}{({\sigma }_{N\,\mathrm{lower},i})}^{2}\right)}^{1/2}\end{eqnarray} \tag{ 2 }$$

and likewise for σ_upper. We find that our column density uncertainties can be dominated either by ${\sigma }_{{T}_{\mathrm{ex}}}$ or by the inherent noise in the data itself, depending on the data noise level. Fields with higher median σ (see Table 2) tend to be noise-dominated, and those with lower median σ tend to be dominated by the uncertainty in T_ex.

3.4.2. Derived Masses, Momenta, and Kinetic Energies

The total gas mass in each channel can be calculated following

$$\begin{eqnarray}&&{M}_{\mathrm{gas}}(i)=N(i){\mu }_{g}{m}_{{{\rm{H}}}_{2}}{\rm{\Omega }}{D}^{2}/\left[\displaystyle \frac{\mathrm{SiO}}{{{\rm{H}}}_{2}}\right]\end{eqnarray} \tag{ 3 }$$

where μ_g = 1.36 is the total gas mass relative to H₂. ${m}_{{{\rm{H}}}_{2}}$ is the mass of a hydrogen molecule, Ω is the solid angle subtended by a single pixel in our image cubes, D is the distance to the source, N(i) is given by Equation (1), $\left[\tfrac{\mathrm{SiO}}{{{\rm{H}}}_{2}}\right]$ is the fractional abundance of SiO relative to molecular hydrogen.

We adopt a flat SiO-to-H₂ abundance ratio of 10^−8.5 (or, 3.16 × 10⁻⁹) for all outflow candidates in our sample, taking into consideration the wide range of abundance ratios in astrochemical shock models (Schilke et al. 1997; Gusdorf et al. 2008), intraoutflow abundance variations (Bally 2016), and typical abundance ratio values reported in the literature (see, e.g., Codella et al. 1999; Lu et al. 2021). A detailed discussion of our choice of $\left[\tfrac{\mathrm{SiO}}{{{\rm{H}}}_{2}}\right]$ can be found in Appendix B. In short, these theoretical and observational studies have shown abundance ratios ranging from 10⁻¹¹ to 10⁻⁶ both within individual outflows and sometimes between outflows. Such large variations are not guaranteed to occur within any given outflow, but they are possible across a population. This means that our adopted ratio of 10^−8.5 could potentially vary by up to 2 orders of magnitude from the true abundance at any single location within an individual outflow candidate.

However, because we are unable to measure the SiO abundance ratio directly (see Appendix B, Section 3), we also do not have a measurement uncertainty for our assumed abundance ratio. For the purposes of error propagation, we therefore treat the fractional abundance as definitional (i.e., σ = 0). This assumed fractional abundance may be an overestimate for some sources and an underestimate for others. This would increase the scatter in our data at the population level, but it is unlikely to change underlying fundamental correlations or distributions unless there is a trend in this over- or underestimation. We compare our data to the SiO-derived outflow masses in the literature (including Lu et al. 2021) in Section 4.1 and find no evidence for such a trend.

With the adoption of $\left[\tfrac{\mathrm{SiO}}{{{\rm{H}}}_{2}}\right]$ = 10^−8.5, Equation (3) becomes

$$\begin{eqnarray}&&{M}_{\mathrm{gas}}(i)=\left(1\times {10}^{8.5}\right)N(i){\mu }_{g}{m}_{{{\rm{H}}}_{2}}{\rm{\Omega }}{D}^{2}\end{eqnarray} \tag{ 4 }$$

alternately, substituting Equation (1) for N(i) gives us

$$\begin{eqnarray}&&{M}_{\mathrm{gas}}(i)=C\displaystyle \sum _{i}{T}_{{\rm{B}}}(i)\end{eqnarray} \tag{ 5 }$$

where T_B(i) is the aperture-integrated brightness temperature of a single channel and

$$\begin{eqnarray}&&C=2.08\times {10}^{-8}{\rm{\Delta }}v{\rm{\Omega }}{D}^{2}\end{eqnarray} \tag{ 6 }$$

is the constant of proportionality between brightness temperature and mass for SiO. In Equation (6), the leading numerical factor has units of [g s cm⁻³ K⁻¹], Δv is in centimeters per second, Ω is in steradians, and D is in centimeters.

The total, velocity- and spatially integrated mass is then

$$\begin{eqnarray}&&{M}_{\mathrm{out}}=\displaystyle \sum _{i}{M}_{\mathrm{gas}}(i)\end{eqnarray} \tag{ 7 }$$

where the summation is performed over all channels in the outflow's unique velocity range (see the complete version of Table 3 for velocity ranges).

We then calculate outflow momentum as

$$\begin{eqnarray}&&{P}_{\mathrm{out}}=\displaystyle \sum _{i}M(i){v}_{i}\end{eqnarray} \tag{ 8 }$$

and outflow kinetic energy as

$$\begin{eqnarray}&&{E}_{\mathrm{out}}=\displaystyle \frac{1}{2}\displaystyle \sum _{i}M(i){v}_{i}^{2}\end{eqnarray} \tag{ 9 }$$

where v_i is the velocity of channel i relative to local V_LSR in both Equations (8) and (9). As discussed in Section 3.3.4, we do not assume an inclination angle when deriving these properties.

Figure 3 shows the log-space distribution of these properties for the full sample (all 315 candidates, with the red and blue lobes of bipolar candidates counted separately for a total of 355 data points in the plotted bins). All three panels show stacked histograms.

Zoom In Zoom Out Reset image size

In general, the distributions have well-defined peaks but are broad, spanning >3.5 orders of magnitude in mass, >4 orders of magnitude in momentum, and >5 orders of magnitude in kinetic energy. Minimum, maximum, and median values for mass, momentum, and energy for the full sample are listed in the upper section of Table 4. These values are not adjusted for inclination angle.

Table 4. Full Sample Mass, Momentum, and Energy Statistics

Property	Min	Max	Median a
Mblue [M⊙]	0.003	10.0	0.25 (0.15)
Mred [M⊙]	0.003	14.1	0.34 (0.23)
Mtot [M⊙]	0.005	14.1	0.30 (0.20)
Pblue [M⊙ km s−1]	0.019	402	3.4 (2.4)
Pred [M⊙ km s−1]	0.030	538	4.7 (3.5)
Ptot [M⊙ km s−1]	0.031	538	4.1 (2.9)
Eblue [erg]	1.1 × 1042	2.3 × 1047	5.8 × 1044 (4.8)
Ered [erg]	1.9 × 1042	3.4 × 1047	8.3 × 1044 (7.1)
Etot [erg]	1.9 × 1042	3.4 × 1047	6.8 × 1044 (5.6)
Inclination-adjusted (i = 573)
Pblue/cos i	0.035	744	6.2 (4.3)
Pred/cos i	0.056	996	8.7 (6.5)
Ptot/cos i	0.057	996	7.6 (5.3)
Eblue/cos2 i	3.8 × 1042	8.0 × 1047	2.0 × 1045 (1.6)
Ered/cos2 i	6.5 × 1042	1.2 × 1048	2.8 × 1045(2.4)
Etot/cos2 i	6.5 × 1042	1.2 × 1048	2.3 × 1045 (1.9)

Notes.

^a The min, max, and median values reported in this table are calculated from the full sample of 315 candidates. ^b Uncertainties on the medians are the scaled MAD, listed in parentheses. For values listed in scientific notation, the order of magnitude of the uncertainty is the same as that of the reported median. ^c The values in the upper section of the table are not adjusted for inclination angle. The values in the lower section assume a uniform inclination angle of 573 for all candidates.

Download table as:  ASCIITypeset image

We derive inclination-adjusted mass, momentum, and energy statistics for the full sample assuming a uniform inclination angle of ∼573, following the method of Bontemps et al. (1996). These statistics are listed in the lower section of Table 4, as are the inclination correction factors for each property. The inclination-corrected statistics are not used in our analysis unless specifically noted.

Table C1 shows our derived median column density, mass, momentum, and kinetic energy for the red- and blueshifted components of each outflow candidate. The first 10 lines are shown here. The full, machine-readable version of this table, including uncertainties for each value, is available and in ESCV format on Zenodo at doi:10.5281/zenodo.8350595.

Overall, we find no statistically significant difference between mass, momentum, or energy values for the red versus blue outflow lobes; this is consistent with a lack of any strong detection bias toward strong or weak emission with lobe color. Our high energy maxima can be attributed to those outflows that have both bright emission at all velocities and strong emission at ∣V − V_LSR∣ ≥60 km s⁻¹. Because energy goes as v², emission at high velocities has an outsized effect on the total derived energy. In most cases, this high-velocity gas is all part of the outflow, but in a small subset of cases, this "high-velocity" emission is due to hot-core line emission contaminating the outflow aperture. This appears as "high-velocity" emission because the lines are at different rest frequencies from SiO 5–4, and so this hot-core line contamination has a strong effect on the derived energies in these few cases. There are 12 outflow candidates within the sample with significant contamination from hot-core lines: G351.77 Candidate #3, W43-MM1 Candidates #16, #17, and #27, W43-MM2 Candidates #14 and #15, W51-E Candidates #19 and #20, and W51-IRS2 Candidates #10, #28, #38, and #40. The derived properties of these candidates should therefore be considered upper limits. The hot cores in each field are explored separately in Brouillet et al. (2022) and Bonfand et al. (2023).

3.4.3. Derived Dynamical Times and Mass, Momentum, and Energy Rates

In order to determine mass flow rate $\dot{M}$, momentum supply rate $\dot{P}$ (alternately mechanical force, F_m ), and outflow power $\dot{E}$ (alternately mechanical luminosity L_m ), most teams measure the distance between an outflow and its driving source and divide this by outflow velocity in order to determine a rough outflow dynamical time. This approach requires the identification of a continuum driving source for each candidate. Since we do not (and in many cases cannot) assign our candidates to specific driving sources in this work, we cannot use this approach. Instead, we use the PV path length (see Figures 1 and 2) as a proxy for outflow size. Our PV path length is the same in all channels, so a channel-by-channel calculation of outflow dynamical time is not possible with this approach. Instead, we divide the outflow path length by the median relative velocity (median intensity-weighted velocity minus V_LSR) to determine a single t_dyn for each candidate.

We assign an uncertainty of ±15% to all outflow dynamical times. This uncertainty is a consequence of our use of path length as a proxy for outflow size. When creating the PV paths, we sometimes extend the path beyond the end of the outflow in order to include baseline zero-emission regions in the PV diagram. Likewise, in very crowded regions, path lengths are truncated slightly to avoid confusion with nearby features. In all cases, the difference between the path marked in CARTA as the "true" outflow size (identified by eye in the moment maps) is ≲15%. Therefore, this is the uncertainty we adopt for the path length and outflow dynamical times.

Our median dynamical time from this method is 6000 ± 2800 yr (where the uncertainty is the scaled MAD), and our minimum and maximum dynamical times are 510 yr and 70,000 yr, respectively. Outflow mass rate is then M_out/t_dyn, outflow momentum rate is P_out/t_dyn, and outflow energy rate is E_out/t_dyn.

We do not calculate a dynamical time for any candidates classified as "complex or cluster," as for these candidates the path is arbitrary and does not reflect a specific outflow axis. This criterion eliminates six candidates from further analysis. We also do not calculate dynamical time or associated rates for any candidates whose path length is less than twice the length of the beam major axis, i.e., candidates whose largest axis remains unresolved. Typically, these are candidates suspected of having a face-on orientation. This criterion eliminates five more candidates from further analysis. In total, we reduce our total number of candidates to 304 for the analysis of mass, momentum, and energy rates and all associated figures.

Figure 4 shows the log-space distribution of these rates for the full sample. All three panels show stacked histograms. Table 5 shows the minimum, maximum, and median values for each derived rate. The values in the upper section of the Table 5 are not adjusted for inclination angle. Inclination-adjusted values are listed in the lower section of Table 5, as are the inclination correction factors for each property. We assume a uniform inclination angle of 57° for all candidates to derive these values. The inclination-corrected values are not used in our analysis unless specifically noted.

Zoom In Zoom Out Reset image size

Table 5. Full Sample Mass, Momentum, and Energy Rate Statistics

Property a	Min	Max	Median b
tdyn [yr]	510	70,000	6000 (2800)
${\dot{M}}_{\mathrm{blue}}$[M⊙ yr−1]	3.4 × 10−7	0.006	3.9 × 10−5 (2.7)
${\dot{M}}_{\mathrm{red}}$ [M⊙ yr−1]	6.0 × 10−7	0.003	4.7 × 10−5 (3.3)
${\dot{M}}_{\mathrm{tot}}$ [M⊙ yr−1]	1.0 × 10−6	0.006	5.0 × 10−5 (3.2)
${\dot{P}}_{\mathrm{blue}}$ [M⊙ km s−1 yr−1]	1.9 × 10−6	0.22	6.1 × 10−4 (4.9)
${\dot{P}}_{\mathrm{red}}$ [M⊙ km s−1 yr−1]	6.0 × 10−6	0.17	8.0 × 10−4 (6.1)
${\dot{P}}_{\mathrm{tot}}$ [M⊙km s−1 yr−1]	6.0 × 10−6	0.22	7.0 × 10−4 (5.7)
${\dot{E}}_{\mathrm{blue}}$ [L⊙]	0.001	1100	0.9 (0.8)
${\dot{E}}_{\mathrm{red}}$ [L⊙]	0.003	900	1.1 (1.0)
${\dot{E}}_{\mathrm{tot}}$ [L⊙]	0.003	1100	1.1 (1.0)
Inclination-adjusted c (i = 573)
tdyn/tan i	330	45,000	3900 (1800)
${\dot{M}}_{\mathrm{blue}}$ (tan i)	5.3 × 10−7	0.009	6.1 × 10−5 (4.2)
${\dot{M}}_{\mathrm{red}}$ (tan i)	9.3 × 10−7	0.005	7.3 × 10−5 (5.1)
${\dot{M}}_{\mathrm{tot}}$ (tan i)	1.6 × 10−6	0.009	7.8 × 10−5 (5.0)
${\dot{P}}_{\mathrm{blue}}$ (sin i / cos2 i)	5.5 × 10−6	0.63	0.0018 (0.0014)
${\dot{P}}_{\mathrm{red}}$ (sin i / cos2 i)	1.7 × 10−5	0.49	0.0023(0.0018)
${\dot{P}}_{\mathrm{tot}}$ (sin i / cos2 i)	1.7 × 10−5	0.55	0.0020 (0.0016)
${\dot{E}}_{\mathrm{blue}}$ (sin i / cos3 i)	0.005	5900	4.8 (4.3)
${\dot{E}}_{\mathrm{red}}$ (sin i / cos3 i)	0.018	4800	5.9 (5.3)
${\dot{E}}_{\mathrm{tot}}$ (sin i / cos3 i)	0.018	5900	5.9 (5.3)

Notes.

^a The min, max, and median values reported in this table are calculated from the full sample of 315 candidates. ^b Uncertainties on the medians are the scaled MAD, listed in parentheses. For values listed in scientific notation, the order of magnitude of the uncertainty is the same as that of the reported median. ^c The values in the upper section of the table are not adjusted for inclination angle. The values in the lower section assume a uniform inclination angle of 573 for all candidates.

Download table as:  ASCIITypeset image

As in Section 3.4.2, we find no significant differences between the rates derived for blueshifted outflow lobes and those derived for redshifted lobes. Likewise, we find that each rate has a reasonably well-defined peak but that the distributions are again broad, spanning 3.7–4.2 orders of magnitude in mass rate, 4.6–5 orders of magnitude in momentum rate, and 5.5–6 orders of magnitude in kinetic energy rate.

In this section, we compare the physical properties derived for our sample to those derived in the literature for similar high-mass samples. Specifically, we compare to lower-resolution, single-dish SiO (Csengeri et al. 2016) and CO (Maud et al. 2015) observations of larger massive protocluster samples than ours, and to a smaller sample of protoclusters detected in SiO with similar spatial resolution to our data (Lu et al. 2021). This comparison allows us to test both methodological effects and the impact of interferometric versus single-dish data, as well as placing our findings in broader context with the literature. We find good agreement between our typical (median, minimum/maximum) outflow derived properties and those reported in the literature for these other samples, suggesting that any methodological or observational bias in our candidate identification procedure or derivation of physical properties has not had a significant impact on our results. We describe our comparison to each sample in greater detail in the following paragraphs.

We first compare our derived outflow column densities to those derived by Csengeri et al. (2016) for a sample of massive clumps selected from the ATLASGAL survey. Csengeri et al. (2016) observed 430 sources with the IRAM 30 m telescope at 84 ∼ 115 GHz (∼26'' beam), and a subsample of 128 sources with the APEX telescope at 217 GHz (29'' beam). For their full sample, Csengeri et al. (2016) derived column densities of 1.6 × 10¹²–7.9 × 10¹³ cm⁻² using SiO J = 2−1 data and assuming an LTE approximation. For their subsample of 128 sources measured with both SiO J = 2−1 and SiO J = 5–4, Csengeri et al. (2016) found a column density range of 9.6 × 10¹¹–1.4 × 10¹⁴ cm⁻² using RADEX modeling, depending on the treatment of the beam filling factor. Our median column densities for any individual field range from 5.0 × 10¹³ cm⁻² to 3.5 × 10¹⁴ cm⁻², and the median column density across all 315 candidates ranges from 9 × 10¹² to 1.2 × 10¹⁵ cm⁻², with a median of 1 × 10¹⁴ cm⁻². These values are reported in Table 6. Our SiO column densities overlap with those of these ATLASGAL sample, but trend ∼1 order of magnitude higher overall. This trend is likely due to differences in angular resolution, as the Csengeri et al. (2016) single-dish data are more strongly affected by beam dilution than our interferometric data.

Table 6. Field-aggregated Outflow Properties

Field	NSiO,median a	Mout b	Pout	Eout	${\dot{M}}_{\mathrm{out}}$	${\dot{P}}_{\mathrm{out}}$	${\dot{E}}_{\mathrm{out}}$	M ${}_{\mathrm{cores},\mathrm{iso}}$ c
	(1013 cm−2)	(M⊙)	(M⊙ km s−1)	(1046 erg)	(10−4 M⊙ yr−1)	(10−2 M⊙ km s−1 yr−1)	(L⊙)	(M⊙)
G008.67	9.5	${2.52}_{-0.04}^{+0.03}$	${52.4}_{-0.8}^{+0.6}$	${1.62}_{-0.02}^{+0.02}$	${2.4}_{-0.1}^{+0.1}$	${0.5}_{-0.03}^{+0.03}$	${14}_{-1}^{+1}$	127 (4)
G010.62	6.5	${8.74}_{-0.1}^{+0.09}$	${141}_{-1}^{+1}$	${4.26}_{-0.06}^{+0.05}$	${7.0}_{-0.3}^{+0.3}$	${1.42}_{-0.08}^{+0.08}$	${45}_{-4}^{+4}$	202 (5)
G012.80	7.0	${9.56}_{-0.08}^{+0.1}$	${184}_{-1}^{+2}$	${5.81}_{-0.06}^{+0.05}$	${16.2}_{-0.8}^{+0.8}$	${4.0}_{-0.2}^{+0.2}$	${128}_{-8}^{+8}$	278 (5)
G327.29	10.5	${7.54}_{-0.09}^{+0.07}$	${125}_{-1}^{+1}$	${2.79}_{-0.03}^{+0.03}$	${19.6}_{-0.8}^{+0.7}$	${3.4}_{-0.2}^{+0.1}$	${67}_{-4}^{+4}$	1525 (5)
G328.25	20.0	${1.40}_{-0.04}^{+0.04}$	${31.8}_{-0.9}^{+0.8}$	${1.03}_{-0.03}^{+0.03}$	${3.4}_{-0.4}^{+0.4}$	${0.8}_{-0.1}^{+0.1}$	${22}_{-3}^{+3}$	130 (1)
G333.60	5.0	${7.5}_{-0.1}^{+0.1}$	${128}_{-1}^{+2}$	${3.43}_{-0.04}^{+0.05}$	${15.0}_{-0.7}^{+0.7}$	${3.2}_{-0.2}^{+0.2}$	${85}_{-6}^{+6}$	448 (5)
G337.92	6.0	${8.7}_{-0.2}^{+0.2}$	${117}_{-2}^{+2}$	${2.30}_{-0.04}^{+0.04}$	${9.4}_{-0.5}^{+0.5}$	${1.68}_{-0.08}^{+0.08}$	${37}_{-2}^{+2}$	507 (9)
G338.93	10.0	${11.5}_{-0.2}^{+0.2}$	${215}_{-5}^{+5}$	${5.7}_{-0.1}^{+0.1}$	${10.5}_{-0.5}^{+0.5}$	${2.1}_{-0.1}^{+0.1}$	${51}_{-3}^{+3}$	512 (3)
G351.77	25.0	${11.3}_{-0.5}^{+0.6}$	${239}_{-8}^{+10}$	${8.3}_{-0.3}^{+0.4}$	${31}_{-3}^{+3}$	${6.4}_{-0.5}^{+0.5}$	${170}_{-10}^{+10}$	515 (6)
G353.41	7.0	${4.2}_{-0.1}^{+0.1}$	${104}_{-3}^{+3}$	${3.5}_{-0.1}^{+0.1}$	${10.3}_{-0.5}^{+0.5}$	${2.7}_{-0.2}^{+0.2}$	${82}_{-5}^{+5}$	172 (3)
W43−MM1	10.5	${44.6}_{-0.3}^{+0.3}$	${1039}_{-5}^{+9}$	${36.0}_{-0.2}^{+0.3}$	${89}_{-3}^{+3}$	${22.4}_{-0.9}^{+0.9}$	${680}_{-30}^{+30}$	1683 (12)
W43−MM2	10.5	${11.52}_{-0.1}^{+0.08}$	${238}_{-2}^{+1}$	${7.83}_{-0.08}^{+0.06}$	${15.3}_{-0.6}^{+0.5}$	${3.5}_{-0.1}^{+0.1}$	${100}_{-6}^{+6}$	582 (4)
W43−MM3	5.5	${4.04}_{-0.06}^{+0.05}$	${51.6}_{-0.7}^{+0.6}$	${0.9}_{-0.01}^{+0.01}$	${2.6}_{-0.1}^{+0.1}$	${0.33}_{-0.02}^{+0.02}$	${4.3}_{-0.3}^{+0.3}$	170 (3)
W51−E	35.0	${56.6}_{-0.4}^{+0.9}$	${1550}_{-10}^{+20}$	${73.3}_{-1.1}^{+0.9}$	${130}_{-10}^{+10}$	${46}_{-5}^{+5}$	${2100}_{-200}^{+200}$	2883 (35)
W51−IRS2	7.5	${46.3}_{-0.4}^{+0.5}$	${990}_{-10}^{+10}$	${33.8}_{-0.4}^{+0.3}$	${79}_{-5}^{+5}$	${17}_{-1}^{+1}$	${500}_{-40}^{+40}$	2473 (20)

Notes.

^a The median SiO column density across all outflow candidates in the given field. Medians are calculated only from pixels meeting the significance threshold within each aperture and channel. ^b The field-total M_out (and P_out, E_out, etc.) is the sum of the derived mass of each individual outflow candidate in the given field. Upper (lower) uncertainties are the square root of the quadrature sum of upper (lower) uncertainties for each individual candidate. ^c The total mass in cores in each field, derived using the flux-density values for each core listed in Appendix D of Louvet et al. (2023) and assuming T = 15 K, τ ≪ 1, and h ν ≪ kT. Uncertainties are the square root of the quadrature sum of the uncertainties of the individual cores.

Download table as:  ASCIITypeset image

We compare our derived outflow masses to those derived by Lu et al. (2021) for their sample of massive star-forming regions in the central molecular zone (CMZ). Lu et al. (2021) targeted 834 molecular clumps with 2000 au resolution and detect 43 outflows. They derived a separate mass for each outflow using each of their molecular tracers; their SiO-derived outflow masses range from a few hundredths to a few tens of solar masses. The masses we derive for our outflow candidates largely fall within this range (see Table 4), but our minimum derived masses extend ∼1 order of magnitude lower than those of Lu et al. (2021). This can likely be attributed to the greater distance to the CMZ (>8 kpc), which will result in decreased mass sensitivity. This consistency is notable considering Lu et al. (2021) derived position-dependent SiO abundance ratios for each outflow. The overall agreement between our mass range and theirs suggests that our choice of $\left[\tfrac{\mathrm{SiO}}{{{\rm{H}}}_{2}}\right]$= 10^−8.5 is a reasonable first-order approximation of SiO abundance at the population level. Though specific abundances may vary within or between individual outflows, the true average value in our data appears to be well-represented by 10^−8.5 to first order.

Lu et al. (2021) found typical outflow velocities of several tens of kilometers per second and an overall range of a few kilometers per second to > 90 km s⁻¹ for their CMZ sample (see Lu et al. 2021), comparable to those we derive for our candidates (see Table 3).

We also compare our outflow properties to those of Maud et al. (2015), who used the JCMT to examine ¹²CO and ¹³CO J = 3−2 emission toward 99 massive young stellar objects drawn from the Red MSX Source survey (RMS). For each outflow, Maud et al. (2015) derived its mass, momentum, kinetic energy, dynamical time, and mass, momentum, and energy rates. As theirs are single-dish data, Maud et al. (2015) only reported a single red and blue outflow for each field.

In order to compare our results to those of Maud et al. (2015), we created "field-aggregated" outflow properties for each protocluster by summing each property (mass, momentum, etc.) across all outflows in each field. These values are reported in Table 6. We find a median field-aggregated outflow mass of ${8.74}_{-0.1}^{+0.09}$ M_⊙, with a minimum of ${1.40}_{-0.04}^{+0.04}$ and a maximum of ${56.6}_{-0.4}^{+0.9}$. We find median, minimum, and maximum field-aggregated outflow momenta of ${141}_{-1}^{+1}$ M_⊙ km s⁻¹, ${31.8}_{-0.9}^{+0.8}$ M_⊙ km s⁻¹, and ${1550}_{-10}^{+20}$ M_⊙ km s⁻¹, and median, minimum, and maximum field-aggregated kinetic energies of ${4.26}_{-0.06}^{+0.05}$ × 10⁴⁶ erg, ${0.9}_{-0.01}^{+0.01}$ × 10⁴⁶ erg, and ${73.3}_{-1.1}^{+0.9}$ × 10⁴⁶ erg, respectively.

We find that our field-aggregated outflow properties typically fall within the ranges observed by Maud et al. (2015), who found outflow masses of ∼0.7 M_⊙–1000 M_⊙, outflow momenta of ∼3 M_⊙ km s⁻¹− 4000 M_⊙ km s⁻¹, and outflow kinetic energies of ∼10⁴⁴ erg–3 × 10⁴⁷ erg for their outflows. We do note that our derived total outflow masses tend toward the lower end of the distribution observed by Maud et al. (2015) for their sample, our momenta are largely in line with the RMS-derived distribution, and our energies tend toward the higher end of the RMS-derived distribution. These trends can be explained by the difference in molecular tracers used and in the difference between interferometric and single-dish angular resolution. CO is more abundant and widespread than SiO and has a longer gas-phase lifetime, so Maud et al. (2015) likely have greater mass sensitivity for their sample than we do for ours. However, their CO-derived outflow velocity ranges are typically narrower than those we observe with SiO by factors of ∼1.5, while we have numerous small, high-velocity bullets that are more easily detected with SiO and interferometric observations. A decreased mass sensitivity but increased sensitivity to higher-velocity gas in our data would explain these comparative trends in both mass (which is velocity-independent) and momentum and energy (which are linearly and quadratically dependent on velocity, respectively).

We find shorter median dynamical times for our candidates (t_dyn = 6000 ± 2800 yr) as compared to Maud et al. (2015; 65,000 ± 34,000 yr). We find median, minimum, and maximum field-aggregated mass flow rates ($\dot{M}$) of ${15.0}_{-0.7}^{+0.7}$ × 10⁻⁴ M_⊙ yr⁻¹, ${2.4}_{-0.1}^{+0.1}$ × 10⁻⁴ M_⊙ yr⁻¹, and ${130}_{-10}^{+10}$ × 10⁻⁴ M_⊙ yr⁻¹. Our median, minimum, and maximum field-aggregated mechanical force values ($\dot{P}$) are ${3.2}_{-0.02}^{+0.02}$ × 10⁻² M_⊙ km s⁻¹ yr⁻¹, ${0.33}_{-0.02}^{+0.02}$ × 10⁻² M_⊙ km s⁻¹ yr⁻¹, and ${46}_{-5}^{+5}$ × 10⁻² M_⊙ km s⁻¹ yr⁻¹. Our median, minimum, and maximum field-aggregated kinetic energy rates ($\dot{L}$) are ${82}_{-5}^{+5}$ L_⊙, ${4.3}_{-0.3}^{+0.3}$ L_⊙, and ${2100}_{-200}^{+200}$ L_⊙, respectively.

Our derived $\dot{M}$ fall within the range observed by ($\dot{M}$ = ∼1 × 10⁻⁵–1 × 10⁻² M_⊙/yr; Maud et al. 2015). Our $\dot{P}$ and $\dot{E}$ values overlap with the ranges observed by Maud et al. (2015; $\dot{P}$ = ∼7 × 10⁻⁵–1 × 10⁻¹ M_⊙ km s⁻¹ yr, $\dot{E}$ = ∼1 × 10⁻² −1 × 10² L_⊙), but our maximum values are a factor of ∼4 and a factor of ∼20 higher than their observed $\dot{P}$ and $\dot{E}$, respectively.

These trends are likely attributable to our higher angular resolution and our use of SiO rather than CO (both of which allow us to detect emission from smaller regions with higher velocities, i.e., smaller t_dyn), and our use of only one t_dyn for each outflow rather than a unique t_dyn,i for each channel.

Overall, we find that the physical properties we derive for our outflow candidates generally fall within the same ranges as those derived for similar high-mass samples at both protostellar (2000 au) and clump (≥0.1 pc) scales. The deviations we note between our results and those in the literature are likely attributable to differences in angular resolution, molecular tracers used (CO versus SiO), and different methods of deriving dynamical times and the values that depend on them ($\dot{M}$, $\dot{P}$, $\dot{E}$).

We further explore our data at the protocluster level by testing for correlations between our field-aggregated outflow properties and clump-scale properties. In particular, we explore the relationship between total outflow mass, momentum, energy, mass rate, mechanical force, and mechanical luminosity in a given protocluster and clump mass (M_clump), clump bolometric luminosity (L_bol), clump luminosity-to-mass ratio (L_bol/M_clump), and total mass in cores (M_cores,Louvet, M_{cores,isotherm}).

The clump masses are derived using the Point Process Mapping tool (PPMAP; Marsh et al. 2015) to fit a modified blackbody function to far-infrared and millimeter data for each field. We use the ALMA-IMF 1.3 mm continuum mosaics (Ginsburg et al. 2022; Motte et al. 2022), Apex Telescope Large Area Survey of the GALaxy 870 μm images (ATLASGAL; Schuller et al. 2009), and 70–500 μm Photodetector Array Camera and Spectrometer (Poglitsch et al. 2010) and Spectral and Photometric Imaging REceiver (SPIRE; Griffin et al. 2010) data from the Herschel telescope (Pilbratt et al. 2010). For the three fields in W43, we specifically use data from the Herschel imaging survey of OB young stellar objects (Motte et al. 2010; Nguyen-Lu'o'ng et al. 2013); for all other fields, we use data from the Herschel infrared Galactic Plane Survey (Hi-GAL; Molinari et al. 2010, 2016). For W51-IRS2 only, we replace the 250 μm Hi-GAL SPIRE data (which are saturated) with 214 μm data from the Stratospheric Observatory for Infrared Astronomy (SOFIA; Temi et al. 2014, 2018) High-resolution Airborne Wideband Camera Plus instrument (Harper et al. 2018). These data were obtained from SOFIA project 05_0038 (Vaillancourt 2016). We calculate the bolometric luminosities by combining the PPMAP-derived modified blackbody functions with Spitzer images (Werner et al. 2004): 3.6, 4.5, 5.8, and 8.0 μm data from the Infrared Array Camera (Fazio et al. 2004) and 24 μm data from the Multiband Imaging Photometer for Spitzer (Rieke et al. 2004). For all fields, we assume κ_{300 μm} = 0.1 cm² g⁻¹ and β = 1.8, and derive background-subtracted M_clump and L_clump for the full 1.3 mm field of view for each field. Further details of the PPMAP procedure can be found in P. Dell'Ova et al. (2024, in preparation).

Most fields contain only one prominent dust clump, as seen in the ATLASGAL 870 μm data. For the three fields that contain two clumps (G008.67, W43-MM3, W51-IRS2), the value of M_clump we use in this analysis is the sum of the two clumps (i.e., the mass within the entire field of view), not the value of one or the other of the ATLASGAL-detected clumps. Readers are advised that this vocabulary differs from some of the other publications in the ALMA-IMF series.

We use two separate values for total mass in cores in order to account for potential biases in the core mass derivation methods. First, we use the core masses derived in Louvet et al. (2023), who used unique temperature values for each core and take free–free emission and optical depth into account (M_cores,Louvet). This method avoids some of the uncertainties associated with assuming a uniform dust temperature and optical depth for all cores. Second, we derive our own total mass in cores using the flux-density values listed in Appendix D of Louvet et al. (2023) and a uniform temperature of 15 K assuming the optically thin Rayleigh–Jeans approximation (M_{cores,isotherm}). This method gives us a more direct comparison between our results and literature samples, as this is the more common method in the literature. By using both methods, we are able to test whether our results are a result of or are robust against methodological differences. The field-aggregated outflow properties are listed in Table 6. All clump properties except M_{cores,isotherm} can be found in Table 1; M_{cores,isotherm} for each field are listed in Table 6.

For each pair of properties, we calculate both the Kendall τ and Spearman ρ correlation coefficients and their associated p-values. These coefficients and p-values are presented as heatmaps in Figure 5. We find that most relationships are positive, but only a few correlations are significant at the >3σ level (excluding autocorrelations and correlations between dependent properties, e.g., M and P = Mv). The aggregated outflow properties are most strongly correlated with total mass in cores (3σ–5σ; τ and ρ values 0.6–0.9), followed by clump mass and bolometric luminosity (1.5σ–2.5σ, τ and ρ values 0.3–0.6), and finally clump L/M (<1σ, τ and ρ values −0.01–0.2).

Zoom In Zoom Out Reset image size

Using the Spearman ρ test, all outflow properties have >3σ correlation with ${M}_{\mathrm{cores},\mathrm{isotherm}}$ and all but M_out and E_out have >3σ correlation with ${M}_{\mathrm{cores},\mathrm{Louvet}}$. Using the Kendall τ test, ${\dot{M}}_{\mathrm{out}}$ and P_out have >3σ correlation with both variations of M_cores, and ${\dot{P}}_{\mathrm{out}}$ and M_out are correlated at >3σ with ${M}_{\mathrm{cores},\mathrm{isotherm}}$ only. The highest correlation coefficients are between M_out and ${M}_{\mathrm{cores},\mathrm{isotherm}}$ in the Kendall τ test and ${\dot{M}}_{\mathrm{out}}$ and ${M}_{\mathrm{cores},\mathrm{isotherm}}$ in the Spearman ρ test (>4.5σ for both). We discuss selected correlations in greater detail in the following subsections.

Our poor correlation with clump bolometric luminosity can be explained if both (a) the field-aggregated outflow properties are the simple sum of individual outflow+protostar pairings, even though we do not associate any candidates with specific driving sources in this work, and (b) the protostellar population in these fields is not purely coeval. Bontemps et al. (1996) found that outflow mechanical force ($\dot{P}$) decreases with protostellar age for their sample of low-luminosity protostars, with Class I sources following a well-defined linear relationship in log–log space and Class 0 sources lying ∼1 order of magnitude above this line. They suggest that the Class 0 sources follow an evolutionary track down to the best-fit line as they age into Class I sources. Bontemps et al. (1996) also found that the relationship between outflow mass rate ($\dot{M}$) and envelope mass for individual protostars (M_env, or M_core in our data) does not appear to change with time. They suggest that individual protostars will occupy different regions of the $\dot{M}\mbox{--}{M}_{\mathrm{env}}$ plot as they evolve (and thus both their envelope mass and accretion rates decrease), but that the relationship between the two values seems to remain approximately log($\dot{P}$) = −(4.15 ± 0.1) + (1.1 ± 0.15) × log(M_env) for both Class 0 and Class I low-mass sources.

In a population of protostars that is not purely coeval, some protostars will be in the Class 0 stage, some in Class I, and some perhaps in Class II. In the Bontemps et al. (1996) interpretation, the mix of protostellar stages will not have a strong impact on correlations with envelope (core) mass for field-aggregated properties, but it will result in increased scatter in correlations with bolometric luminosity. This is consistent with what we see for our data, and therefore a plausible explanation for this difference in correlation strength; our clump sample size (15) is likely small enough to mask a correlation with luminosity (if present) due to small number statistics. An alternate possibility is that we have additional sources of luminosity (e.g., external irradiation) that are contributing to clump luminosity but are not associated with outflows. Additional analysis in which outflow candidates are associated with specific driving sources, and bolometric luminosities for each driving source are derived, will be needed to answer this question definitively. We therefore do not attempt to fit a relationship between our field-aggregated outflow properties and clump L_bol or clump L/M in this work. Instead, we place these results in context with the literature in Sections 4.2.3 and 4.2.4.

4.2.1. Relationship between M_out, M_cores, and M_clump

In this section, we consider the extent to which outflow mass is driven by the total mass in cores rather than the clump mass. Since most of the material observed in outflows is assumed to be entrained (Bally 2016), both possibilities are worth considering. We find that total outflow mass is correlated with both clump mass and total core mass in a region. We compare these correlations and determine that the mass in cores has a stronger effect.

We use python's scipy.optimize.curve_fit function to determine best-fit linear relations between M_out and M_clump, M_out and M_cores, and M_cores and M_clump in log space. We also calculate ρ and τ correlation coefficients and their p-values for all three cases. The parameters for all best-fit lines, and ρ and τ correlation coefficients and p-values, are shown in Table 7. We discuss only the Louvet et al. (2023) values of M_cores in this section, but the results for the isothermally derived M_cores data are similar and are listed in Table 7.

Table 7. Linear Regression, Kendall's τ, and Spearman ρ Results for log(M_clump), log(M_cores), and log(M_out)

Independent Variable	Dependent Variable	Slope a	Intercept a	Spearman ρ	Kendall's Tau
				ρ, p(ρ)	τ, p(τ)
log(Mclump)	log(Mout)	0.76 ± 0.20	−2.08 ± 0.81	0.62, 0.013	0.47, 0.016
log(Mcores) (L) c	log(Mout)	1.06 ± 0.15	−1.53 ± 0.36	0.80, 3.5 × 10−4	0.67, 5.2 × 10−4
log(Mcores) (I) c	log(Mout)	0.88 ± 0.13	−1.38 ± 0.34	0.88, 1.36 × 10−5	0.77, 1.01 × 10−5
log(Mclump)	log(Mcores) (L)	0.73 ± 0.14	−0.53 ± 0.56	0.79, 4.6 × 10−4	0.59, 0.002
log(Mclump)	log(Mcores) (I)	0.71 ± 0.22	−0.14 ± 0.88	0.63, 0.012	0.47, 0.016
log(Mclump)	log(Mout) - f(Mcores) (L)	−0.01 ± 0.13	−0.05 ± 0.53	0.04, 0.89	0.03, 0.92
log(Mclump)	log(Mout) - f(Mcores) (I)	0.13 ± 0.13	−0.53 ± 0.52	0.31, 0.26	0.22, 0.28

Notes.

^a The slopes, intercepts, and uncertainties for each best-fit line are determined by the scipy.optimize.curve_fit ordinary least-squares fitting package. ^b The ρ and τ correlation coefficients and their p-values are calculated with scipy.stats.spearmanr and scipy.stats.kendalltau, respectively. We convert the p-values to σ-values assuming a normal distribution. The conversions are: 3σ = 2.7 × 10⁻³, 3.5σ = 4.7 × 10⁻⁴, 4σ = 6.3 × 10⁻⁵, and 4.5σ = 6.8 × 10⁻⁶. ^c As in Figure 5, we use two values for M_cores: those whose masses were determined with unique temperatures for each core in Louvet et al. (2023), and those assuming a temperature of 15 K for all cores (see Section 4.1). The M_cores values taken from Louvet et al. (2023) are noted with an "(L)," and those calculated assuming T = 15 k are noted with an "(I)." We discuss only the results derived from M_cores (L) in the text, but present the M_cores (I) results here and in Figures 6 and 7 for completeness.

Download table as:  ASCIITypeset image

Our best-fit line for log(M_out) versus log(M_clump) shows a positive relationship, with a slope of 0.76 and intercept of −2. This fit is shown in the left-hand panel of Figure 6. The p-values for both ρ and τ correspond to a significance of ∼2.5σ. Our fitted slope and intercept are largely consistent with previous literature values: both Li et al. (2018) and Beuther et al. (2002) found slopes of 0.8–0.9 and an intercept of −1 for their fits to log(M_out) versus log(M_clump) for their samples of massive outflows.

Zoom In Zoom Out Reset image size

Our best-fit line for log(M_out) versus log(M_cores) is steeper, with a slope of 1.06 and intercept of −1.5 (Figure 6, middle panel, and Table 7). The correlation coefficients for log(M_out) versus log(M_cores) are both larger and more statistically significant than those we derive for log(M_clump), with ρ and τ p-values corresponding to a significance of ∼3.5σ.

Our best-fit line to log(M_cores) versus log(M_clump) has a similar slope to that of log(M_out) versus log(M_clump), but the ρ and τ p-values correspond to a higher degree of significance: ∼3.25σ for log(M_cores) versus log(M_clump) compared to 2.5σ for log(M_out) versus log(M_clump). The parameters of our best-fit line also do not agree within errors with our fit to log(M_out) versus log(M_cores). This suggests that the relationship between total mass in cores and clump mass is separate from the relationship between outflow mass and total mass in cores. These data and our best-fit line are shown in the right-hand panel of Figure 6.

In order to test whether and how these correlations depend on each other, we subtract the best-fit line to log(M_out) versus log(M_cores) from the log(M_out) data and fit the residuals against log(M_clump). This sequence of simple linear regression and subtraction allows us to determine how much of the correlation between M_out and M_clump can be explained by the relationship between M_out and M_cores—if the residuals still have a noticeable trend and strong ρ and τ coefficients, then M_cores cannot fully explain the correlation between M_out and M_clump.

For each data point, we calculate log(M_out,residual) = log(M_out) − f(M_cores), where f(M_cores) is the y-value predicted by the best-fit line to log(M_out) versus log(M_cores) as shown in the middle panel of Figure 6. Then, we plot each log(M_out,residual) value against its corresponding clump mass, and fit a line with scipy.optimize.curve_fit. Our results are shown in Figure 7.

Zoom In Zoom Out Reset image size

Controlling for log(M_out) versus log(M_cores) significantly decreases the correlation with clump mass. The best-fit line is effectively flat, with a slope and intercept of zero within uncertainties. The τ and ρ correlation coefficients also fall to near zero, and their p-values correspond to <0.5σ significance.

We cannot invert this test—subtracting the best-fit line to log(M_out) versus log(M_clump) and fitting the residuals against log(M_cores)—because M_clump by definition includes M_cores. However, we can safely say that our current results are not consistent with a scenario in which clump mass directly determines total mass in outflows, with no alteration by core mass. The dominant correlation in our analysis is between outflow mass and total mass in cores, not either of the relations involving clump mass. We suggest that the total mass in cores is at least mediating the total mass in outflows to a physically significant degree, and may in fact be dominating it.

4.2.2. Dominance of the Most Massive Outflow in Each Field

The typically large distances to massive star-forming regions mean that the spatial resolution of most outflow surveys in such regions is ≥0.1 pc (Beuther et al. 2002; Maud et al. 2015; Liu et al. 2022). In most cases, this is too large to resolve individual outflows and identify their associated driving sources. Many authors therefore assume, as a first-order approximation, that the total outflow mass is dominated by the most massive individual outflow in the field, presumed to be generated by the most massive protostellar core. For example, Maud et al. (2015) used a simulated, coeval Salpeter population of protostars to test the contribution of massive protostars to the total mechanical force ($\dot{P}$) of the protocluster. They concluded that $\dot{P}$ can be entirely explained by outflows from low- and intermediate-mass protostars up to L = 6400L_⊙, and dominated entirely by the massive protostars above this limit.

In order to test this assumption for our own data, we compare the mass of the most massive outflow (M_out,maximum) with total outflow mass (M_out,total) for each field. These data are plotted in Figure 8. The left-hand panel shows the most massive outflow versus total outflow mass, and the right-hand panel shows the percentage of total mass the most massive outflow accounts for, compared to total outflow mass in the field. There is a correlation between M_out,maximum and M_out,total, as expected.

Zoom In Zoom Out Reset image size

We find that, for our SiO-detected outflows, the most massive outflow in each field is typically responsible for only 12%–30% of the measured total outflow mass, regardless of how much material is contained in outflows overall. The most massive outflow is responsible for the majority of outflow mass (>50%) in only one out of our 15 fields. We also find no significant trend in outflow maximum-mass percentage with total outflow mass. We conclude that the total spatially integrated mass is not dominated by the most massive outflows in our sample.

We also explore this trend in outflow momentum and energy, and mass, momentum, and energy rates. The maximum values for an individual outflow in each field are shown in Table 8, along with their percentage contribution to the field total and the ID number of the candidate responsible. We find similar trends as for outflow mass—there are correlations between P_out,maximum and P_out,total, E_out,maximum and E_out,total, etc., but no indication that the strongest outflow is responsible for >50% of each derived property. We further note that the same outflow is not always responsible for the greatest share of every property. In some fields, the same outflow does dominate all six of M_out, P_out, E_out, ${\dot{M}}_{\mathrm{out}}$, ${\dot{P}}_{\mathrm{out}}$, and ${\dot{E}}_{\mathrm{out}}$, but in other fields, different outflows will be responsible for ${M}_{\mathrm{out},\max }$ and ${P}_{\mathrm{out}\,\max }$, ${P}_{\mathrm{out}\,\max }$ and ${\dot{E}}_{\mathrm{out}\,\max }$, etc. Modifying our test to include the two most massive outflows brings the median M_{out,maximum 2}/M_out,total ratio up to 44%.

Table 8. Strongest Outflows: Absolute and Percentage Contributions, and Candidate ID Numbers ^a

Field	Mmax			Pmax			Emax			Field	${\dot{M}}_{\max }$			${\dot{P}}_{\max }$			${\dot{E}}_{\max }$
	(M⊙)	(%)	(ID)	(M⊙ km s−1)	(%)	(ID)	(1046 erg)	(%)	(ID)		(10−4 M⊙ yr−1)	(%)	(ID)	(M⊙ km s−1 yr−1)	(%)	(ID)	(L⊙)	(%)	(ID)
G008.67	1.11	44	1	23.4	45	1	0.79	49	1	G008.67	0.81	34	1	0.0017	34	1, 3	5.4	39	1
G010.62	2.11	24	17	63	45	17	2.80	66	17	G010.62	1.8	26	17	0.0053	37	17	20	44	12
G012.80	1.18	12	35	45.6	25	35	2.29	39	35	G012.80	3.7	23	35	0.014	35	35	60	47	35
G327.29	0.91	12	18	20.4	16	18	0.57	20	18	G327.29	2.6	13	18	0.0058	17	18	14	21	18
G328.25	0.91	65	2	25.4	80	2	0.93	90	2	G328.25	2.4	71	2	0.007	88	2	20	91	2
G333.60	1.02	14	15	18.2	14	5	0.86	25	5	G333.60	2.7	18	5	0.011	34	5	41	48	5
G337.92	3.7	43	13	44	38	13	0.68	30	13	G337.92	2.1	22	12	0.003	18	5	11	30	5
G338.93	5.7	50	9	116	54	9	3.10	54	9	G338.93	2.7	26	9	0.0058	28	9	16	31	13
G351.77	5.3	47	3 b	135	56	3 b	5.90	71	3 b	G351.77	10	32	3 b	0.026	41	3 b	100	59	3 b
G353.41	0.92	22	11	33	32	11	1.42	41	11	G353.41	1.9	18	11	0.007	26	11	25	30	11
W43−MM1	5.0	11	27 b	141	14	24	6.60	18	24	W43-MM1	9	10	18	0.03	13	18	110	16	18
W43−MM2	2.24	19	17	71	30	17	2.96	39	17	W43-MM2	3.4	22	5	0.011	31	5	38	38	5
W43−MM3	1.13	28	6	20	39	6	0.48	53	6	W43-MM3	0.7	27	10	0.0009	27	10	1.4	33	6
W51−E	10.9	19	20 b	538	35	20 b	34.4	47	20 b	W51-E	56	43	19 b	0.22	48	19 b	1100	52	19 b
W51−IRS2	14.1	30	38 b	382	39	38 b	14.4	43	38 b	W51-IRS2	24	30	38 b	0.06	35	38 b	190	38	38 b

Notes.

^a For each derived physical property (M_max, P_max, E_max, ${\dot{M}}_{\max }$, ${\dot{P}}_{\max }$, ${\dot{E}}_{\max }$), we show the absolute value of the strongest outflow in each field (in M_⊙, M_⊙ km s⁻¹, etc.), the fractional contribution that that outflow makes to the field-aggregated total (in percent, rounded to the nearest integer), and the identification number of that candidate (ID). The same outflow is not always the dominant contributor of each property. Likewise, in one case (G008.67 ${\dot{P}}_{\max }$), two candidates are equally strong; we list both outflow ID numbers in the relevant ID column. ^b These candidates have known contamination from hot-core line emission within the velocity range of the SiO emission. Therefore, the physical properties reported for these candidates (especially E and $\dot{E}$) should be treated as upper limits.

Download table as:  ASCIITypeset image

In Section 3.3.2, we note the possibility of "broken" outflows, i.e., instances in which what we identify as two individual candidates are instead both part of the same larger outflow. This phenomenon is more common with fainter outflow candidates, not the most massive ones, but we nonetheless consider what effect this phenomenon would have on our analysis in this section. If our most massive identified outflow were instead one component of a single, larger outflow encompassing additional SiO emission, this would increase the mass of what we identify as the most massive outflow in each region by up to a factor of 2 (assuming the second component is nearly identical in mass, the largest possible case). If this were the case for every field, this would raise our typical M_out,maximum/M_out,total ratios to 24%–60%. Alternately, this analysis treats the red and blue lobes of bipolar outflows separately; combining the masses of each would have the same effect of, at most, doubling the contribution. Regardless, the typical maximum outflow contribution of 12%−30% (upper limit 65%) is not trivial, but is not large enough for observers to safely neglect contributions from lower-mass outflows.

4.2.3. Outflow Mechanical Force versus Clump Bolometric Luminosity

Bontemps et al. (1996) found that the relationship between outflow mechanical force ($\dot{P}$) and source bolometric luminosity for low-mass protostars evolves with time, with Class I protostars falling along a linear correlation in log–log space (log($\dot{P}$) = −(5.6 ± 0.1) + (0.9 ± 0.15) × log(L_bol)) and Class 0 protostars following an evolutionary track a factor of ∼10 above this line. Duarte-Cabral et al. (2013) found similar results for their sample of nine Class 0 high-mass protostars in Cygnus X, as do van der Marel et al. (2013) for their sample of 16 low-luminosity Class I sources in Ophiuchus. Maud et al. (2015) compared their sample of high-mass protoclusters to the individual protostellar samples of these previous papers, and found that their sample is reasonably well fit by the relationship derived by Bontemps et al. (1996) as well. Maud et al. (2015) additionally derived a $\dot{P}$–L_bol relationship using only their RMS-selected data, and found a slightly shallower best-fit line of log($\dot{P}$ = −4.8+0.61 × log(L_bol).

In Figure 9, we plot outflow mechanical force against clump bolometric luminosity for our sample. The left-hand panel shows our field-aggregated outflow mechanical force for each field as dark blue circles, and the right-hand panel shows the mechanical force of only the most powerful outflow in each field as light blue circles. In order to maintain consistency with previous literature samples, we use our inclination-corrected $\dot{P}$ values for this comparison. Correction factors for $\dot{P}$ can be found in the lower section of Table 5, and we assume a uniform inclination angle of 573 for all candidates (see Bontemps et al. 1996; Maud et al. 2015).

Zoom In Zoom Out Reset image size

We find that our field-aggregated mechanical force values agree well with the Bontemps et al. (1996) best-fit line, and our data extend this line up to L > 10⁶ L_⊙ and $\dot{P}$ ∼ 1.0 M_⊙ km s⁻¹ yr⁻¹ for protocluster-aggregated emission. The agreement between our data and prior literature results is consistent with our earlier interpretation of these field-aggregated values being the sum of individual, well-behaved protostellar outflows with little contamination from ambient emission.

The percentage of total $\dot{P}$ that the most powerful outflow is responsible for ranges from 13%–88%, with a median of 34%. The most powerful outflow is responsible for >50% of the total mechanical force in only one case. In other words, the most powerful outflow in a field is typically responsible for a nontrivial portion of mechanical force but not a majority, consistent with our results in Section 4.2.2 for outflow mass.

We find that considering only the most powerful outflow in each field (${\dot{P}}_{\mathrm{out},\max }$) causes our data to deviate from the best-fit line of Bontemps et al. (1996), and from the larger observational trend established in the literature (Figure 9, right-hand panel). We take this as further consistency with our results in Section 4.2.2 and with the broader literature. This supports the picture of the total mechanical outflow feedback in massive star-forming regions being the sum of multiple individual outflows, and which is poorly described by assuming the aggregate outflow properties are reflective of only the most massive or powerful outflow in the protocluster.

4.2.4. Protocluster Outflow Properties and Clump Evolutionary State

In the course of searching for protostellar outflows, we also identified a significant amount of SiO emission with no high-velocity components and no change in velocity structure with position. These regions are typically elongated in shape, similar to outflows or filaments, and their emission is within ±5 km s⁻¹ of the field V_LSR. Their integrated spectra are often Gaussian or sometimes triangular in shape, with line widths <10 km s⁻¹ in all cases and <6 km s⁻¹ in most. They are found both spatially coincident with and entirely independent from high-velocity SiO emission. An isolated example of this emission is shown in Figure A1.

Zoom In Zoom Out Reset image size

Similar narrow-line, low-velocity SiO emission has been previously reported in, e.g., Codella et al. (1999), Motte et al. (2007), Duarte-Cabral et al. (2014), Csengeri et al. (2016), Louvet et al. (2016), and Minh et al. (2016). The origin of low-velocity SiO emission is still a subject of some debate, but the dominant explanations at present are: (1) the emission has a purely low-velocity origin, e.g., cloud-cloud collisions or slow shocks induced by gravitational collapse, or (2) the emission had a high-velocity origin initially (e.g., in protostellar outflows) and has since cooled kinematically but has not yet frozen out of the gas phase. Duarte-Cabral et al. (2014) additionally suggested that it is possible the SiO abundance was initially enhanced in a high-velocity shock but SiO is being maintained in the gas phase by low-velocity shocks alone. A detailed characterization of the low-velocity SiO emission in our sample, including tests of these possibilities, will be presented in our follow-up paper, A. P. M. Towner et al. (2024, in preparation).

In the field G351.77, we find additional large-scale, high-velocity emission that does not appear to trace individual protostellar outflows. This emission ranges from −94 to +56 km s⁻¹, and spans more than half the field of view. The morphology is bidirectional; blueshifted emission occurs predominantly northeast and north of field center and redshifted emission predominantly southwest and west, but there is little collimation in either direction. The velocity of this gas typically increases with distance from the field center, i.e., appears to exhibit Hubble flow; this trend is especially pronounced in the redshifted emission to the southwest and west. We do not include this high-velocity emission in our catalog. G351.77 is the only field with this exception. A three-color RGB image of G351.77 is shown in Figure A2.

Zoom In Zoom Out Reset image size

We suggest three primary possibilities for this emission, which will be explored in detail in a separate paper. First, this region may contain an "explosive outflow," akin to the explosive event in OMC1 (Bally et al. 2017). In OMC1, the explosion is defined by a high-velocity, roughly spherically symmetric Hubble flow spanning ∼50'' as traced by ¹²CO J = 2 − 1. Here, we find lower velocities than the OMC1 explosion and less spherical symmetry, but there does still appear to be semicoherent gas motion on the cloud scale. The second possibility is that a single massive protostar near the cloud center has recently undergone a significant episodic accretion event (e.g., Caratti o Garatti et al. 2017; Hunter et al. 2017) that ejected material at high enough velocities to produce (at least temporarily) a Hubble flow. The third possibility is that what we interpret as large-scale high-velocity emission does actually originate from individual protostellar outflows whose axes are aligned with each other. We do not favor this latter possibility at present due to the unlikelihood of both outflow axis alignment and redshift/blueshift alignment, but this scenario cannot yet be ruled out completely.

We detect arched or looping structures in the PV diagrams of several of our candidates. These features are telltale signs of bowshock/backsplash in outflows colliding with the ambient medium (Bally 2016); an example from our data set is shown in Figure A3. These features are most common in our 27 bipolar outflow candidates, but do appear in monopolar candidates as well. While we do not examine these structures in detail in this work, the data set presented herein is one of the largest homogeneous interferometric data sets examining outflows in the literature to date. It may therefore be a useful starting point for studies of small-scale outflow physics in the future, particularly intraoutflow structure studies.

Zoom In Zoom Out Reset image size

Both Equations (B1) and (1) assume that the SiO emission is optically thin (τ ≪ 1). This is a common assumption for SiO emission (see, e.g., Lu et al. 2021), especially for higher-energy transitions such as the J = 5 − 4 line. In some cases where direct derivations of optical depth have been done (e.g., Codella et al. 1999, using multiple J-transitions of SiO), SiO lines have also been directly shown to be optically thin. However, as demonstrated by the models of Gusdorf et al. (2008), SiO 5-4 may become optically thick for at least a portion of the shock lifetime in some cases (e.g., 100 yr ≲ t ≲4000 yr for initial shock parameters of 30 km s⁻¹ and n_H = 10⁵ cm⁻³; see Gusdorf et al. 2008). Because we do not have multiple SiO lines available in our data set, we cannot directly solve for τ_SiO for our data. We therefore assume that the lines are optically thin in all cases. This is consistent with our by-eye examination of the SiO data cubes, in which we see no clear evidence of self-absorption for any candidates.

The assumption of a single excitation temperature at all locations is unlikely to be truly physical, and so we investigate the effect of varying T_ex on the resulting SiO column density. We find that column density is not strongly sensitive to excitation temperature in the range 30 K < T_ex < 130 K; in this range, for a given brightness temperature, N_SiO varies by 0.25 dex (a factor of 1.78) for our typical channel width (Δv) of 0.339 km s⁻¹. In Figure B1, we plot this relationship for brightness temperatures between 0.1 and 30 K, in half-dex increments, and find the same results for all T_B tested, as expected. We therefore adopt a flat excitation temperature of 50 K for all outflow candidates, with uncertainties of −20 K and +30 K. This range in temperature translates to changes of −0.03 dex (factor of 1.08, or 8%) and +0.09 dex (factor of 1.22, or 22%), respectively, in the SiO column density.

Zoom In Zoom Out Reset image size

Equation (3) comes from Equation (A5) in Maud et al. (2015), adapted for SiO. The fractional abundance of SiO in Equation (3) presents a particular problem, as it can vary significantly with many factors: the density of the preshock medium, the initial velocity of the shock, and assumptions made about the initial form of elemental Si and the subsequent astrochemical reactions, which create or destroy SiO (Gusdorf et al. 2008). Indeed, Si and SiO have multiple pathways both into and out of the gas phase in the interstellar medium, including direct release of SiO from dust-grain ice mantles or grain cores, release of elemental Si followed by a sequence of astrochemical reactions, and destruction of SiO through the creation of SiO₂ (Schilke et al. 1997). SiO fractional abundance also varies with time, depending on whether a particular parcel of gas is yet to be shocked, experiencing partial or maximum compression from the shock, or undergoing postshock cooling (Schilke et al. 1997; Gusdorf et al. 2008). Because protostellar outflows are typically resolved phenomena, this variation with time also translates to a variation with position within an outflow (Bally 2016).

As we have no independent probe of H₂ column density, we cannot derive SiO fractional abundance directly. Instead, we turn to the shock-chemistry models and previous observational studies as a guide. The models of Schilke et al. (1997) suggest that SiO fractional abundance can range from a few times 10⁻¹¹ to nearly 10⁻⁶, depending on postshock time and the initial location of Si atoms (assuming a shock speed of 30 km s⁻¹). Likewise, the models of Gusdorf et al. (2008) found that SiO fractional abundance can range from as low as a few times 10⁻¹² (if the initial abundance of O₂ ice is assumed to be negligible, and for a shock velocity of 25 km s⁻¹ and a preshock density of 10⁵ cm⁻³) to nearly 10⁻⁶ (assuming an initial abundance of O₂ ice of 1.3 × 10⁻⁵, a preshock density of 10⁴ cm⁻³, and shock velocity ≥35 km s⁻¹). In general, SiO fractional abundance in the ISM has been observed to fall range from 10⁻¹¹–10⁻⁶, with lower values being typical in the ambient medium and values above ∼10⁻⁹ typical in outflows; fractional abundances at or near 10⁻⁶ have been observed toward high-velocity bullets (see Schilke et al. 1997; Codella et al. 1999; Gusdorf et al. 2008, and references therein). For their multispecies analysis of protostellar outflows in massive protoclusters in the central molecular zone, Lu et al. (2021) derived fractional abundances for all species based on H₂ column densities, using HC₃N as an anchor molecule and where the H₂ column densities were derived from dust continuum emission. Their derived SiO abundances range from a few times 10⁻¹⁰–10⁻⁸, with a mean value of 2.05 × 10⁻⁹.

We therefore adopt a flat SiO-to-H₂ abundance ratio of 10^−8.5 (or, 3.16 × 10⁻⁹) for all candidates. This is the midpoint (in log space) of the SiO abundances in both theoretical and observational literature, and within 3% of the mean abundance (in log space) observed by Lu et al. (2021).