© The Authors 2023

1 Introduction

High-mass stars have a significant impact on their environments and on galaxy evolution globally through their ionising radiation, stellar winds, and their deaths in supernova explosions (Zinnecker & Yorke 2007). Already during the earliest stages of their formation, massive protostars might inject significant amounts of energy and momentum into the interstellar medium (ISM) in the form of outflows, capable of disrupting clumps and cores (Beuther et al. 2002; Bally 2016).

Outflows, a ubiquitous phenomenon in both low- and highmass star-forming regions, play an essential role in transporting angular momentum and regulating the star-forming process across multiple spatial scales (Bally & Lada 1983; Evans 1999). Both the dissipation of the envelope material and mass loss via the outflows lower the core-to-star formation efficiency (Krumholz et al. 2014; Offner & Chaban 2017). At cluster and clump scales, outflows drive turbulence that provides additional support against gravitational collapse (Frank et al. 2014).

Outflows are typically detected using low-J (J ≲ 5) velocity-resolved rotational lines of carbon monoxide (CO), which is the second most abundant molecule in the ISM (CO/H₂ = 1.2 × 10⁻⁴, Frerking et al. 1982). The low-lying rotational levels of CO are easily collisionally excited even at low densities and can be observed readily at millimetre wavelengths. These lines constitute a useful diagnostic of the gas kinetic temperature of outflows (Bally & Lada 1983; Yıldız et al. 2015). An extensive search for outflows traced by such low-J CO lines towards a total of 2052 massive star-forming clumps that were identified in the APEX Telescope Large Area Survey of the Galaxy (ATLASGAL, Schuller et al. 2009), provided an overall outflow detection rate of 58% (Yang et al. 2018a, 2022).

Observations of high-J CO (J ≳ 10) lines provide an opportunity to study denser and warmer parts of star-forming clumps and the outflows that arise in them. Recent surveys with the Herschel Space Observatory (Pilbratt et al. 2010) found that CO lines account for the bulk far-infrared (IR) gas cooling in both low- and high-mass star-forming regions (Karska et al. 2013, 2014, 2018; van Dishoeck et al. 2021). The velocity-resolved profiles of high-J CO towards low-mass protostars revealed a significant contribution of high-velocity (υ ~ 20–30 km s⁻¹) gas to the total far-IR line emission (San José-García et al. 2013; Yıldız et al. 2013), and a similarity to the H₂O emission likely arising from the same gas (San José-García et al. 2016; Kristensen et al. 2017). Single-pointing observations towards high-mass sources have also revealed broad, outflow wings in high-J CO line profiles, but they have been limited to just a few sources: W3 IRS5 (San José-García et al. 2013), AFGL 2591 (Kaźmierczak-Barthel et al. 2014), Orion KL, Orion S, Sgr B*, and W49N (Indriolo et al. 2017). Complementary observations have been obtained with the German REceiver for Astronomy at Terahertz frequencies¹ (GREAT, Risacher et al. 2018) on board the Stratospheric Observatory For Infrared Astronomy (SOFIA, Young et al. 2012). High-resolution spectroscopy of far-IR CO lines from an intermediate-mass protostar Cep E revealed an extremely high-velocity (EHV) gas (υ up to ~140 kms⁻¹), tracing shocks associated with the jet and intermediate-to-high velocity gas (υ from 50 to 100km s⁻¹) associated with outflow cavities and a bow shock (Gómez-Ruiz et al. 2012; Lefloch et al. 2015; Gusdorf et al. 2017). The line profiles of CO 16−15 towards two highmass sources, however, lacked the EHV component and revealed broad line wings extending up to υ ~ 50km s⁻¹ (Leurini et al. 2015; Gusdorf et al. 2016). Other surveys, conducted with the PACS and SPIRE instruments aboard the Herschel Space Telescope (Pilbratt et al. 2010), lacked the high spectral resolution necessary to disentangle the envelope and outflow emission in the spectra (Karska et al. 2014; Goicoechea et al. 2013, 2015).

For this study, we used SOFIA/GREAT to quantify high-J CO emission towards 13 high-mass star-forming clumps with the aim to isolate the contribution from the outflows and estimate excitation conditions associated with the line wing emission. We also examined how the high-J CO emission varies as a function of clump properties and evolutionary stages.

The paper is organised as follows: Sect. 2 describes the source sample, observations with SOFIA, and the complementary CO observations with the APEX telescope and Herschel.

In Sect. 3, we present line profiles of high-J CO transitions (Sect. 3.1) and decompose the emission that belongs to the line wings (Sect. 3.2). In addition, we study the correlations of velocity-integrated emission with source properties, and those of the fraction of wing emission with source evolutionary stages (Sect. 3.3). Subsequently, we analyse the excitation of high-J CO lines using LTE and non-LTE approaches (Sects. 3.4 and 3.5). Section 4 consists of the discussion of our results in the context of previous studies and Sect. 5 presents a summary and our conclusions.

2 Observations

2.1 Sample

All sources have been selected from the ATLASGAL survey covering 420 deg² of the inner Galactic plane in the 870 µm dust continuum (Urquhart et al. 2014; König et al. 2017). The latest version of the ATLASGAL source catalogue contains 5007 clumps spanning a wide range of masses (M_clump) and luminosities (L_bol), and it is divided into four evolutionary stages – quiescent, protostellar, young stellar objects (YSOs), and H II regions (H II). For more details, readers can refer to Urquhart et al. (2022).

For this work, we originally selected a representative sample of 20 sources grouped within four star-forming regions in the Galactic plane. Among them, 13 sources within three regions were successfully observed with SOFIA. Table 1 shows the final list of sources with the overview of their properties and evolutionary stages. The sample consists of three protostellar (24d), seven YSOs (IRb), and three H II regions (HII), with L_bol from 1.6 × 10³ to 4.6 × 10⁵ L_⊙ and M_clump from 1.6 × 10² to 2.3 × 10³ M_⊙ (König et al. 2017; Urquhart et al. 2019, 2022).

2.2 SOFIA observations and data reduction

Observations of the CO 11−10 and 16−15 lines were collected using the SOFIA/GREAT (Heyminck et al. 2012; Risacher et al. 2016) and upGREAT receivers (Risacher et al. 2018). Our programme ‘Probing high-J CO through the evolution of high-mass star-forming clumps’ (project IDs 02_0102 and 03_0103; PI: F. Wyrowski) ran during Cycle 2 (2014 May) and Cycle 3 (2016 May).

GREAT was a high-resolution, dual-colour spectrometer (R ≥ 10⁷) initially designed for single-beam observations. In 2014, we used its L1 and L2 channels to obtain simultaneous coverage of bands in the 1.25–1.52 THz and 1.80–1.90 THz windows, respectively. In 2016, we combined GREAT's L1 channel with the upGREAT Low Frequency Array (LFA), which covered the 1.83–2.07 THz window in two polarisations. The 7-pixel hexagonal setup of the LFA provided spatial information about the line emission, whereas each pixel had a full-width at half-maximum (FWHM) beam size of 14.8″ on the sky². The corresponding beam size in the L1 channel was 19.9″ in 2014 and 19.1″ in 2016. The higher frequency L2 channel provided a Field of View (FoV) of 14.1″. The adopted main beam efficiencies (η_MB) are 0.7 (in 2014) and 0.66 (in 2016) for the L1 channel, 0.65 for L2 channel, and 0.65 for the central spaxel of LFA. Data were processed and reduced by SOFIA/GREAT staff and released at product level 3 where first order baselines have been subtracted. Most of the data are ready to use, except for CO 16−15 spectra of G13.66−0.6 where an additional third order baseline was subtracted. Spectral resolutions are presented in Table 2. To perform the analyses without any spectral resolution bias, all spectra were smoothed to a common resolution of 1.0km s⁻¹.

For G12.81−0.2 and G351.25+0.7, additional line observations were collected using the SOFIA/4GREAT receiver (Duran et al. 2021). The observations were carried out in 2019 under project ‘high-J CO observations towards high-mass star-forming clumps’ (project ID 83_0711; PI: H. T. Dat). 4GREAT was a single-beam system with four sub-receivers (4G-1 to 4G-4) and could observe four spectral windows simultaneously. The 4G-3 and 4G-4 modules, which cover the 1.24−1.52 THz and 2.49− 2.69 THz windows, were tuned to map the CO 13−12 and 22−21 transitions. The maps were scanned in 5 × 5 grids with centres 6″ away from each other. The typical beam sizes for the 4G-3 and 4G-4 modules are 20″ and 10.5″, respectively (Duran et al. 2021). Main beam efficiencies are 0.7 for 4G-3 and 0.57 for 4G-4. Observations of the CO 22−21 line are affected by instrumental standing waves that make it difficult to confidently detect line emission. The noise levels of averaged spectra range from 0.40 K to 0.88 K at Δυ of 0.6km s⁻¹.

Data reduction for the CO 13−12 line was performed with the CLASS programme, which is part of the GILDAS³ software developed by the Institut de Radioastronomie Millimétrique (IRAM). A second order baseline was subtracted, and the spectra were also smoothed to an adequate 1.0km s⁻¹. For this study, we extracted averaged spectra within a beam of 20″.

2.3 Additional observations and ancillary data

Additional single pointing observations of ¹³CO10-9 and C¹⁸O 9−8 were conducted with the Herschel-Heterodyne Instrument for the Far-Infrared (HIFI, de Graauw et al. 2010) on board of the Herschel Space Telescope. Observations for ten sources (Appendix D) were obtained as part of project ‘A Water survey of massive star-forming clumps in the inner Galaxy’ (project ID OT2_fwyrowsk_3, PI: F. Wyrowski). In addition, archival data for G351.44+0.7 using Herschel/HIFI were taken from the ‘Water in star-forming regions with Herschel’ programme (San José-García et al. 2013; van Dishoeck et al. 2021). Data from the H and V polarisations of the wide-band spectrometer were averaged. Baselines lower than third order were also subtracted. The spectra were converted to a T_MB scale using a forward efficiency of 0.96 and a main beam efficiency of 0.64 for the ¹³CO 10−9 line and 0.74 for the C¹⁸O9−8 line, respectively. Finally, the spectra were smoothed to 1.0 km s⁻¹. Angular and original spectral resolutions are listed in Table 2.

We also used high spectral resolution ¹³CO6−5 and C¹⁸O6−5 data from Dat et al. (in prep.) and CO 6−5 from Navarete et al. (2019). All three transitions were observed with the CHAMP⁺ receiver (Kasemann et al. 2006; Güsten et al. 2008) at the Atacama Pathfinder Experiment (APEX) 12 m sub-millimetre telescope (Güsten et al. 2006). The on-the-fly (OTF) scans resulted in datacubes of 80″ × 80″ with angular resolutions of ~9″. For comparisons with the higher-J CO observations, averaged spectra with an effective beam size of 20″ around the sources were extracted and then smoothed to 1.0 km s⁻¹ for all three lines.

3 Results and analysis

3.1 Line detections

In this section, we examine detection rates of high-J CO lines towards high-mass clumps from ATLASGAL. We also present their line profiles and quantify the correlations of integrated intensities with the sources’ properties.

Figure 1 shows the spectra of high-J CO lines towards the central position of high-mass clumps from our sample (see also Table 1). The pattern of emission is generally compact, based on additional observations offset from the clump centres towards four sources (see Appendix A).

The CO 11−10 line is detected at 3σ or higher levels towards all sources, which span a broad range of evolutionary stages and have diverse properties. The CO 16−15 line, however, is firmly detected towards ten out of 13 clumps; G13.66−0.6 and G34.41+0.2 show only 2σ peaks and G14.19−0.2 shows a non-detection. In addition, the CO 13−12 line was successfully observed and detected towards G12.81−0.2 and G351.25+0.7. In Appendix B, the peak and integrated intensity of the detected lines are given.

The line profiles of clump central positions exhibit a broad line wing emission, suggesting the presence of outflows (Fig. 1). The median full width at zero power⁴ (FWZP) of 45km s⁻¹ was measured for the CO 11−10 line and 33 km s⁻¹ for the CO 16−15 line (Appendix C). The broadest profile, with FWZP of 165 km s⁻¹, is seen towards G351.77−0.5 where high-velocity gas has been detected in CO 2−1 and 6−5 lines (Leurini et al. 2009). However, multiple pointing observations show a lack of a EHV gas component towards the central source; it is only detected at offset outflow positions (Leurini et al. 2009), consistent with the analysis of the outflow emission from an intermediate mass protostar Cep E (Gómez-Ruiz et al. 2012; Lefloch et al. 2015; Gusdorf et al. 2017). The lack of clear evidence of EHV gas towards our sources may also result from the beam dilution, and could only be addressed using high angular resolution observations (e.g. Cheng et al. 2019).

The velocity ranges of high-J CO lines resemble those detected in CO 6−5 towards the same sources (Fig. 2). Self-absorption features are seen in the CO 11−10 line profiles towards G351.25+0.7 and G351.77−0.5. In addition, G12.81−0.2 and G35.20−0.7 have tilted peaks, which could be an indication of self-absorption. The latter source also shows a sign of self-absorption in the CO 16−15 line. Other profile asymmetries, in particular the triangular blue-wing shape of G351.16+0.7, resemble those of high-J CO emission from a photodissociation region in M17 SW (Pérez-Beaupuits et al. 2015).

For G34.26+0.15, an additional narrow peak is seen at ~38km s⁻¹ in both the CO 11−10 and CO 16−15 spectra. This feature is an artefact due to over-corrected mesospheric CO, which shows the limitations of the adopted atmospheric model (see also, Gusdorf et al. 2016).

For G34.40-0.2, the line profiles of high-J CO lines seem to be shifted by ~1km s⁻¹ from the source velocity obtained from the C¹⁸O9−8 line (Fig. 1). The uncertainty of the Gaussian fit to the C¹⁸O line is smaller than 0.25 km s⁻¹, and thus cannot account for the observed shift, suggesting that it may be caused by self-absorption. Small velocity shifts are also present in the line profiles of other objects, for example G34.26+0.15 and G351.25+0.7.

We calculated CO line luminosities, L_CO, as , where D is the distance to the source (Table 1) and is the velocity integrated flux in Wm⁻¹. The flux conversion from Kkm s⁻¹ to Wm⁻¹ followed Eq. (1) in Indriolo et al. (2017). Figure 3 shows the correlations between L_CO and source properties (Table 1). The significance of the correlations was quantified by the Pearson correlation coefficient r, which also depends on the number of data points N (Marseille et al. 2010).

Both CO 11−10 and CO 16−15 line luminosities show weak correlations (r of 0.63−0.66) with the clump mass, M_clump, primarily tracing a cold gas and dust reservoir (König et al. 2017). Stronger correlations (r of 0.85–0.95) are found for the high-J CO line luminosities and clump bolometric luminosities, L_bol, in line with previous studies using CO 10−9 (see Sect. 4). It is noteworthy that clumps at different evolutionary stages do not show any clear trend in Fig. 3, suggesting that similar underlying physical processes are responsible for high-J CO emission from all sources in the sample.

In summary, high-J CO emission is detected in high-mass clumps and correlates most strongly with clump bolometric luminosity. The line shapes show that high-velocity gas is most likely associated with the outflows.

3.2 Profile decomposition

We used mid-J (6 ≲ J ≲ 10) CO rare isotopologue lines to subtract the envelope component from the line profiles of CO 11−10 and CO 16−15. This way, we isolated the high-velocity emission associated with the line wings.

The emission in the line wings is characterised using a decomposition method which is described in detail in Appendix D. Briefly, the decomposition procedure aims to subtract the contribution from the envelope, as traced by rare isotopologue emission, resulting in the residual outflow component (Codella et al. 2004; van der Walt et al. 2007; de Villiers et al. 2014; Yang et al. 2018a). This method was initially used for kinematical studies of methanol masers, and subsequently adopted for low-J CO line profiles. Here, we used the version described in Yang et al. (2018a), which does not account for the opacity broadening because the high-J CO lines are likely optically thin. Rare isotopologue lines are used as a proxy for the envelope emission; here, depending on data availability and detection, we used the emission of the C¹⁸O9−8 line for eight sources, the ¹³CO 10−9 line for three sources, the C¹⁸O6−5 line for one source, and the ¹³CO6−5 line for one source (see Appendix D).

We identified line wing emission in the CO 11−10 line towards all sources except G13.66−0.6 (Table 3). The wings in the CO 16−15 line are seen only towards eight out of ten sources with the 3σ line detection. Properties and profiles of all wing emission are shown in Appendix D.

The ubiquity of line wings is consistent with previous detections of the outflows towards the same sources using lower-J lines of CO and SiO (Table 3). In particular, all sources from our sample show line wings in the CO 6−5 line (Navarete et al. 2019). The non-detection of the CO 11−10 line wing in G13.66−0.6 could be either due to a low signal-to-noise ratio (S/N; Fig. 1) or a lack of recent heating of the outflow gas due to shocks (Karska et al. 2013; Kristensen et al. 2017). The ¹³CO2−1 wings have only been seen towards G351.77−0.5, G12.81−0.2, and G14.19−0.2 (Yang et al. 2022) due to limited sensitivity, illustrating the difficulty in detecting line wings in rare CO isotopologues (see also, Stephens et al. 2018, 2019). Finally, SiO 2−1 have been observed towards our sources (Urquhart et al. 2019; Csengeri et al. 2016) and line wings have been detected in six of them. All the non-detections, in fact, show line wings in the high-J CO lines (Table 3), indicating that additional factors play a role in the excitation of SiO and CO lines.

Detecting outflows towards distant star-forming clumps is often hampered by confusion. Background and foreground galactic sources along the line of sight might contribute to the wing emission, which may result in false outflow detections. We note, however, that the high detection statistics (>60%) of line wings and the wings’ smooth shapes in our source sample are very unlikely to be explained by source confusion. The high-J CO emission is typically well confined to the regions with active star formation.

In conclusion, our decomposition method results in the estimate of line wing emission towards 12 and eight sources in the CO 11−10 and CO 16−15 lines, respectively.

Fig. 1 SOFIA line profiles of CO J = 11−10 (black), 16−15 (blue), and 13−12 (bottom right) transitions. All spectra were resampled to a common spectral resolution of 1.0 km s−1. Black vertical lines show values of Vlsr (see Table 1 ). Green horizontal lines show baselines.

Fig. 2 SOFIA/GREAT line profiles of the CO 11−10 and 16−15 lines as well as the CO 6−5 line. Source velocities (Vlsr) are shown with vertical lines. The lines were smoothed to a common bin of 1.0 km s−1.

Fig. 3 Line luminosities of CO 11−10 and 16−15, as a function of Mclump and Lbol IR-weak (24d) sources are shown with blue circles, IR-bright (IRb) sources with red triangles, and HII regions (HII) with green squares. The linear regression fit with Markov chain Monte Carlo is shown with dashed black lines and yellow shades. The linear log-log and Pearson correlation coefficients, r, are presented in each plot. Objects with self-absorption are shown with an upward arrow, indicating the lower limit for calculated luminosities.

3.3 CO line wing emission

Decomposition of the line profiles allowed us to quantify the amount of high-J CO emission in the line wings, and its contribution to the entire line profiles. Furthermore, the ratio of the two CO transitions can be studied as a function of gas velocity.

The fraction of emission in the line wings of the CO 11−10 transition ranges from ~29 to 73%, whereas the mean fraction of each evolutionary stage is ~50%, suggesting that there is no dependence with the source evolution (Fig. 4). The fraction of emission in CO 16−15 line wings is similar to the CO 11−10 transition, and ranges from ~28 to 76%. These results are consistent with the fraction of line wing emission measured towards two high-mass protostellar objects: AFGL 2591 in both CO 11−10 (~37%) and CO 16−15 (~34%) from van der Wiel et al. (2013), and W3 IRS5 in CO 10−9 (~50%) from San José-García et al. (2016).

The fraction of high-J CO emission has also been estimated for several high-mass YSOs by subtracting the envelope contribution from the total, unresolved line profiles. The Herschel/HIFI observations of rare isotopologues of CO were used to constrain models of CO (main isotopologue) emission arising from envelopes (e.g. Herpin et al. 2012, 2016; Karska et al. 2014; Jacq et al. 2016). In the case of NGC 7538 IRS1, 70–100% of velocity-unresolved CO J = 15−14 emission and 322% of CO J = 22−21 was attributed to the envelope (Karska et al. 2014). Thus, the contribution of the outflow component was predicted to increase with the rotational level of the CO line.

The increase in the relative contribution of the wing emission from J_up = 11 to 16 was indeed measured for four out of eight sources in our sample, for which outflow wings are detected in both CO transitions. The fraction of wing emission increases from ~42% (CO 11−10) to 50% (CO 16−15) for G34.26+0.15, from 57% to 68% for G351.16+0.7, from 69% to 76% for G351.44+0.7, and from 62% to 70% for G351.58−0.4. The increase is therefore not as sharp as for the CO 15−14 and CO 22−21 lines of NGC 7538 IRS1, but consistent with a rising contribution of wing emission in higher-J lines.

The amount of emission in the wings of higher-J CO lines allows us to study the gas excitation conditions in the outflowing gas. Assuming that emission in the line wings is optically thin and thermalised, the higher CO line ratios would correspond to higher gas kinetic temperatures, T_kin (see Sects. 3.4 and 3.5). Figure 5 shows the observed ratio of CO 16−15 and 11−10 lines in the red and blue wings as a function of the absolute offset from the source velocity. The ratio was calculated in steps of 1.0 km s^-1, avoiding the line centres (±5km s⁻¹ ), and presented for channels where the S/N is above 2.

The ratio of CO 16−15 and CO 11−10 increases as a function of velocity for at least a few sources, for example, the red wing of G12.81−0.2, G351.16+0.7, and G351.77−0.5, and the blue wing of G35.20−0.7 (Fig. 5). In most of those cases, the highest-velocity emission is stronger in CO 16−15 than in CO 11−10. Such trends are consistent with similar studies using CO 3−2,10− 9, and 16−15 towards a sample of low- to high-mass protostars (see Sect. 4.2).

In summary, we find a lack of correlation between the fraction of high- J CO integrated emission in the line wings and the clump evolutionary stage. Even so, the fraction increases with the CO rotational level in half of the sources. The ratio of the wing emission in the CO 16−15 and CO 11−10 lines increases with the offset velocity in several sources.

Fig. 4 Fraction of the integrated emission in the line wings of CO 11−10 (left) and CO 16−15 (right), as a function of Lbol /Mclump (Table 1). Dashed horizontal lines show the mean wing fraction at each evolutionary stage. The colour-coding is the same as in Fig. 3.

Fig. 5 Ratio of the line wing emission in CO 16−15 and 11−10 transitions as a function of the absolute velocity offset from the source velocity. The red-shifted emission is shown with red squares, and the blue-shifted emission with blue circles. The dashed horizontal line presents the level above which CO 16−15 is greater than CO 11−10.

3.4 Molecular excitation in LTE (full profile + wings)

Detection of at least two CO lines allowed us to determine the rotational temperature of the outflowing gas detected in the line wings under the assumption of LTE. For comparisons with previous studies with Herschel/PACS, the calculations were also performed for the velocity-integrated line profiles.

Emission line fluxes of CO 11−10 and CO 16−15 were used to calculate the number of emitting molecules, 𝒩_u, for each molecular transition as follows:

(1)

where L_CO refers to the line luminosity of the CO line at wavelength λ, A to the Einstein coefficient, c to the speed of light, and h to Planck’s constant. We note that for two sources, G12.81−0.2 and G351.25+0.7, additional observations of CO 13−12 are included. The number of emitting molecules, 𝒩_u, was used instead of column densities, because the size of the emitting region is unresolved.

The relation between 𝒩_u and the total number of emitting molecules, 𝒩_tot, follows the equation

(2)

where g_u is the statistical weight of the upper level, E_u the energy of the upper level, k_b the Boltzmann constant, T_rot the rotational temperature, and Q(T_rot) the partition function at the temperature T_rot·

The rotational temperature was calculated from slope b of the linear fit (y = ax + b) to the data in the natural logarithm units, T_rot = −1/a. The total number of emitting molecules, 𝒩_tot, was determined from the fit intercept b as follows:

(3)

Figure 6 shows example Boltzmann diagrams for G351.25+0.7 and G12.81−0.2, constructed using the velocity-integrated emission of CO (full profile). Table 4 shows T_rot and 𝒩_tot for all sources with at least two CO line detections separately for the integrated-profile emission and the line wings (see Sect. 3.2).

The two sources with three CO line detections are characterised by T_rot of ~170 K using the integrated line emission. The remaining sources show T_rot in the range from ~110 K to 200 K, with a mean of 152 K. Similar temperatures were obtained for the wing emission tracing outflow gas, with a mean T_rot of 167 K. While the wing emission is often responsible for the bulk of the total emission, hot core emission, with typical temperatures of ~ 100–200 K (Fontani et al. 2007; Taniguchi et al. 2023), might also contribute to the far-IR emission at the source velocity.

For G34.26+0.15, T_rot of ~150K is significantly lower than 365 ± 15 K obtained from Herschel/PACS (Karska et al. 2014). We note, however, that the latter temperature was obtained using CO lines with J_u from 14 to 30, sensitive to both ‘warm’ and ‘hot’ gas components (Karska et al. 2018). If CO transitions with J_u from 14 to 16 are used instead, T_rot of 244 ± 45 K is obtained (adopting values from Table C.1. in Karska et al. 2014). Even lower T_rot is expected when the CO 11−10 line, tracing a colder gas component, is used in the calculation, in line with the results obtained for G34.26+0.15. It is noteworthy that it is essential to have many CO transitions to determine all of the underlying physical conditions.

Finally, we note that the ratio of the total number of emitting molecules (𝒩_tot) in the line wings and the total line profile ranges from 40% to 79% (Table 4), consistent with the overall fraction of wing emission (Sect. 3.2). In absolute terms, log₁₀𝒩_tot ranges from 51.7 to 53.6, consistent with the average 52.4(0.1) ± 0.5 measured for high-mass protostars with Herschel/PACS (Karska et al. 2014).

Fig. 6 Rotational diagrams of CO for G12.81−0.2 and G351.25+0.7, which are based on the observations of full profile CO transitions with Ju of 11, 13, and 16. The natural logarithm of the number of emitting molecules from a level u, 𝒩u (dimensionless), divided by the degeneracy of the level, gu, is shown as a function of the upper level energy, Eu /kB, in Kelvins. Detections are shown as blue circles. Dashed orange lines show linear regression fits to the data; the resulting rotational temperatures are provided in the plots with the associated errors from the fit.

Fig. 7 CO excitation conditions from RADEX models versus observations. The plots present models at a different N(CO) of 1017 and 1018 cm−2. The models are presented with empty circles, and their colours correspond to a different hydrogen volume density, , between 103 and 107 cm−3. At each volume density level, five temperatures – 150, 250, 500, 1000, and 3000 K – were sampled. Observations assuming a beam filling factor of 1 are shown in crosses, while observations assuming a tiny source of 2″, which correspond to an extreme case of a small beam filling factor, are shown with triangles.

3.5 Molecular excitation in non-LTE (wings)

Due to relatively low densities in the regions of the ISM where outflows propagate, the LTE assumption may not hold. Non-LTE modelling is therefore necessary to determine the physical conditions responsible for the observed line emission. For this study, we used the well-established code RADEX (van der Tak et al. 2007) to estimate gas temperatures, densities, and CO column densities, which reproduce the observed line wing emission of three mid- and high-J CO lines: CO 6−5, 11−10, and 16−15.

We calculated model grids for a range of kinetic temperatures, T_kin, from 150 to 3000 K, H₂ number densities, , from 10³ to 10⁷ cm⁻³, and CO column densities, N(CO), of 10¹⁶,10¹⁷, and 10¹⁸ cm⁻². We assumed H₂ as the only collision partner, and a background temperature of 2.73 K. The linewidths of all lines were fixed at 19km s⁻¹, based on the observations of CO 6-5 (Appendix B and Navarete et al. 2019).

For comparisons of models with observations, we used peak intensities obtained from RADEX, since the wing emission does not follow a simple Gaussian; we also converted the observations from T_MB to T_r through T_r = T_MB/η. Because we did not spatially resolve the line emitting regions, we considered two cases during the computation of the beam filling factor: (i) the source that fills the entire beam (η = 1); and (ii) the source size of 2″ or η ~ 8 × 10⁻³− 5 × 10⁻², consistent with the size of a source in our sample, G34.26+0.15, which was measured from the Spitzer/IRAC 3.6 µm image. The spatial extent of CO 6−5 emission is ~4 times larger than the APEX/CHAMP⁺ beam, according to previous observations (Navarete et al. 2019).

Figure 7 shows the comparison of non-LTE radiative transfer models with line wing observations of high-J CO lines⁵ (Sect. 3.2). The ratio of CO 16−15 and CO 11−10 depends on both T_kin and , and shows a spread of 3 orders of magnitude. On the other hand, the intensity of CO 6−5 is most sensitive to the assumed N(CO), and increases by 2–3 orders of magnitude between 10¹⁶ and 10¹⁸ cm⁻². The impact of the assumed beam filling factor is almost negligible to the CO 16−15 / CO 11−10 ratio.

The models match the observations best for the assumed CO column densities of 10¹⁷ and 10¹⁸ cm⁻² (Fig. 7). The solutions for temperature and density are degenerate and can be split into two regimes: (i) a lower-density scenario with of 10³– 10⁴ cm⁻³ and T_kin of at least 1000 K, and (ii) a high-density, moderate-temperature scenario with of 10^5–10⁷ cm⁻³ and T_kin between 150 and 500 K. The ratio of high-J CO lines can be well reproduced for both considered filling factors in the scenario (ii); the best-matching source size is likely larger than 2″ but depends on the assumed column density. In scenario (i), the ratio of high-J CO lines can be reproduced for a small fraction of our sample assuming T_kin of 1000 K. Much higher temperatures would be required to match observations of the majority of targets. In general, the CO 6−5 peak intensity increases with gas density: for example, for of 10³ cm⁻³ and N(CO) of 10¹⁸ cm⁻², models would match the observations assuming the filling factor of 1, whereas of 10⁴ cm⁻³ and N(CO) of 10¹⁷−10¹⁸ cm⁻² point at smaller filling factors.

In conclusion, only scenario (ii) can explain the observations of all targets. The T_kin range in this scenario is also in better agreement with T_rot estimated under the LTE condition in Sect. 3.4, and it is consistent with detections of molecular species excited exclusively in high-density environments towards other high-mass clumps (e.g. van der Tak et al. 2013, 2019). On the other hand, scenario (i) requires temperatures in excess of 3000 K to explain the observed CO lines (E_up < 800 K) at more than half of our targets; such temperatures are too high even for the outflows from high-mass stars. Therefore, we prefer the high-density, moderate-temperature scenario to describe the physical conditions towards our source sample. We note, however, that our models constrain only the ranges of temperature and density, as we cannot fully break the degeneracy between different models.

To summarise, non-LTE radiative transfer models provide support to the LTE excitation of high- J CO emission in the highmass clumps. The best match with observations was obtained for gas densities of 10⁵–10⁷ cm⁻³, T_kin between 150 and 500K, and CO column densities of 10¹⁷ and 10¹⁸ cm⁻². Such conditions are consistent with CO excitation in outflows and this is discussed further in Sect. 4.2.

4 Discussion

High spectral resolution observations from SOFIA/GREAT allowed us to disentangle dynamical properties of high-J CO emission towards high-mass star-forming clumps. The excitation conditions have been studied in the high-velocity gas component assuming both LTE and non-LTE regimes, supporting the origin in moderate-temperature, high-density gas associated with the outflows (Sects. 3.4–3.5). Here, we discuss our results in the context of previous observations of high-mass protostars with Herschel and SOFIA.

4.1 High-J CO emission in high-mass clumps

The high-J (J ≳ 10) CO emission in high-mass star-forming regions has been attributed to gas cooling of several physical components, including (i) a warm, dense envelope of central protostars (Ceccarelli et al. 1996; Doty & Neufeld 1997), (ii) UV-irradiated outflow cavity walls (Bruderer et al. 2009; San José-García et al. 2016), (iii) currently shocked gas in the outflows (van der Wiel et al. 2013; Karska et al. 2014), and (iv) photodissociation regions (Lane et al. 1990; Ossenkopf et al. 2010, 2015; Stock et al. 2015). A similarity of CO to H₂O, both in spatial extent and line shapes, supported the scenario of shock excitation in similar layers composing the outflow cavity walls (see e.g. San José-García et al. 2016; Kristensen et al. 2017; van Dishoeck et al. 2021).

Broad line profiles of high-J CO lines provide a solid evidence of the outflow origin of a part of CO emission in high-mass clumps from the ATLASGAL survey (Sect. 3.1, see also San José-García et al. 2013; Indriolo et al. 2017). It is noteworthy that the fraction of CO emission in the line wings with respect to the total line emission is not sensitive to the evolutionary stage of the clumps (Sect. 3.3). In fact, a significant fraction of CO 11−10 is detected in the line wings of clumps at very early evolutionary stages (up to 76%, Sect. 3.3).

The signposts of outflow activity in the IR-weak clumps are in agreement with the ubiquitous detection of broad line wings in the SiO 2−1 line towards ATLASGAL sources spanning all evolutionary stages, including 25% of IR-quiet clumps (Csengeri et al. 2016). Molecular outflows are indeed also commonly detected towards 70 µm dark clumps using other tracers (Urquhart et al. 2022; Yang et al. 2022).

The integrated high-J CO emission shows a strong correlation with the clump bolometric luminosity (Sect. 3.1). The correlation extends even to low- and intermediate-mass protostars (Fig. 8), suggesting a similar physical mechanism operating over a few orders of magnitude different spatial scales. In deeply embedded low-mass objects, L_bol is dominated by accretion luminosity, which in turn is closely related to the amount of mass ejected in the outflows (Frank et al. 2014). Thus, the tight correlation of high-J CO emission with L_bol for ATLASGAL clumps suggests an equally high contribution of accretion luminosity in high-mass regions.

The velocity-resolved SOFIA spectra provide strong support for an origin of the bulk high- J CO emission in outflows, during all evolutionary stages of high-mass clumps. The correlation of CO line fluxes with bolometric luminosity suggest common physical conditions and processes leading to high-J CO emission from low- to high-mass star-forming regions.

Fig. 8 Velocity-integrated CO line luminosity of 11−10 and 16−15 transitions versus source bolometric luminosity from low- to high-mass star-forming regions. The dashed lines show a linear fit obtained using only the sources from our study, which are shown in blue empty circles. Blue stars show observations of other high-mass protostars from Karska et al. (2014), Indriolo et al. (2017), and Kaźmierczak-Barthel et al. (2014); orange squares present emission from intermediate-mass objects (Matuszak et al. 2015); and red triangles show data for Class 0 protostars (Kristensen et al. 2017).

4.2 Excitation conditions

Observations of multiple CO lines allow one to study gas excitation across various source properties and evolutionary stages. In combination with other far-IR lines, they also constrain the properties of shocks responsible for the emission in broad line wings of high- J CO lines.

The rotational temperatures in the high-velocity gas in the ATLASGAL clumps range from ~120K to 219 K, and are similar to the temperatures obtained from the full line profiles (Sect. 3.4). Five IR-bright clumps show a mean T_rot of 169 ± 30 K, whereas two HII clumps are characterised by T_rot of 147 ± 16 K. Thus, a possible decrease in the gas temperature in the outflows as the clumps evolve might be present, but for the sources in our sample the difference is not significant.

Rotational temperatures of ~200–210 K, consistent with our measurements, have been estimated in the line wings of the high-mass protostar DR21(OH) assuming LTE (Leurini et al. 2015). Non-LTE modelling of multiple CO lines indicated T_kin of 60–200 K in the outflow gas component of another high-mass source, AFGL 2591 (van der Wiel et al. 2013). In W3 IRS5, excitation temperatures of ~ 100–210 K were measured in a range of CO 10−9 and 3−2 lines’ velocities covered by an outflow in this region (from 5 to 20 km s⁻¹), with the highest temperatures corresponding to the highest velocities (San José-García et al. 2013). A similar trend of increasing gas temperature with velocity was also clearly detected in the outflow wing emission from ATLAS-GAL clumps observed with SOFIA/GREAT (Sect. 3.3), and in the SiO survey of ~430 clumps observed with the IRAM 30 m telescope (Csengeri et al. 2016).

Comparisons of gas excitation using high-J CO lines can be extended to a larger number of sources once the emission in the full line profiles is considered. The velocity-integrated Herschel/PACS detections of CO transitions from J_u of 14 to 30 towards ten high-mass protostars provided an average T_rot of ~300(23)±60K (Karska et al. 2014). Protostars with detections of higher-J lines were generally characterised by a higher T_rot, suggesting the possible presence of an additional T_rot ≳ 700 K gas component detected towards low-mass protostars that appears to be ‘hidden’ in their high-mass counterparts, possibly due to the a small beam filling factor of such emission and/or optically thick continuum emission (Manoj et al. 2013; Green et al. 2013; Karska et al. 2013, 2018). Clearly, any comparisons of T_rot should consider the similar J levels for their calculation (Sect. 3.4, and e.g. Neufeld 2012; Jiménez-Donaire et al. 2017; Yang et al. 2018b).

Table 5 compares T_rot measurements for several high-mass protostars with the data of the same or similar CO transitions to our SOFIA/GREAT survey (Sect. 3.4). All sources with at least three observed transitions show T_rot from 130 K to 220 K, consistent with the values determined for ATLASGAL clumps and hot cores. The relatively narrow range of T_rot is qualitatively similar to that of the universal ‘warm’, ~300 K gas component based on CO 14−13 to 25−24 transitions towards low-, intermediate-, and high-mass protostars (Karska et al. 2014, 2018; Matuszak et al. 2015; van Dishoeck et al. 2021). The CO 11−10 transition in low-mass protostars is typically associated with a ‘cool’ gas component with T_rot ~l·00K (e.g. Yang et al. 2018b), and its inclusion in the fit causes the lower values of T_rot (< 300 K).

The CO rotational temperature depends on the gas density and kinetic temperature, and it can be characterised with both (i) low-density, high-temperature (Neufeld 2012; Manoj et al. 2013; Yang et al. 2018b) and (ii) high-density, low-temperature regimes (Karska et al. 2013, 2018; Green et al. 2013; Kristensen et al. 2017). The first scenario requires cm⁻³ and T_kin ≳ 2000 K, and has the advantage of reproducing the positive curvature of the CO diagrams over a broad range of energy levels (Neufeld 2012). In contrast, the second scenario accounts for the similarity of J ≳ 14 CO and H₂O emission most evidently seen in low-mass protostars with high-densities required for H₂O excitation (Karska et al. 2013; van Dishoeck et al. 2021).

Non-LTE modelling of massive clumps provides strong support for a high-density scenario of CO excitation (Sect. 3.5). In the regime of moderate gas temperatures (T_kin from 150 to 500 K), gas densities of 10⁵−10⁷ cm⁻³ match the data best. Such physical conditions are fully consistent with the modelling of high-J CO and H₂O emission towards high-mass protostars (San José-García et al. 2016; van Dishoeck et al. 2021). They are also comparable to the physical conditions determined in the jet, terminal shock, and cavities of the intermediate-mass protostar Cep E (Lefloch et al. 2015).

The underlying mechanism behind the highly excited CO gas has been investigated for both Cep E and its high-mass counterpart Cep A. Detailed comparisons of CO, in combination with [O I] and OH, suggest that the origin is in dissociative or UV-irradiated shock models with pre-shock densities above 10⁵ cm⁻³ (Gusdorf et al. 2016, 2017). Assuming a compression factor of ~100, typical for dissociative shocks (Karska et al. 2013), such models would also be in agreement with radiative-transfer modelling for high-mass clumps (Sect. 3.5). However, a fraction of high-J CO emission detected at source velocity could also originate from the central hot core.

5 Conclusions

We have characterised the SOFIA/GREAT line profiles observed towards 13 high-mass protostars selected from the ATLASGAL survey, which significantly increases the number of high-mass objects that have velocity-resolved high-J CO lines. The velocity information enabled us to quantify the line components and the properties of their emitting sources. We summarise and draw the following conclusions:

• CO 11−10 emission was detected towards all of the sources, as early as the 24d stage. Ten out of 13 clumps also show a clear detection of CO 16−15. Additionally, CO 13−12 was detected towards two sources. The lines exhibit broad line wing emission typical for outflows from YSOs;

• We detected wing emission in the CO 11−10 line from 12 clumps and in the CO 16−15 line from eight clumps, implying that the highly excited CO lines originate in outflows. The wing fraction is similar for all clump evolutionary stages. On the other hand, we found no signatures of high-velocity gas (i.e. bullets) in the far-IR spectra;

• Under the LTE assumption, we found T_rot of 110–200 K for the entire line profiles and 120–220 K for the wing component. Such temperatures are in agreement with gas densities of 10⁵−10⁷ cm⁻³, moderate temperatures of 150K and 500K, and CO column densities of 10¹⁷ and 10¹⁸ cm⁻² obtained from the non-LTE models;

• Significant correlations between high-J CO emission and bolometric luminosities suggest similar underlying physical processes and conditions across all evolutionary stages of high-mass clumps. The correlation also extends to low-mass protostars, where high-J CO emission originate in outflow shocks, consistent with our study.

High angular resolution maps of high-mass clumps would be necessary to better characterise the physical structure of the regions with strong high-J CO emission and spatially disentangle outflows and hot cores (Goicoechea et al. 2015). The MIRI instrument on board the James Webb Space Telescope could pinpoint the spatial extent of shocked gas in high-mass star-forming clumps.

Acknowledgements

The authors thank the anonymous referee for detailed comments that have helped us improve this paper. We thank SOFIA/GREAT staff for collecting and reducing the data. We also thank Dr. Helmut Wiesemeyer for his support in reducing data from the SOFIA/4GREAT observations. A.K. acknowledges support from the Polish National Agency for Academic Exchange grant no. BPN/BEK/2021/1/00319/DEC/1. A.Y.Y. acknowledges support from the National Natural Science Foundation of China (NSFC) grants no. 11988101 and no. 12303031. This work is based on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA). SOFIA is jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NNA17BF53C, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 2002 to the University of Stuttgart. Herschel was an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. This publication is based on data acquired with the Atacama Pathfinder Experiment (APEX). APEX is a collaboration between the Max-Planck-Institut fur Radioastronomie, the European Southern Observatory, and the Onsala Space Observatory.

Appendix A Spatial extent of line emission

For four sources, the GREAT/LFA provided additional CO 16−15 spectra towards six positions that are at 31″.8 (Risacher et al. 2018) from the source centre. We detected no extended emission towards G34.40−0.2 and G35.20−0.7. On the other hand, G34.41+0.2 and G34.26+0.15 show a weak line emission at one offset position each (Fig. A.1). For G34.41+0.2, we found similar line intensity peaks at the central and extended positions. However, G34.26+0.15 shows ~5 times weaker emission at the offset positions than in the centre. At both objects, we observed no signs of line wing emission at the offset positions.

Fig. A.1 Spectra of the CO 16−15 emission arising from extended positions away from source centres. The offsets in right ascension and declination are shown in each spectrum. Black vertical lines mark the position of V1sr. Baselines are shown with green horizontal lines. The spectra were smoothed to 1.0 km s−1.

Appendix B Line parameters for emission at central position

Table B.1 presents peak and integrated intensities, and line luminosities of the high-J CO lines. Table B.2 shows results of the multiple Gaussian fit for sources with self-absorption.

Appendix C Full velocity spectra at central position

Figure C.1 and C.2 present the CO 11−10 and CO 16−15 spectra from Section 3.1, but covering a velocity range of ±150 km s⁻¹. The spectra illustrate a lack of EHV emission or any other significant high-velocity gas components towards the high-mass clumps. A 2-σ feature at -125 km s⁻¹ appears to be tentatively detected towards G351.77−0.5, which is known to drive multiple outflows associated with EHV gas detected in CO 2 − 1 and 6 − 5 (Leurini et al. 2009). It is noteworthy that the EHV gas component was exclusively detected towards the outflow positions of G351.77−0.5, and absent from the on-source spectra (Leurini et al. 2009), consistent with observations towards an intermediate-mass protostar Cep E (Gómez-Ruiz et al. 2012; Lefloch et al. 2015; Gusdorf et al. 2017). Thus, we do not consider this feature as real.

Fig. C.1 CO 11−10 spectra on a 300 km s−1 window, in which source velocities were shifted to 0 km s−1 and are marked by the solid vertical line. Horizontal dashed lines present the zero intensity level.

Fig. C.2 CO 16−15 spectra on a window of 300 km s−1. The source velocities were shifted to 0 km s−1. Vertical and horizontal lines are similar to those in Fig. C.1.

Appendix D Details of the profile decomposition method

Appendix D.1 Line wing emission in high-J CO lines

Assuming that the high- J CO lines are optically thin, the detailed steps to identify wing emission are as follows:

(1) Scale the isotopologue line that is used as a proxy for envelope emission (called proxy line hereafter) to the height of the high- J CO profiles (see Fig. D.1, top). Table D.2 presents the envelope tracer used for each source.

(2) Fit a Gaussian to the peak of the scaled proxy line (see Fig. D.1, middle) and use it as a model for the envelope emission. In some cases, proxy lines have non-Gaussian shapes with additional wing features. The Gaussian fit helps remove contamination from the wing emission. The fit is iterated following van der Walt et al. (2007), with high velocity channels being gradually removed until a reasonable fit is obtained with a significant improvement of reduced-χ².

(3) Shift the envelope model to the velocity of the high-J CO peaks. This step is skipped for the CO 11−10 line of G12.81−0.2 and G035.20−0.7 as their line peaks are tilted, suggesting that a part of the line profile is attenuated.

(4) Subtract the envelope model from the high-J CO profiles (see Fig. D.1, bottom), leaving wing emission in the residual. If the residual shows significant emission above three times the spectral noise, a wing detection is confirmed.

Properties of line wing emission are presented in Table D.3. From step (2), the Gaussian fit to the scaled proxy lines provides source radial velocities, V_lsr. Although V_lsr was previously determined with other tracers (e.g. C¹⁷O 3−2 in Giannetti et al. (2014)), low opacity proxy lines provide more reliable estimations as they probe the central gas envelope thoroughly.

The use of different envelope tracers in Table D.2 could pose a bias on the extracted outflow emission. To assess this issue, we decomposed the CO 11−10 line of G351.16+0.7 using all four envelope tracers (see Fig. D.3). Table D.1 shows that wings extracted by using the C¹⁸O lines and the ¹³CO 10−9 agree well with each other with differences less than 10%. Using the ¹³CO 6−5 line, however, returns much less wing emission. The mid-J ¹³CO line might have a higher opacity than the others, and thus having a wider width and smaller wing component at the CO 11−10 line. For this work, we used ¹³CO 6−5 for only one source. Therefore, the bias has been minimised. In fact, wing emission derived for that sources is a lower limit for the actual outflow emission, thereby the wing detection statistics is not affected.

To extract wing emission from self-absorbed lines, we need an adjusted approach (Fig D.2). The line peak is underestimated because the emission is absorbed by cooler gas in front of the

Fig. D.1 Wing extraction steps for an example high-J CO line from top to bottom.

bulk emission. To recover the line peak, we fitted a multiple-Gaussian profile to the line with the absorbed parts being masked out. The fitting process was adopted from Navarete et al. (2019). At first, a single Gaussian was fit to the full profiles. A second component was added if there was significant emission above 3σ in the residual. Subsequently, a third Gaussian would be added if the residual still showed emission. The fit was set to not adding further Gaussian component after the third one. Results of the multiple-Gaussian fit are shown in Table B.2. Adopting the fit peak, we could carry out line decomposition following the four steps above. In step (3), the envelope model was shifted to the velocity of the fit peak.

Appendix D.2 Wing emission in CO 6−5

Navarete et al. (2019) found broad components associated with outflows on CO 6−5 lines of our sources. For this study, we attempted to extract wing emission from the CO 6−5 line of seven sources in our sample (Table D.3) by applying the decomposition method in (Yang et al. 2022) (see Fig. D.8). This method is similar to the decomposition method we used for the high- J CO lines, but an additional step was added after step (2) to account for the opacity broadening effect of the optically thick CO 6−5.

Fig. D.2 Line wing extraction for a self-absorption line in which the line peak was determined from a multiple Gaussian fit.

Fig. D.3 CO 11−10 wing emission of G351.16+0.7 obtained by using different envelope tracers.

To correct for the opacity broadening, one needs to derive opacity in each frequency channel:

(D.1)

where τ₀ and v₀ are the central opacity and line frequency, and σ is the gas intrinsic velocity dispersion, namely, the FWHM of the Gaussian fit to the scaled proxy line. The CO 6−5 central opacity was estimated by multiplying the optical depth of ¹³CO 6−5 (Dat et al., in preparation) with the abundance ratio of¹²CO relative to ¹³CO, , which was calculated using the galactic radius relation derived in Giannetti et al. (2014) and the galactic distance. The broaden line, which can be used as a model for the envelope emission, then, could be obtained by the radiative

transfer equation:

(D.2)

where , and T_ex and T_bg are the excitation and background temperature, respectively. The term [J_v(T_ex)−J_v(T_bg)] was directly obtained from the peak of the Gaussian fit to the scaled proxy line.

The very high optical depths of the CO 6−5 line (up to 122) result in broad lines with flat peaks (see Fig. D.8, middle). Such line profiles are not seen in observed spectra, which suggests that our envelope models are not perfect. However, the wing emission retained after subtracting such envelope models from the full line profiles can still serve as a lower limit for the outflow emission. In total, we could determine CO 6−5 wing emission for all seven sources.

Fig. D.4 CO 11−10 wing emission (dashed red profile) extracted from the full line (solid black profile). The dotted magenta line is a model of envelope emission. The detection 3σ threshold is shown with a dashed blue line.

Fig. D.5 CO 13−12 wing emission (dashed red profile). The colour-coding and line styles are the same as in Fig. D.4.

Fig. D.6 CO 16−15 wing emission (dashed red profile). The colour-coding and line styles are the same as in Fig. D.4.

Fig. D.7 Ratio of line wing emission in the CO 16−15 and 11−10 transitions versus the absolute velocity offset from the source velocity. The red-shifted emission is shown with the red squares, while the blue-shifted emission is shown with the blue circles. The dashed horizontal line presents the level above which CO 16−15 is greater than CO 11−10.

Fig. D.8 Wing extraction steps for the CO 6−5 line of G351.16+0.7 from left to right (see Appendix D).

Notes. V_range is the velocity range of line wing. T_peak and S_int are peak and integrated intensity of the line wing emission, respectively.

References

All Tables

All Figures

	Fig. 1 SOFIA line profiles of CO J = 11−10 (black), 16−15 (blue), and 13−12 (bottom right) transitions. All spectra were resampled to a common spectral resolution of 1.0 km s−1. Black vertical lines show values of Vlsr (see Table 1 ). Green horizontal lines show baselines.
In the text

	Fig. 2 SOFIA/GREAT line profiles of the CO 11−10 and 16−15 lines as well as the CO 6−5 line. Source velocities (Vlsr) are shown with vertical lines. The lines were smoothed to a common bin of 1.0 km s−1.
In the text

Fig. 3 Line luminosities of CO 11−10 and 16−15, as a function of Mclump and Lbol IR-weak (24d) sources are shown with blue circles, IR-bright (IRb) sources with red triangles, and HII regions (HII) with green squares. The linear regression fit with Markov chain Monte Carlo is shown with dashed black lines and yellow shades. The linear log-log and Pearson correlation coefficients, r, are presented in each plot. Objects with self-absorption are shown with an upward arrow, indicating the lower limit for calculated luminosities.

In the text

	Fig. 4 Fraction of the integrated emission in the line wings of CO 11−10 (left) and CO 16−15 (right), as a function of Lbol /Mclump (Table 1). Dashed horizontal lines show the mean wing fraction at each evolutionary stage. The colour-coding is the same as in Fig. 3.
In the text

	Fig. 5 Ratio of the line wing emission in CO 16−15 and 11−10 transitions as a function of the absolute velocity offset from the source velocity. The red-shifted emission is shown with red squares, and the blue-shifted emission with blue circles. The dashed horizontal line presents the level above which CO 16−15 is greater than CO 11−10.
In the text

Fig. 6 Rotational diagrams of CO for G12.81−0.2 and G351.25+0.7, which are based on the observations of full profile CO transitions with Ju of 11, 13, and 16. The natural logarithm of the number of emitting molecules from a level u, 𝒩u (dimensionless), divided by the degeneracy of the level, gu, is shown as a function of the upper level energy, Eu /kB, in Kelvins. Detections are shown as blue circles. Dashed orange lines show linear regression fits to the data; the resulting rotational temperatures are provided in the plots with the associated errors from the fit.

In the text

Fig. 7 CO excitation conditions from RADEX models versus observations. The plots present models at a different N(CO) of 1017 and 1018 cm−2. The models are presented with empty circles, and their colours correspond to a different hydrogen volume density, , between 103 and 107 cm−3. At each volume density level, five temperatures – 150, 250, 500, 1000, and 3000 K – were sampled. Observations assuming a beam filling factor of 1 are shown in crosses, while observations assuming a tiny source of 2″, which correspond to an extreme case of a small beam filling factor, are shown with triangles.

In the text

Fig. 8 Velocity-integrated CO line luminosity of 11−10 and 16−15 transitions versus source bolometric luminosity from low- to high-mass star-forming regions. The dashed lines show a linear fit obtained using only the sources from our study, which are shown in blue empty circles. Blue stars show observations of other high-mass protostars from Karska et al. (2014), Indriolo et al. (2017), and Kaźmierczak-Barthel et al. (2014); orange squares present emission from intermediate-mass objects (Matuszak et al. 2015); and red triangles show data for Class 0 protostars (Kristensen et al. 2017).

In the text

	Fig. A.1 Spectra of the CO 16−15 emission arising from extended positions away from source centres. The offsets in right ascension and declination are shown in each spectrum. Black vertical lines mark the position of V1sr. Baselines are shown with green horizontal lines. The spectra were smoothed to 1.0 km s−1.
In the text

	Fig. C.1 CO 11−10 spectra on a 300 km s−1 window, in which source velocities were shifted to 0 km s−1 and are marked by the solid vertical line. Horizontal dashed lines present the zero intensity level.
In the text

	Fig. C.2 CO 16−15 spectra on a window of 300 km s−1. The source velocities were shifted to 0 km s−1. Vertical and horizontal lines are similar to those in Fig. C.1.
In the text

	Fig. D.1 Wing extraction steps for an example high-J CO line from top to bottom.
In the text

	Fig. D.2 Line wing extraction for a self-absorption line in which the line peak was determined from a multiple Gaussian fit.
In the text

	Fig. D.3 CO 11−10 wing emission of G351.16+0.7 obtained by using different envelope tracers.
In the text

	Fig. D.4 CO 11−10 wing emission (dashed red profile) extracted from the full line (solid black profile). The dotted magenta line is a model of envelope emission. The detection 3σ threshold is shown with a dashed blue line.
In the text

	Fig. D.5 CO 13−12 wing emission (dashed red profile). The colour-coding and line styles are the same as in Fig. D.4.
In the text

	Fig. D.6 CO 16−15 wing emission (dashed red profile). The colour-coding and line styles are the same as in Fig. D.4.
In the text

	Fig. D.7 Ratio of line wing emission in the CO 16−15 and 11−10 transitions versus the absolute velocity offset from the source velocity. The red-shifted emission is shown with the red squares, while the blue-shifted emission is shown with the blue circles. The dashed horizontal line presents the level above which CO 16−15 is greater than CO 11−10.
In the text

	Fig. D.8 Wing extraction steps for the CO 6−5 line of G351.16+0.7 from left to right (see Appendix D).
In the text