© ESO 2022

1 Introduction

Stars form deep inside molecular clouds by accreting their surrounding matter through a circumstellar disk. This process gives rise to strong bipolar jets, which entrain the surrounding envelope material forming extended outflows (e.g., Bally 2016). These protostellar outflows have often been observed, mostly using low-J CO line emission (e.g., Bally & Lada 1983; Bachiller et al. 1990; Bontemps et al. 1996; Maud et al. 2015b; Mottram et al. 2017), which allows for an estimate of their properties, such as outflow velocity and extent. One of the most interesting characteristics of an outflow, which can be inferred from such observations, is the momentum rate or outflow force (F_CO). The outflow force gives the rate at which the outflow injects momentum into the surrounding envelope, and it is often used to describe the strength of an outflow. Correlations have been detected, connecting F_CO to the properties of the source (e.g., luminosity and mass) of the outflow for both protostellar (e.g., Bontemps et al. 1996) and extragalactic (e.g., Lutz et al. 2020) sources. The questions that naturally arise are, whether the physical mechanisms behind these two outflow types are the same (Romero 2021), whether outflows that vary in scale and energetics by many orders of magnitude can indeed be launched by the same physical processes, and how statistical methods can be used to compare and thereby infer similarities and differences between these types of outflows.

The first of the correlations regarding F_CO was initially reported by Bontemps et al. (1996) for a sample of low-mass, primarily Class I protostellar sources, where they found that F_CO scales with the mass of the surrounding envelope (M_env). In the same work, the observed decrease in outflow energetics between Class 0 and Class I sources is attributed to evolutionary effects, due to the dispersion of the surrounding envelope and the reported correlation between F_CO and M_env. The correlation was later expanded to the high-mass protostellar regime by Beuther et al. (2002), using a sample of 26 relatively isolated high-mass sources. The high-mass sample was later significantly strengthened by Maud et al. (2015b), who examined outflows from 99 massive young stellar objects (MYSOs) and very compact H II regions taken from the Red MSX source Survey (RMS, Lumsden et al. 2002, 2013) sample, reporting a similar correlation. Beuther et al. (2002) also reported that the correlation appears tighter between low-mass Class I sources and highmass sources, tentatively suggesting the absence of a proper Class 0 phase in high-mass protostars. The low-mass sample has been expanded in more recent works (Yildiz et al. 2015; Mottram et al. 2017), where it is also found that the commonly used low-J CO observations better describe the time-averaged behavior of an outflow, while other tracers are needed in order to look into current outflow behavior (for example high-J CO transitions, H₂O, or O I). The existence of this correlation points to a connection between M_env and the mass accretion rate (M_acc), which in turn is connected to the outflow (Duarte-Cabral et al. 2013; Mottram et al. 2017).

The second correlation is with the bolometric luminosity of the source (L_bol) (Cabrit & Bertout 1992). In this correlation, a split between Class 0 and Class I sources has been found by Bontemps et al. (1996), with Class 0 sources systematically showing stronger outflows than Class I sources of the same L_bol. Similar to the F_CO - M_env correlation, the F_CO–L_bol correlation was also expanded to the high-mass regime (Beuther et al. 2002; Maud et al. 2015b) and the low-mass sample was expanded by Yildiz et al. (2015); Mottram et al. (2017). In this case though, Maud et al. (2015b) reports a shallower slope for the high-mass sources when compared to the low-mass sources, which, if extrapolated, falls between the Class 0 and Class I sources.

Apart from protostars, galaxies are also often seen driving outflows (e.g., Stuber et al. 2021), albeit on a much larger scale. The predominant mechanisms powering extragalactic outflows are believed to be accretion onto a central supermassive black hole in galaxies hosting active galactic nuclei (AgN) or a combination of stellar processes such as supernovae and star formation in the case of starbursts (see, e.g., Veilleux et al. 2020, for a recent in-depth review on galactic-scale outflows, the energy sources behind them and their driving mechanisms). Extra-galactic outflows are most commonly observed through P-Cygni profiles in OH emission lines (e.g., Sturm et al. 2011; González-Alfonso et al. 2017) or from broad line wings detected in CO line emission (e.g., Cicone et al. 2014; Fluetsch et al. 2019) which allow for a calculation of outflow forces in a way similar to that of the protostellar outflows. Such measurements have been carried out extensively for low-redshift sources (e.g., Cicone et al. 2014; González-Alfonso et al. 2017; Fluetsch et al. 2019; Lutz et al. 2020) and also more recently for higher redshift sources (Spilker et al. 2020). Examining the outflow forces for scaling relations, Lutz et al. (2020) detected a correlation between F_CO and L_bol similar to the one found for protostellar ouflows. A possible connection between protostellar and galactic-scale outflows is further emphasized in Aalto et al. (2020), where the magneto-hydrodynamic (MHD) launch mechanism suggested for the AGN-driven narrow wind and molecular jet of NGC 1377, is expected to be similar to the suggested launch mechanism of protostellar winds.

In this work, we examine 10 sources in the Cygnus-X molecular cloud. Cygnus-X is one of the most intense nearby star-forming regions, hosting a significant number of protoclusters and high-mass protostellar objects deeply embedded within massive, dense, dusty cores (Motte et al. 2007), a large number of ultra-compact H II regions (Cyganowski et al. 2003) as well as a significant population of OB stars (Wright et al. 2010), including Cyg OB2 (Knödlseder 2000), one of the largest known OB associations in our Galaxy. It is estimated that Cygnus-X has a total mass of ≈4 × 10⁶ M_⊙, extends over 100 pc in diameter and the distance to its various substructures ranges from 1.3 kpc for W75/N30 up to 3.3 kpc for AFGL 2591/CygX-S26 (Rygl et al. 2012).

In this paper, we present CO J =3–2 and SiO J = 8–7 emission from the ten sources of the PILS-Cygnus survey (van der Walt et al. 2021), allowing for detection of outflows and measurements of outflow properties. The outline of the paper is as follows: in Sect. 2 we describe the set-up used for the observations, in Sect. 3 we present our results which we discuss further in Sect. 4. Finally in Sect. 5 we present a summary and our conclusions.

2 Observations

The ten sources examined in this work (CygX-N12, N30, N38, N48, N51, N53, N54, N63, S8, S26, where the nomenclature is from Motte et al. 2007, see Table 1) were targeted for the PILS-Cygnus survey (van der Walt et al. 2021), a large frequency-range line survey, covering 32 GHz of continuous frequency observations in the 345 GHz atmospheric window. All the observations were carried out with the SubMillimeter Array (SMA) on Mauna Kea, Hawaii.

The sources were selected due to their location in the same molecular cloud structure and their proximity to the Cyg-OB2 association, which makes them ideal targets to examine the effects of clustered star formation on protostellar outflows. In addition, being classified as intermediate- to high-mass sources they offer a great opportunity to bridge the gap between low- and high-mass sources in protostellar outflow studies, connecting the two samples.

The observations for the PILS-Cygnus survey, from which the data presented in this paper are taken, were performed with the SMA using the SWARM (SMA Wideband Astronomical ROACH2 [second generation Reconfigurable Open Architecture Computing Hardware] Machine) correlator. The observations were carried out with a combination of Compact and Extended configurations with baselines from 7 to 211 meters. All ten sources were observed for approximately the same time over ten tracks, five tracks in the compact configuration and five tracks in the extended.

The compact configuration observations were performed between 27 June and 7 August, 2017, while the extended configuration observations were taken between the 20 of October and the 10 of November 2017. During this time between 6 and 8, antennas where available in the array.

For all the observations, MWC349A was used as the complex gain calibrator, the quasars 3c273, 3c454.3 and 3c84 for bandpass calibration (depending on the exact time of the observations), and Neptune, Titan, Callisto, and Uranus for the flux calibration. We refer to van der Walt et al. (2021) for details on the observations and calibration.

All initial observations had uniform spectral channel widths of 140 kHz, as provided by the SWARM correlator, but were then re-binned into a spectral resolution of 560 kHz (≈0.48 km s⁻¹), in order to improve noise levels. All data calibrations were performed through the Common Astronomy Software Applications (CASA) v4.7 (McMullin et al. 2007).

Subsequent cleaning of the data was also performed through CASA v4.7, using the interactive mode of the clean command. The manual selection of the areas to be cleaned was preferred in this case as to better trace the extended and varying CO line emission. In addition, throughout the cleaning process, natural weighting was used, by applying the Briggs weighting function with a robust parameter of 2.

The resulting clean data have an image size of 256 × 256 pixels and a pixel size of 0″.2. The spectral resolution for the CO data is kept at the initial 0.48 km s⁻¹ but for SiO it is further re-binned to ~ 1.93 km s⁻¹ in order to increase the signal-to-noise ratio (S/N). The rms noise (σ_rms) and the beam sizes of the final cleaned data are shown in Table 2.

3 Results

3.1 Outflow maps

We detected CO 3–2 emission toward all ten sources of our sample. The line profiles show high-velocity wings extending out to >30 km s⁻¹ from the source velocity in the most extreme cases. Furthermore, the profiles are characterized by missing emission at velocities close to that of the source, which is likely a combined effect of optically thick emission as well as spatial filtering by the interferometer. The outflow maps are presented in Figs. 1 and 2 as integrated intensity contours plotted over the corresponding continuum emission.

The velocity ranges used for the integration of the images are given in Table 3 and are relative to the source velocity. For the CO observations, we refer to the velocities closest to the source velocity as “inner” velocities (ν_r,in and ν_b,in for the red and blue-shifted range, respectively), while the highest velocities are called “outer” (v_r,out and v_b,out). For SiO emission, a single range is used in the integration and thus only the outer limits are given, with v_r being the maximum red-shifted velocity and v_b the maximum blue-shifted velocity.

To estimate the velocity ranges used for the CO integrated emission maps a series of different spectra was used. First, the spectrum from the location of peak continuum emission was used for an initial velocity estimate. This velocity estimate was done as to best cover the entirety of the detected line-wings of the CO spectrum. Next, based on the initial velocity estimates, the spectra from the location of peak blue- and red-shifted emission were extracted and used for a new velocity range estimate. This was in some cases necessary as the initial spectra did not always show clear line-wings for both blue- and red-shifted emission. For a final check on the velocity ranges, the spectra from the location of peak integrated blue- and red-shifted emission where used. An example of these spectra, with the adopted velocity ranges marked in red and blue for the red- and blue-shifted part of the outflow respectively is shown in Figs. 3 and 4 for N12. In addition all the spectra used for determining the velocity ranges of the ten Cygnus-X sources can be found in Appendix E.

For the SiO emission, a single velocity range was used covering continuously the entire range from red- to blue-shifted emission. This velocity range was determined so as to maximize the S/N in each case and to give the best possible detection at the 3σ level. As the SiO emission detected was overall weak and not very extended it was only used as an additional tool assisting in determining the origin of the CO outflows.

A clear bipolar outflow structure is visible toward the majority of the sources, and in most cases this appears to be originating from the continuum sources with a few noteworthy exceptions. First, toward N38 only a small outflow is detected that appears to originate from a weak, secondary continuum peak, located slightly north of the primary peak. The presence of relatively strong SiO emission from the same location strengthens the argument that this is indeed the actual source of the outflow. N30 also appears irregular, showing a significant asymmetry between the two lobes. The red-shifted part of the outflow appears stronger and closer to the continuum peak, which we identify as the source of the outflow, while the blue-shifted part of the outflow is small and further away from the continuum peak. This we believe to be a result of missing blue-shifted emission near the source velocity, seen on the spectrum used to determine the velocity integration limits. This pushed the inner blue-shifted velocity significantly further away from the source velocity comparatively with the red-shifted one (For N30 υ_b,in = −24 km s⁻¹ while υ_r,in = 6 km s⁻¹). In the case of N48, an elongated blue-shifted structure is detected slightly to the northeast of the continuum source. This structure is, as mentioned in Duarte-Cabral et al. (2013), the outflow of IRS-1 and therefore not taken into account for this study. Turning to N53, it is clear that, even though there are two distinct continuum peaks, the CO outflow is only detected from the western peak, which might suggest that the second core is not currently driving an outflow. Finally, for the case of S8, we consider the outflow to originate from the area of extended continuum emission under the primary peak, and traveling close to the plane of the sky. In this case, one of the outflow lobes is seen in both red and blue-shifted emission near the image center, while the second lobe is only visible in the red-shifted CO emission extending toward the east of the image. In addition, a fairly strong red-shifted structure appears north of the main continuum peak, but lacks a matching blue-shifted lobe and is therefore not considered an outflow and is not used for the subsequent outflow force calculation.

Fig. 1 Continuum emission of the sources of the PILS-Cygnus survey (Pitts et al. 2022) plotted in gray-scale ranging from −0.1 to 0.5 Jy beam−1 in all maps. Displayed in red contours is the red-shifted CO emission, in blue contours is the blue-shifted CO emission, and in green contours is the integrated SiO emission, when detected. All contour levels are set at 3, 6, 9, 12, 15σrms noise of the corresponding integrated emission. The yellow star marks the assumed origin of the outflow and the light blue and orange lines mark the extend of the outflow lobes. Finally, the beam is marked with the black contour at the bottom left.

Fig. 2 As Fig. 1 but for the remaining sources.

Fig. 3 Spectrum of N12 from the location of the peak of the integrated red-shifted emission. The red marked area corresponds to the velocity range used in the integration of the red-shifted emission.

Fig. 4 Spectrum of N12 from the location of the peak of the integrated blue-shifted emission. The blue marked area corresponds to the velocity range used in the integration of the blue-shifted emission.

3.2 Outflow forces

Properly calculating F_CO is important for constraining the F_CO −M_env and F_CO − L_bol correlations. Recently van der Marel et al. (2013) made a comparison of the 7 most commonly used F_CO calculation methods, where an inferred uncertainty of up to a factor of 6 was found. In the same work, van der Marel et al. (2013) suggests the use of the so-called “separation method” as the one least affected from observational uncertainties and it is the one used here.

In this method, the outflow momentum (P_out) and the dynamical time (t_d) are treated as separate parameters. The dynamical time is calculated by using the maximum outflow velocity (υ_max) and, at the same time, the momentum is calculated from the velocity-weighted integrated intensity. This approach is chosen as it offers a better estimate of the actual outflow momentum since it takes into account the kinematic structure of the outflow (Downes & Cabrit 2007; Curtis et al. 2010). The outflow force using the separation method is calculated as: $$F_{\text{CO}}=c_3\times \frac{K\left({\displaystyle \sum _{\text{j}}{\left[{\displaystyle {\int }_{{\upsilon }_{\text{in}}}^{{\upsilon }_{\text{out,j}}}T\left({\upsilon }^{\prime }\right){\upsilon }^{\prime }\text{d}{\upsilon }^{\prime }}\right]}_{\text{j}}}\right){\upsilon }_{\max }}{R_{\text{lobe}}}.$$(1)

Here c₃ is a factor that corrects for inclination effects shown in Table 4, K accounts for the mass (see Eq. (A.5)) for an assumed excitation temperature of 50 K, and the integral ${\displaystyle {\int }_{{\upsilon }_{\text{in}}}^{{\upsilon }_{\text{out,j}}}T\left({\upsilon }^{\prime }\right){\upsilon }^{\prime }\text{d}{\upsilon }^{\prime }}$ is the velocity-weighted integrated intensity. Finally, R_lobe is the length of the outflow lobe Table 5. The calculation is performed for each outflow lobe individually and then summed up for the total F_CO. For a detailed explanation of the calculation method and the parameters used, see Appendix B. Representative results for each source, calculated for an inclination angle of 30°, are presented in Table 6 and the full results are shown in Table C.1.

Apart from the separation method we used two additional methods to calculate F_CO for all sources in our sample. These were the υ_max method (Method 1) and the 〈υ〉 method (Method 3; van der Marel et al. 2013). The υ_max method was selected because it is the most commonly used method in the literature (van der Marel et al. 2013) therefore presenting F_CO measurements of our sample with the υ_max method might prove useful for further comparisons with similar works. The (υ) method, on the other hand, was selected because it is used by Maud et al. (2015b) for the calculation of the outflow forces of their high-mass protostellar sample. As we aim to compare the results of our work with this sample, having the measurements in both methods will allow us to estimate if any observed difference in F_CO is significant, or is inferred by the calculation method used.

The calculation of the outflow force for both methods is very similar to the separation method so we do not go into detail here. In short, for the υ_max method, the force is calculated as: $$F_{\text{CO}}=c_1\times \frac{K\left({\displaystyle \sum _{\text{j}}{\left[{\displaystyle {\int }_{{\upsilon }_{\text{in}}}^{{\upsilon }_{\text{out,j}}}T\left({\upsilon }^{\prime }\right)\text{d}{\upsilon }^{\prime }}\right]}_{\text{j}}}\right){\upsilon }_{\max }^2}{R_{\text{lobe}}},$$(2)

with c₁ the corresponding inclination angle correction given in Table 4 and the rest of the parameters calculated as in the separation method (see Appendix B). For the 〈υ〉 method, the force is: $$F_{\text{CO}}=c_2\times \frac{K{\left({\displaystyle \sum _{\text{j}}{\left[{\displaystyle {\int }_{{\upsilon }_{\text{in}}}^{{\upsilon }_{\text{out,j}}}T\left({\upsilon }^{\prime }\right){\upsilon }^{\prime }\text{d}{\upsilon }^{\prime }}\right]}_{\text{j}}}\right)}^2}{R_{\text{lobe}}{\displaystyle \sum _{\text{j}}{\left[{\displaystyle {\int }_{{\upsilon }_{\text{in}}}^{{\upsilon }_{\text{out,j}}}T\left({\upsilon }^{\prime }\right)\text{d}{\upsilon }^{\prime }}\right]}_{\text{j}}}},$$(3)

with c₂ again being the corresponding inclination correction factor, given in Table 4. The rest of the calculation is performed as for the separation method, with the only difference that υ_max is not needed in this case. The outflow forces, as measured with the υ_max and 〈υ〉 methods, are presented in Appendix C (Tables C.2 and C.3, respectively).

In Fig. 5, the ratio of outflow forces as measured using the different calculation methods for the 10 sources of the PILS-Cygnus survey are plotted as a function of their envelope masses. This comparison reveals small differences between the different methods, with the results varying by a factor of <3. This result falls well within the factor of 6 that van der Marel et al. (2013) reports among all calculation methods discussed therein.

Focusing on the comparison between the separation method and the 〈υ〉 method reveals a difference between the methods that is up to a factor of ~2. The very small difference found between the two methods is important as it means that a direct comparison between results from Maud et al. (2015b) and other works, using the separation method, is possible and will not be significantly affected by the choice of outflow force calculation method.

Regarding the different inclination angles, for all three methods, we find that the difference, even between the most extreme cases of going from i = 10° to i = 70°, is smaller than an order of magnitude. Therefore, even though constraining the precise inclination angle of an outflow is difficult, this uncertainty does not affect the measured outflow force significantly.

Fig. 5 Comparison of different outflow force calculation methods. Top: ratio of the outflow forces for the 10 sources of the PILS-Cygnus survey, as measured using the υmax over the separation method. Bottom: ratio of the outflow forces for the 10 sources of the PILS-Cygnus survey, as measured using the 〈υ〉 over the separation method. In both cases, red points mark the results for i = 10°, blue for i = 30°, green for i = 50° and orange for i = 70°. Also plotted in black and gray dotted lines are the 1:1 and 2:1 ratios respectively.

3.2.1 Uncertainties on the outflow forces

It is important to highlight that the outflow force measurements for the PILS-Cygnus sources likely are lower limits of the true F_CO. This is a result of the lower velocity limits used for the integration in Eq. (1) (the inner velocities in Table 3). Due to the optical thickness of the CO emission in combination with spatial filtering effects arising from the nature of interferometric observations, a velocity range of ≈10 up to ≈30 km s⁻¹ surrounding υ_source is excluded from the calculation. Therefore, the low-velocity outflowing gas, that holds a significant part of the outflow mass (Dunham et al. 2014a), is not taken into account in measuring F_CO. The outflow force scales with the square of the outflow velocity υ², but depends only linearly on the outflow mass. As a result, the higher-velocity material contributes significantly more to the total outflow force compared to the slow moving material. Thus, not taking into account the smaller velocities of the outflow in the calculation of the outflow force likely has a smaller impact on the resulting outflow forces. More precisely, Offner et al. (2011) found that integrating over velocities >2 km s⁻¹ leads to an underestimate of the outflow mass by a factor of 5–10 from simulations, a result that was confirmed observationally later by Dunham et al. (2014a). Interestingly, N30, S26, and N63, the sources driving the strongest outflows, are also the sources with the largest masked-out velocity range. This is most likely due to, based on the F_CO − M_env correlation, these sources having larger envelopes and therefore being susceptible to spatial filtering effects at higher velocities, which in turn forces the inner velocities limits to higher velocities (in absolute for the blue-shifted ones).

We do not expect the difference in F_CO to be significant, but even if our measurements underestimate the true F_CO by as much as an order of magnitude, the results for the PILS-Cygnus sources would still fall within the scatter of the high-mass sample from Maud et al. (2015b). Therefore, although the outflow measurements presented in Table C.1 are lower limits of the true F_CO, our results regarding clustered star formation remain unchanged.

3.2.2 Sample biases

An important issue that needs to be addressed when studying samples such as here, are the selection criteria used for the observations of the high- and low-mass protostellar samples (Figs. 6 and Fig. 7). The sample of Yildiz et al. (2015) was specifically chosen to have strong and bright outflows (the original sample was from a Herschel key program; see van Dishoeck et al. 2011, for more details). Similarly, the sample from Mottram et al. (2017) was based on a luminosity cut-off to target the brightest protostars in nearby clouds, which presumable coincides with the sources driving the strongest outflows. As a result, these samples are most likely representative of the upper limit of the actual low-mass distribution. In contrast, the sample presented in Maud et al. (2015b) has no such biases. Including a more representative sample of low-mass sources might allow for a better connection of the high- and low-mass samples, which currently appear separate. In addition, such a sample would of course yield more accurate statistics for the correlations followed by the protostellar samples.

4 Discussion

To better be put into context, the resulting outflow forces are compared with the correlations of the outflow force with the envelope mass and bolometric luminosity already established for protostellar sources (e.g., Bontemps et al. 1996; Beuther et al. 2002). Specifically, this comparison is first performed for the F_CO − envelope mass correlation for protostars and then it is expanded to include extragalactic sources when comparing to the luminosity.

4.1 M_env determination

As demonstrated by Bontemps et al. (1996), outflow forces correlate with the mass of the protostellar envelope. The envelope masses used for the Cygnus-X sources, shown in Table 1, are mostly taken from Pitts et al. (2022) who used SED fitting of Herschel observations assuming a distance of 1.4 kpc, the average distance to the Cygnus-X sources (Rygl et al. 2012), with the exception of S26 that was initially calculated from Pitts et al. (2022) for the distance of 3.3 kpc (Rygl et al. 2012), also assumed in this work. For sources N38 and N54, no envelope masses are given in Pitts et al. (2022), therefore previous measurements from Motte et al. (2007) are used. These are scaled to the assumed distance of 1.4 kpc.

Pitts et al. (2022) used SCUBA 450 µm observations to determine the size of the envelopes in order to measure their mass. These observations have a resolution of 14 arcesonds, significantly larger than the <2 arcseconds of the SMA observations used to map the CO outflows. These SMA observations reveal multiple substructures that are included in the same envelope, but are not all driving outflows. In such cases, the entire envelope mass is attributed to the single observed outflow source even though it is likely to break down into multiple sources driving separate outflows. A clear example of this is N53, where the continuum emission (Fig. 1) shows two separate cores, but the CO contours reveal that only one of them actually drives an outflow and is therefore assigned the entire envelope mass. Similar cases are the outflows of S8, N48 and N38. Especially for N38, the origin of its outflow is located at a secondary continuum peak north of the primary continuum peak. This structure remains unresolved in the observations used by Motte et al. (2007) to estimate the envelope mass of the source. It is therefore reasonable to assume that the envelope mass found in Motte et al. (2007) is mostly representative of the primary continuum peak, and not of the actual outflow source. In an attempt to use a more representative estimate of the actual envelope mass of the source behind the detected outflow we scaled down the mass found by Motte et al. (2007) by a factor of 4, and used this new value (see Table 1) for this work. The scaling factor of 4 was determined by the ratio of the intensity of continuum emission of the primary and secondary peak. We recognize of course that this estimate is very uncertain, and that the actual envelope mass of N38 as a result is highly uncertain.

4.2 Correlation with envelope mass

To explore the F_CO–M_env correlation we plot in Fig. 6 the outflow forces versus the envelope mass for the ten sources in the Cygnus-X molecular cloud. These represent the extent of the resulting outflow forces for the four different inclination correction factors. Alongside our results, the low-mass protostellar sample is plotted. Finally, the high-mass sample from Maud et al. (2015b) is plotted. All four of these studies use CO J =3–2 observations, similar to this study. We also note that Yildiz et al. (2015), Mottram et al. (2017) and van der Marel et al. (2013) use the separation method for calculating F_CO but Maud et al. (2015b) uses the 〈υ〉 method. As previously discussed, the 〈υ〉 method, on average, varies from the separation method by up to a factor of 2. This uncertainty is significantly smaller than the inherent uncertainties in the F_CO measurements therefore a direct comparison of the different measurements is reasonable.

Finally, the blue dashed line represents the best-fit on the correlation as calculated by van der Marel et al. (2013) for a sample of low-mass protostars.

The PILS-Cygnus survey sources are, in their majority, tightly packed in the same molecular cloud structure and located in close proximity to the Cyg-OB2 association, one of the largest OB associations in our Galaxy; they therefore are a prime example of clustered star formation. On the other hand, Maud et al. (2015b) specifically targeted more isolated high-mass sources. The good agreement between the results of the two studies for sources of similar envelope masses, as can be seen in Fig. 6, suggests that clustered star formation does not significantly impact the outflow activity of protostars. An exception to this may be N38, which appears to have a slightly weaker outflow than sources of similar mass. As discussed in Sect. 3.1, N38 has only a compact outflow, that originates from a secondary continuum peak close to the main peak which causes great uncertainty for its true envelope mass. Even so, given the large scatter of the high-mass sample, N38 largely still follows the same correlation as the rest of the sample.

It is believed that the outflow force is directly connected to the mass accretion rate of the driving source (e.g., Bontemps et al. 1996; Duarte-Cabral et al. 2013). Therefore, the previous result, that clustered star formation does not significantly affect protostellar outflow activity, suggests that the accretion process itself is not affected by the surrounding environment of the protostar but is rather dominated by local processes. Of course, the sample of clustered sources examined is very limited, containing only 10 sources all taken from the same star-forming region. In order to verify this result further studies of clustered protostellar sources are required.

Additionally, in Fig. 6, a clear distinction can be seen between the low- and high-mass protostellar sources, highlighted by the lack of measurements for sources with M_env ≈ 10 M_⊙. The few intermediate-mass sources plotted from van Kempen et al. (2009) seem more in line with the low-mass sample but the sample consists of too few sources to allow for any concrete conclusions. Properly connecting the two samples will assist significantly in determining if they actually follow the same correlation or not.

Fig. 6 Outflow forces for the ten sources of the PILS Cygnus survey over their envelope mass. The outflow forces are plotted using colored rectangles covering the range of possible results for the different inclination correction factors. Plotted in the same figure are, Class 0 sources from Yildiz et al. (2015) in red triangles, Class I sources from Yildiz et al. (2015) in blue triangles, Class 0 sources from Mottram et al. (2017) in red circles, Class I sources from Mottram et al. (2017) in blue circles, Class I sources from van der Marel et al. (2013) in blue stars, intermediate-mass sources from van Kempen et al. (2009) in dark gray diamonds and high mass sources from Maud et al. (2015b) in light gray squares. The blue dashed line shows the best fit as found in van der Marel et al. (2013) for CO J = 3–2 observations of a low mass protostellar sample.

4.3 Correlation with luminosity

To explore the F_CO–L_bol correlation we used a sample of pro-tostellar and extragalactic sources collected from the literature in which we added the outflow forces measured in this work for the Cygnus-X sources. The protostellar sample consists of low-mass protostellar sources from van der Marel et al. (2013), Yildiz et al. (2015) and Mottram et al. (2017) and high-mass protostars from Maud et al. (2015b) and this work, while the extragalactic sample is taken from Lutz et al. (2020) with additional sources taken from Fluetsch et al. (2019). This sample includes both AGN galaxies, split into Seyfert I, II, and LINERs, and starburst (SB) galaxies. Both Lutz et al. (2020) and Fluetsch et al. (2019) used the time-averaged thin-shell approach for calculating the outflow force. In this approach, a constant outflow rate over the entire life-time of the outflow (t_flow = R_out/υ_out) is assumed, yielding a final density profile of ρ ∝ r⁻². The outflow rate is then given as: $${\dot{M}}_{\text{out}}=\frac{M_{\text{out}}{\upsilon }_{\text{out}}}{R_{\text{out}}},$$(4)

and the outflow force is calculated as: $$F_{\text{CO}}={\dot{M}}_{\text{out}}{\upsilon }_{\text{out}}.$$(5)

This approach is very similar to the υ_max method (van der Marel et al. 2013) often used to calculate protostellar outflows where a unique velocity is assumed for the entirety of the outflow. Figure 7 shows the outflow forces over the corresponding bolometric luminosities for all sources. The L_bol values used for the Cygnus-X sample, presented in Table 1, are mostly from Pitts et al. (2022), with the exception of N38 and N54 which are taken from San José-García (2015) and Kryukova et al. (2014) respectively.

In Fig. 7, the same characteristics displayed in Fig. 6 (i.e., good agreement between the Cygnus-X and the Maud et al. (2015b) sources and a clear separation between the low- and high-mass sample) can be seen. We note though that the separation between the low- and high-mass samples appears more pronounced in the F_CO–L_bol correlation (Fig. 7) than in the F_CO–M_env correlation (Fig. 6). Therefore, overall this comparison reinforces the points made in Sect. 4.1.

The most important questions arising from Fig. 7 is how well the F_CO–L_bol correlations for the different object types compare, and if it is at all possible to assume a common correlation for all three (low-mass, high-mass and extragalactic) samples. In order to address these questions, we calculated the best-fit of the F_CO–L_bol correlation for each of the samples. In addition to these, we also examined the AGN- and SB-driven outflows separately, as well as the Class 0 and Class I sources. The least-squares method is used to find the best-fits, which are of the form: $$\log \left(F_{\text{CO}}\right)=a\log \left(L_{\text{bol}}\right)+b.$$(6)

The a and b parameters describing the fits for the various sub-samples are shown in Table 7 and the correlations themselves are shown in Fig. D.1.

To properly compare the various fits a statistical approach is needed. As a first step, we performed an F-test on the variances of the samples, comparing them in pairs. Briefly, to perform an F-test we need to calculate the F-value, which is given as: $$F=\frac{s_1}{s_2},$$(7)

where s₁ and s₂ are the variances to be compared, with s₁ > s₂. The F-value is then compared with a critical value (F_crit), which depends on the degrees of freedom (DF) of both samples. These are calculated as N – 1 for each sample, with N being the number of sources in the sample. If the F-value calculated is smaller than F_crit, then the two variances are consistent with being similar and the two samples can be compared. In this work we use the F-table for α = 0.05, which corresponds to a 95% confidence level.

For the pairs with comparable variances we then performed two-tailed t-tests in order to examine if the parameters of their best-fit, the slope (a) and the intercept (b), are similar or not. For a t-test the initial hypothesis is that the parameters are equal (in this case the intercept, H₀: b₁ = b₂) and we calculate the t-value as: $$t=\frac{b_1-b_2}{\sqrt{s_{b_1}+s_{b_2}}},$$(8)

where b₁ and b₂ are the compared parameters of the two fits, and $s_{b_1}$ and $s_{b_2}$ the corresponding variances. The t-value is then used alongside the degrees of freedom, which is in this case is N₁ + N₂ − 4 with N₁ and N₂ the sizes of the two samples whose fits we are comparing, in order to calculate the P-value.

If the resulting P-value is lower than the confidence level α then the initial hypothesis can be rejected and the compared parameters are not equal. On the other hand, if the P-value exceeds α then the initial hypothesis cannot be rejected and we assume that the two parameters are equal. The confidence level is again α = 0.05 resulting in a 95% level of confidence in our results.

We present in Table 8 the F-values calculated for all comparisons among the seven samples as well as the degrees of freedom corresponding to each sample. Comparing these values with the appropriate critical values from the F-table we were able to determine which samples are comparable with each other and are marked in boldface in Table 8. Some of the most interesting results though are that the variances for the three extragalactic samples (the entire sample, AGN, and SB) are comparable with one another, that the high-mass protostellar sample variance is only comparable with the SB sample, and also that the Class 0 sample is comparable with all three extragalactic samples.

We performed t-tests for the best-fits of the sample pairs found having comparable variances and present the resulting P-values along with the corresponding DF in Table 9. The P_slope value corresponds to the comparison of the slopes of the fits, while P_inter to the comparison of the intercepts. All pairs marked in boldface in Table 9 can statistically be considered to follow the same correlation.

At this point, it is also important to note that two methods are often used for calculating outflow forces from extragalactic sources. In the first, the outflow force is calculated assuming a time-averaged thin shell approximation (e.g., Lutz et al. 2020; Fluetsch et al. 2019) (described earlier in Sect. 4.3). In the second method, the outflowing gas is assumed to fill the entire cavity of the outflow with a constant density (e.g., Cicone et al. 2014; Fiore et al. 2017). This alternate calculation method is found to differ by a factor of three compared to the thin shell approach (Lutz et al. 2020). Both extragalactic samples used in this work use the thin shell method, therefore in order to test if the selected method of outflow-force calculation has any impact on our results, we increased the outflow forces of the extragalactic sample by a factor of three, and performed again the statistical analysis for the new best-fits. The resulting differences in F-values where minimal (< 10⁻⁴) but the P-values concerning the comparison of the estimated intersections varied by up to 30% in some cases. Even the largest of these changes did not have a statistically significant effect on any of our results. Therefore, at this point, there is no apparent preference toward either of the methods used for calculating outflow forces, and the following conclusions appear to stand for both scenarios.

The first interesting result of this statistical comparison is that, while Maud et al. (2015b) argue that outflow forces from low-mass Class I sources and high-mass sources follow the same correlation with bolometric luminosity, based on an extrapolation of the two correlations, here it is found that the two samples, Class I and high-mass outflows, show significantly different variances, and thus they are likely drawn from two different populations. Going beyond the Class I outflows, the variance is also significantly different between the Class 0 and total low-mass outflows when compared to the high-mass outflows. Specifically, the low-mass outflows show larger variances compared to their high-mass counterparts. This discrepancy is unlikely to come from the use of different calculation methods for the outflow force, as this only adds an offset in a log-log plot. Instead, it is possible that a significant mass fraction of the high-mass outflows is CO-dark, which means that it is not traced by CO lines. This could be in the form of C⁺ (Leurini et al. 2015; Gusdorf et al. 2015), where the ionization is driven by the higher UV-luminosity of high-mass sources. An alternative is that the low-mass outflows inherently show more variation, for example due to episodic accretion (e.g., Audard et al. 2014). However, exploring these scenarios in detail is beyond the scope of this work.

The lack of agreement between the low- and high- mass samples found here is in agreement with the apparent separation of the two samples detected in both Figs. 6 and 7 and discussed in Sect. 4.1.

The most interesting result arising from this comparison though is that all the different extragalactic sources appear to follow the same correlation as low-mass Class 0 protostellar sources. A result that hints toward a connection between the two otherwise very different samples.

Additionally, an intriguing question that comes from the apparent agreement between high-mass protostellar outflows and extragalactic outflows from starburst galaxies, found in the previous statistical comparison, is whether this is due to the dominant role of high-mass protostellar outflows in the total extragalactic outflow activity. To address this, an order-of-magnitude calculation was done, assuming an initial mass function (IMF) of (Kroupa 2002): $$\xi \left(M_{*}\right)=\{\begin{array}{ll}{\left(M_{*}/0.08M_{\odot }\right)}^{-0.3},\hfill & 0.01\le M_{*}/M_{\odot }<0.08\hfill \\ {\left(M_{*}/0.08M_{\odot }\right)}^{-1.3},\hfill & 0.08\le M_{*}/M_{\odot }<0.5\hfill \\ 0.09{\left(M_{*}/0.5M_{\odot }\right)}^{-2.3},\hfill & 0.5\le M_{*}/M_{\odot }<1.0\hfill \\ 0.02{\left(M_{*}/0.5M_{\odot }\right)}^{-2.7},\hfill & 1.0\le M_{*}/M_{\odot }\hfill \end{array},$$(9)

where ξ is the number of stars for a given stellar mass M_*. The IMF is defined for main-sequence stellar masses whereas the masses shown in Fig. 6 are envelope masses. In order to translate stellar masses to envelope masses a core-efficiency of ~30% is assumed (e.g., André et al. 2010), that is, approximately 30% of the envelope or core mass ends up in the star. While this assumption comes with a number of caveats (see, e.g., Offner et al. 2014, for a full discussion), it is the best available at this moment. To get the contribution to the overall outflow force from protostellar sources with M_env we weighted the outflow force with the number of stars at a given envelope mass, where the outflow force is estimated from the best-fit presented in van der Marel et al. (2013): $$\log \left(F_{\text{CO}}\right)=-4.4+0.86\log \left(M_{\text{env}}\right).$$(10)

We selected this fit due to the good agreement it shows with the high-mass sample, as can be seen in Fig. 6. To get the total outflow contribution for the two samples, the result is integrated over 0.03 ≤ M_env/M_⊙ < 30 for the low-mass sample and over 30 ≤ M_env/M_⊙ < 3 × 10⁴ for the high-mass sample. The resulting outflow forces show that low-mass protostellar sources are responsible for ~95% of the total protostellar outflow force, while high-mass sources are responsible for the remaining ~5%.

The above result is based on a comparison between outflow forces and envelope masses, but a similar correlation exist between the outflow force and the luminosity. The question then is whether this correlation will lead to the same conclusion. Pitts et al. (2022) compared envelope masses and bolometric luminosities for protostars with luminosities in the range of <1 L_⊙ to >10⁵ L_⊙, and found that they correlate with envelope mass as $L_{\text{bol}}\propto M_{\text{env}}^{1.2}$. Thus, since the power-law index is shallower than that of the IMF at M > 0.5 M_⊙, we would also conclude that the low-mass protostars drive the outflow contribution.

The correlations are based on observations of a number of high-mass star-forming regions, which likely contain both high-and low-mass protostars. In these cases, the regions contain a single dominant outflow that is driven by the high-mass protostar, a result that is common to other studies (e.g., Beuther et al. 2002). This would suggests that the outflows from high-mass protostars dominate the output from clusters, in contradiction with the result above. However, the time over which high-mass sources accrete and drive outflows is significantly shorter than that of their low-mass counterparts (~10⁵ vs. 5 × 10⁵ years; e.g., Motte et al. 2018; Dunham et al. 2014b, and references therein). Furthermore, it is unclear when low-mass stars form in clusters compared to the high-mass stars, more precisely whether they form at the same time or if the low-mass stars form before or after the high-mass stars. In this analysis, we are looking at the integrated effect of the outflows from all stars formed in a cluster, implying that at certain points in time the outflow from a single high-mass protostar may very well dominate the outflow output.

Even so, this result, although an order-of magnitude estimate, suggests that it might be low-mass protostars that dominate the total extragalactic outflow activity especially in the case of the starburst galaxies, where the outflows are believed to originate from forming stars. Therefore, the detected agreement of starburst galaxies and high-mass protostellar sources cannot be easily attributed to the dominance of high-mass protostellar outflows when averaging over a large number of clusters in the total outflow activity of starburst galaxies, unless the initial mass function is significantly more top-heavy in these galaxies.

Fig. 7 Outflow forces over bolometric luminosities for a sample of low-mass protostars, a sample of high-mass protostars, and extragalactic sources. Red points mark Class 0 low-mass protostars, with triangles taken from Yildiz et al. (2015) and circles from Mottram et al. (2017). Blue points mark Class I low-mass protostars, again with triangles from Yildiz et al. (2015), circles from Mottram et al. (2017), and stars from van der Marel et al. (2013). Gray squares mark the high-mass protostellar sample from Maud et al. (2015b) while gray crosses mark the intermediate mass sources described in van Kempen et al. (2009). The Cygnus-X sources, measured in this work, are shown with orange rectangles covering the extent of their outflow forces for the different available inclination correction factors. For the extragalactic sample, AGN sources are marked in magenta, with left-looking triangles for Seyfert I, right-looking triangles for Seyfert II and pentagons for LINER sources. Purple diamonds mark the starburst sources. Finally, filled points represent measurements taken from Lutz et al. (2020) and empty ones from Fluetsch et al. (2019).

5 Summary and conclusions

This paper presents interferometric observations of CO J =3–2 and SiO J = 8–7 emission lines of ten high-mass protostellar sources in the Cygnus-X molecular cloud, taken from the PILS-Cygnus survey and carried out with the SubMillimeter Array. In addition, an extended sample of low-mass protostellar, high-mass protostellar and extragalactic outflows with measured outflow forces, was gathered from the literature. Using these data, the main conclusions arising from this work, are as follows:

• Using integrated intensity maps of the CO J = 3–2 line emission, molecular outflows were detected originating from all ten sources of the PILS-Cygnus survey, and their respective outflow forces were successfully measured using both the separation, the 〈υ〉, and the υ_max methods;

• The close agreement found between the F_CO, measured for the ten Cygnus-X sources, with results from the literature for sources with similar M_env and L_bol, indicate that clustered star formation has no significant effect on the outflow behavior of protostars. As the outflow activity is tied to the accretion process, this result suggests that protostellar accretion is not significantly affected by the surrounding environment, but it is rather tied to local processes, determined by the protostar itself;

• Examination of the F_CO–M_env and F_CO–L_bol correlations for the protostellar sources, reveals clear differences between the low- and high-mass protostellar sources both in the correlation itself as well as in the observed scatter;

• The statistical comparison of the best-fits to the F_CO – L_bol correlation, for the different types of sources in the literature sample, revealed, with 95% confidence, that extragalactic and low-mass Class 0 protostellar outflows follow the same correlation, thus suggesting the presence of a common underlying launch mechanism;

• An order-of-magnitude calculation, comparing the total outflow force contribution of low- and high-mass protostars, suggests that low-mass protostars dominate the total proto-stellar contribution in extragalactic starburst outflows.

The results presented in this paper allow for some interesting conclusions, but they also highlight the need for additional research and observations on the matter. More specifically, additional and more accurate observations are certainly required in order to constrain the observed protostellar correlations. Constraining these correlations and increasing the number of sources from regions of clustered star formation would solidify the conclusions presented in this work. This need is even more apparent in the case of intermediate-mass sources, required to properly connect the low- and high-mass protostellar samples, as such observations would allow for a better understanding of the observed differences between the various subsamples and the mechanisms behind it thus offering a more clear view of the overall behavior of outflows on all scales.

At this point though the data clearly point toward a correlation between the physical parameters of Galactic and extra-galactic outflows strongly supporting the existence of a common launching mechanism for both types of outflows. As such, this paves the way for developing a universal outflow launch model.

Acknowledgments

Acknowledgements. We would like to thank the anonymous referee for suggestions that greatly helped in clarifying the presentation of the results, and for constructive comments on the manuscript in general. We would also like to thank Susanne Aalto for discussing extragalactic outflows, and Sarel van der Walt for assistance with the SubMillimeter Array data. The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica. We acknowledge and thank the staff of the SMA for their assistance and continued support. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have had the opportunity to conduct observations from this mountain. The research of IMS and LEK is supported by a research grant (19127) from VILLUM FONDEN.

Appendix A Outflow mass derivation

In order to measure the mass of the outflow, it is standard to assume that the emission from the outflow is optically thin (Dunham et al. 2014a) and that the CO level populations in the outflowing gas are in local thermodynamic equilibrium (LTE) (see Appendix A in Maud et al. 2015a)). In this case, the column density of the upper level of the observed transition is: $$N_{\text{u}}=\frac{8\pi k_{\text{B}}}{h{c}^3}\frac{v^2\int T_{\text{B}}\text{d}\upsilon }{A_{\text{ul}}},$$(A.1)

where k_B is Boltzmann’s constant, h is Planck’s constant, c the speed of light, v the frequency of the emission line, and T_B the brightness temperature, which is connected to the observed intensity as I_v = 2kT_Bv²/c², and is integrated over the velocity range covered by the outflow. The total column density is then given by: $$N_{\text{CO}}=\frac{Q\left(T_{\text{ex}}\right)}{g_{\text{u}}}{\text{e}}^{\left(E_{\text{u}}/k_{\text{B}}{T}_{\text{ex}}\right)}{N}_{\text{u}},$$(A.2)

with Q(T_ex) the partition function, g_u the degeneracy of the upper level, E_u the energy of the upper level, and T_ex the excitation temperature of the observed gas. Multiplying with the abundance ratio of CO to H₂, [H₂/CO], gives the total column density for H₂ and then the total mass can be calculated as: $$M=\mu m_{\text{H}}{N}_{{\text{H}}_2}A,$$(A.3)

where μ is the mean molecular weight, m_H the hydrogen mass, and A is the observed area. By combining Eq. A.1, A.2, A.3 we get that the mass is: $$M=\mu m_{\text{H}}A\left[\frac{{\text{H}}_2}{\text{CO}}\right]\frac{Q\left(T_{\text{ex}}\right)}{g_{\text{u}}}{\text{e}}^{\left(E_{\text{u}}/k_{\text{B}}{T}_{\text{ex}}\right)}\frac{8\pi k_{\text{B}}{v}^2}{h{c}^3{A}_{\text{u1}}}{\displaystyle \int T_{\text{B}}\text{d}\upsilon }.$$(A.4)

All the constants can be gathered in a conversion factor K used to convert integrated intensity into mass. K is then defined as: $$K=\mu m_{\text{H}}A\left[\frac{{\text{H}}_2}{\text{CO}}\right]\frac{Q\left(T_{\text{ex}}\right)}{g_{\text{u}}}{\text{e}}^{\left(E_{\text{u}}/k_{\text{B}}{T}_{\text{ex}}\right)}\frac{8\pi k_{\text{B}}{v}^2}{h{c}^3{A}_{\text{u1}}}.$$(A.5)

In this work, for calculating K, a mean molecular weight μ = 2.8 was used in order to also take into account He (Kauffmann et al. 2008), a CO to H₂ abundance ratio of [H₂/CO] = 1.2 × 10⁴ (Frerking et al. 1982) and an excitation temperature T_ex = 50 K, as used in van der Marel et al. (2013); Yıldız et al. (2012). The values used for the partition function is Q = 18.5 (corresponding to T_ex of 50 K), the upper level degeneracy g_u = 7, the upper level energy E_u = 33.2 K, and the Einstein A-coefficient A_ul = 10^−5.6 s⁻¹; these values are taken from Splatalogue¹ and correspond to the CO J = 3−2 transition examined in this work. The selected temperature is found to not significantly affect the results as, for example, using T_ex = 100 K only leads to an increase of the final outflow forces by a factor of ~1.4 Lastly, A is now the physical area covered by one pixel and is calculated assuming a distance of 1.4 kpc, except for S26 for which a distance of 3.3 kpc is used. The distances used are consistent with those used in Pitts et al. (2022) and are based on measurements from Rygl et al. (2012).

Appendix B The separation method

In the separation method, the outflow force is calculated by: $$F_{\text{CO}}=c_3\times \frac{K\left({\displaystyle \sum _{\text{j}}{\left[{\displaystyle {\int }_{{\upsilon }_{\text{in}}}^{{\upsilon }_{\text{out,j}}}T\left({\upsilon }^{\prime }\right){\upsilon }^{\prime }\text{d}{\upsilon }^{\prime }}\right]}_{\text{j}}}\right){\upsilon }_{\max }}{R_{\text{lobe}}}.$$(B.1)

In this equation, c₃ is a correction factor for the inclination of the outflow and can take the values shown in Table 4 according to the inclination angle of the outflow. Due to the large uncertainties in determining the actual inclination angle of the outflows from their morphology c3 takes four different values representing the inclination angle ranges of 0° − 30°, 30° − 50°, 50° − 70°, and 70° − 90° respectively. For this work, the correction factors are calculated, following the method presented in van der Marel et al. (2013), by taking the mean of the $t_{{\upsilon }_{\max }}$ / True age values for mass density contrasts of 0.1 and 1 presented in Table 6 in Downes & Cabrit (2007) for velocities of 7 − 40 km s⁻¹. Contrary to van der Marel et al. (2013), who used the table values for a velocity range of 2 − 40 km s⁻¹, the 7 − 40 km s⁻¹ velocity range is deemed more appropriate here, where most of the emission at low velocities is missing. In this work, the outflow forces are calculated using all four values of c₃ as to both make up for the uncertainty inherent in determining the inclination angle of an outflow but also to examine the effect these uncertainties actually have on the final result.

Next, the sum of the integrated intensities is performed over all pixels contributing to the outflow to get the total velocity-weighted integrated intensity of the outflow. In order to properly determine the pixels contributing to the outflow, we isolate the outflow area, aiming to avoid including unrelated structures in the measurements, and then sum up all pixels with an integrated intensity higher than the 3σ_rms (σ_red and σ_blue for the red and blue-shifted lobe respectively) level used in the contour maps of the outflows (Fig. 1) while the velocity-weighted data cube is created by multiplying each channel with the corresponding velocity.

The integration itself is performed for each pixel separately, with velocity limits v_in and v_out,j. The lower velocity limit v_in is the same for all pixels but v_out,j is calculated individually for each pixel as the highest velocity at which the intensity remains above the 1 σ_rms level. For the lower limit v_in, the inner velocities v_b,in and v_r,in shown in table 3 are used, while the upper limit v_out,j is calculated as described previously using as σ_rms limit the single channel rms noise of each source, presented in Table 2.

Then, v_max is the maximum outflow velocity, used for the calculation of t_dyn, and is calculated as the maximum (in absolute value for the blue-shifted lobes) of the v_out,j for each outflow lobe. Finally, R_lobe, the projected extent of the outflow lobes is measured directly from the contour maps. The values of v_max and R_lobe measured for our sample are given in table 5.

Appendix C Outflow forces for different inclination angles and calculation methods

Outflow forces as calculated using the Separation, v_max and <v> methods for all four different possible inclination angles.

Appendix D Outflow force - bolometric luminosity correlation best fits

F_CO - M_env correlation plots, with the calculated best-fits of the various subsamples plotted with dashed lines.

Fig. D.1 Outflow forces over bolometric luminosities for a sample of low-mass protostars, a sample of high-mass protostars and extragalactic sources. In each plot, the sample, from which the plotted best-fit of the FCO - Lbol was taken, is highlighted. (1) Best-fit to all extragalactic sources. (2) Best-fit to AGN extragalactic sources. (3) Best-fit to starburst extragalactic sources. (4) Best-fit to high-mass protostellar sources. (5) Best-fit to all low-mass protostellar sources. (6) Best-fit to Class I low-mass protostellar sources. (7) Best-fit to Class 0 low-mass protostellar sources.

Appendix E Spectra for determining integration velocity ranges.

All the spectra used for determining the velocity ranges used to create the CO integrated emission maps of the Cygnus-X sources presented in Section 3.1.

E.1 N12

The spectra of N12 used to determine the velocity ranges for the creation of the integrated CO emission maps.

Fig. E.1 Spectrum of N12 from the location of peak continuum emission.

Fig. E.2 [Left]: Spectrum of N12 from the location of peak red-shifted emission. [Right]: Spectrum of N12 from the location of peak blue-shifted emission.

Fig. E.3 [Left]: Spectrum of N12 from the location of the peak of the integrated red-shifted emission. [Right]: Spectrum of N12 from the location of the peak of the integrated blue-shifted emission.

E.2 N30

The spectra of N30 used to determine the velocity ranges for the creation of the integrated CO emission maps.

Fig. E.4 Spectrum of N30 from the location of peak continuum emission.

Fig. E.5 [Left]: Spectrum of N30 from the location of peak red-shifted emission. [Right]: Spectrum of N30 from the location of peak blue-shifted emission.

Fig. E.6 [Left]: Spectrum of N30 from the location of the peak of the integrated red-shifted emission. [Right]: Spectrum of N30 from the location of the peak of the integrated blue-shifted emission.

E.3 N38

The spectra of N38 used to determine the velocity ranges for the creation of the integrated CO emission maps.

Fig. E.7 Spectrum of N30 from the location of peak continuum emission.

Fig. E.8 [Left]: Spectrum of N38 from the location of peak red-shifted emission. [Right]: Spectrum of N38 from the location of peak blue-shifted emission.

Fig. E.9 [Left]: Spectrum of N38 from the location of the peak of the integrated red-shifted emission. [Right]: Spectrum of N38 from the location of the peak of the integrated blue-shifted emission.

E.4 N48

The spectra of N48 used to determine the velocity ranges for the creation of the integrated CO emission maps.

Fig. E.10 Spectrum of N48 from the location of peak continuum emission.

Fig. E.11 [Left]: Spectrum of N48 from the location of peak red-shifted emission. [Right]: Spectrum of N48 from the location of peak blue-shifted emission.

Fig. E.12 [Left]: Spectrum of N48 from the location of the peak of the integrated red-shifted emission. [Right]: Spectrum of N48 from the location of the peak of the integrated blue-shifted emission.

E.5 N51

The spectra of N51 used to determine the velocity ranges for the creation of the integrated CO emission maps.

Fig. E.13 Spectrum of N51 from the location of peak continuum emission.

Fig. E.14 [Left]: Spectrum of N51 from the location of peak red-shifted emission. [Right]: Spectrum of N51 from the location of peak blue-shifted emission.

Fig. E.15 [Left]: Spectrum of N51 from the location of the peak of the integrated red-shifted emission. [Right]:Spectrum of N51 from the location of the peak of the integrated blue-shifted emission.

E.6 N53

The spectra of N53 used to determine the velocity ranges for the creation of the integrated CO emission maps.

Fig. E.16 Spectrum of N53 from the location of peak continuum emission.

Fig. E.17 [Left]: Spectrum of N53 from the location of peak red-shifted emission. [Right]: Spectrum of N53 from the location of peak blue-shifted emission.

Fig. E.18 [Left]: Spectrum of N53 from the location of the peak of the integrated red-shifted emission. [Right]: Spectrum of N53 from the location of the peak of the integrated blue-shifted emission.

E.7 N54

The spectra of N54 used to determine the velocity ranges for the creation of the integrated CO emission maps.

Fig. E.19 Spectrum of N54 from the location of peak continuum emission.

Fig. E.20 [Left]: Spectrum of N54 from the location of peak red-shifted emission. [Right]: Spectrum of N54 from the location of peak blue-shifted emission.

Fig. E.21 [Left]: Spectrum of N54 from the location of the peak of the integrated red-shifted emission. [Right]: Spectrum of N54 from the location of the peak of the integrated blue-shifted emission.

E.8 N63

The spectra of N63 used to determine the velocity ranges for the creation of the integrated CO emission maps.

Fig. E.22 Spectrum of N63 from the location of peak continuum emission.

Fig. E.23 [Left]: Spectrum of N63 from the location of peak red-shifted emission. [Right]: Spectrum of N63 from the location of peak blue-shifted emission.

Fig. E.24 [Left]: Spectrum of N63 from the location of the peak of the integrated red-shifted emission. [Right]: Spectrum of N63 from the location of the peak of the integrated blue-shifted emission.

E.9 S8

The spectra of S8 used to determine the velocity ranges for the creation of the integrated CO emission maps.

Fig. E.25 Spectrum of S8 from the location of peak continuum emission.

Fig. E.26 [Left]: Spectrum of S8 from the location of peak red-shifted emission. [Right]: Spectrum of S8 from the location of peak blue-shifted emission.

Fig. E.27 [Left]: Spectrum of S8 from the location of the peak of the integrated red-shifted emission. [Right]:Spectrum of S8 from the location of the peak of the integrated blue-shifted emission.

E.10 S26

The spectra of S26 used to determine the velocity ranges for the creation of the integrated CO emission maps.

Fig. E.28 Spectrum of S26 from the location of peak continuum emission.

Fig. E.29 [Left]: Spectrum of S26 from the location of peak red-shifted emission. [Right]: Spectrum of S26 from the location of peak blue-shifted emission.

Fig. E.30 [Left]: Spectrum of S26 from the location of the peak of the integrated red-shifted emission. [Right]: Spectrum of S26 from the location of the peak of the integrated blue-shifted emission.

References

All Tables

All Figures

Fig. 1 Continuum emission of the sources of the PILS-Cygnus survey (Pitts et al. 2022) plotted in gray-scale ranging from −0.1 to 0.5 Jy beam−1 in all maps. Displayed in red contours is the red-shifted CO emission, in blue contours is the blue-shifted CO emission, and in green contours is the integrated SiO emission, when detected. All contour levels are set at 3, 6, 9, 12, 15σrms noise of the corresponding integrated emission. The yellow star marks the assumed origin of the outflow and the light blue and orange lines mark the extend of the outflow lobes. Finally, the beam is marked with the black contour at the bottom left.

In the text

	Fig. 2 As Fig. 1 but for the remaining sources.
In the text

	Fig. 3 Spectrum of N12 from the location of the peak of the integrated red-shifted emission. The red marked area corresponds to the velocity range used in the integration of the red-shifted emission.
In the text

	Fig. 4 Spectrum of N12 from the location of the peak of the integrated blue-shifted emission. The blue marked area corresponds to the velocity range used in the integration of the blue-shifted emission.
In the text

Fig. 5 Comparison of different outflow force calculation methods. Top: ratio of the outflow forces for the 10 sources of the PILS-Cygnus survey, as measured using the υmax over the separation method. Bottom: ratio of the outflow forces for the 10 sources of the PILS-Cygnus survey, as measured using the 〈υ〉 over the separation method. In both cases, red points mark the results for i = 10°, blue for i = 30°, green for i = 50° and orange for i = 70°. Also plotted in black and gray dotted lines are the 1:1 and 2:1 ratios respectively.

In the text

Fig. 6 Outflow forces for the ten sources of the PILS Cygnus survey over their envelope mass. The outflow forces are plotted using colored rectangles covering the range of possible results for the different inclination correction factors. Plotted in the same figure are, Class 0 sources from Yildiz et al. (2015) in red triangles, Class I sources from Yildiz et al. (2015) in blue triangles, Class 0 sources from Mottram et al. (2017) in red circles, Class I sources from Mottram et al. (2017) in blue circles, Class I sources from van der Marel et al. (2013) in blue stars, intermediate-mass sources from van Kempen et al. (2009) in dark gray diamonds and high mass sources from Maud et al. (2015b) in light gray squares. The blue dashed line shows the best fit as found in van der Marel et al. (2013) for CO J = 3–2 observations of a low mass protostellar sample.

In the text

Fig. 7 Outflow forces over bolometric luminosities for a sample of low-mass protostars, a sample of high-mass protostars, and extragalactic sources. Red points mark Class 0 low-mass protostars, with triangles taken from Yildiz et al. (2015) and circles from Mottram et al. (2017). Blue points mark Class I low-mass protostars, again with triangles from Yildiz et al. (2015), circles from Mottram et al. (2017), and stars from van der Marel et al. (2013). Gray squares mark the high-mass protostellar sample from Maud et al. (2015b) while gray crosses mark the intermediate mass sources described in van Kempen et al. (2009). The Cygnus-X sources, measured in this work, are shown with orange rectangles covering the extent of their outflow forces for the different available inclination correction factors. For the extragalactic sample, AGN sources are marked in magenta, with left-looking triangles for Seyfert I, right-looking triangles for Seyfert II and pentagons for LINER sources. Purple diamonds mark the starburst sources. Finally, filled points represent measurements taken from Lutz et al. (2020) and empty ones from Fluetsch et al. (2019).

In the text

Fig. D.1 Outflow forces over bolometric luminosities for a sample of low-mass protostars, a sample of high-mass protostars and extragalactic sources. In each plot, the sample, from which the plotted best-fit of the FCO - Lbol was taken, is highlighted. (1) Best-fit to all extragalactic sources. (2) Best-fit to AGN extragalactic sources. (3) Best-fit to starburst extragalactic sources. (4) Best-fit to high-mass protostellar sources. (5) Best-fit to all low-mass protostellar sources. (6) Best-fit to Class I low-mass protostellar sources. (7) Best-fit to Class 0 low-mass protostellar sources.

In the text

	Fig. E.1 Spectrum of N12 from the location of peak continuum emission.
In the text

	Fig. E.2 [Left]: Spectrum of N12 from the location of peak red-shifted emission. [Right]: Spectrum of N12 from the location of peak blue-shifted emission.
In the text

	Fig. E.3 [Left]: Spectrum of N12 from the location of the peak of the integrated red-shifted emission. [Right]: Spectrum of N12 from the location of the peak of the integrated blue-shifted emission.
In the text

	Fig. E.4 Spectrum of N30 from the location of peak continuum emission.
In the text

	Fig. E.5 [Left]: Spectrum of N30 from the location of peak red-shifted emission. [Right]: Spectrum of N30 from the location of peak blue-shifted emission.
In the text

	Fig. E.6 [Left]: Spectrum of N30 from the location of the peak of the integrated red-shifted emission. [Right]: Spectrum of N30 from the location of the peak of the integrated blue-shifted emission.
In the text

	Fig. E.7 Spectrum of N30 from the location of peak continuum emission.
In the text

	Fig. E.8 [Left]: Spectrum of N38 from the location of peak red-shifted emission. [Right]: Spectrum of N38 from the location of peak blue-shifted emission.
In the text

	Fig. E.9 [Left]: Spectrum of N38 from the location of the peak of the integrated red-shifted emission. [Right]: Spectrum of N38 from the location of the peak of the integrated blue-shifted emission.
In the text

	Fig. E.10 Spectrum of N48 from the location of peak continuum emission.
In the text

	Fig. E.11 [Left]: Spectrum of N48 from the location of peak red-shifted emission. [Right]: Spectrum of N48 from the location of peak blue-shifted emission.
In the text

	Fig. E.12 [Left]: Spectrum of N48 from the location of the peak of the integrated red-shifted emission. [Right]: Spectrum of N48 from the location of the peak of the integrated blue-shifted emission.
In the text

	Fig. E.13 Spectrum of N51 from the location of peak continuum emission.
In the text

	Fig. E.14 [Left]: Spectrum of N51 from the location of peak red-shifted emission. [Right]: Spectrum of N51 from the location of peak blue-shifted emission.
In the text

	Fig. E.15 [Left]: Spectrum of N51 from the location of the peak of the integrated red-shifted emission. [Right]:Spectrum of N51 from the location of the peak of the integrated blue-shifted emission.
In the text

	Fig. E.16 Spectrum of N53 from the location of peak continuum emission.
In the text

	Fig. E.17 [Left]: Spectrum of N53 from the location of peak red-shifted emission. [Right]: Spectrum of N53 from the location of peak blue-shifted emission.
In the text

	Fig. E.18 [Left]: Spectrum of N53 from the location of the peak of the integrated red-shifted emission. [Right]: Spectrum of N53 from the location of the peak of the integrated blue-shifted emission.
In the text

	Fig. E.19 Spectrum of N54 from the location of peak continuum emission.
In the text

	Fig. E.20 [Left]: Spectrum of N54 from the location of peak red-shifted emission. [Right]: Spectrum of N54 from the location of peak blue-shifted emission.
In the text

	Fig. E.21 [Left]: Spectrum of N54 from the location of the peak of the integrated red-shifted emission. [Right]: Spectrum of N54 from the location of the peak of the integrated blue-shifted emission.
In the text

	Fig. E.22 Spectrum of N63 from the location of peak continuum emission.
In the text

	Fig. E.23 [Left]: Spectrum of N63 from the location of peak red-shifted emission. [Right]: Spectrum of N63 from the location of peak blue-shifted emission.
In the text

	Fig. E.24 [Left]: Spectrum of N63 from the location of the peak of the integrated red-shifted emission. [Right]: Spectrum of N63 from the location of the peak of the integrated blue-shifted emission.
In the text

	Fig. E.25 Spectrum of S8 from the location of peak continuum emission.
In the text

	Fig. E.26 [Left]: Spectrum of S8 from the location of peak red-shifted emission. [Right]: Spectrum of S8 from the location of peak blue-shifted emission.
In the text

	Fig. E.27 [Left]: Spectrum of S8 from the location of the peak of the integrated red-shifted emission. [Right]:Spectrum of S8 from the location of the peak of the integrated blue-shifted emission.
In the text

	Fig. E.28 Spectrum of S26 from the location of peak continuum emission.
In the text

	Fig. E.29 [Left]: Spectrum of S26 from the location of peak red-shifted emission. [Right]: Spectrum of S26 from the location of peak blue-shifted emission.
In the text

	Fig. E.30 [Left]: Spectrum of S26 from the location of the peak of the integrated red-shifted emission. [Right]: Spectrum of S26 from the location of the peak of the integrated blue-shifted emission.
In the text