Velocity-resolved high-J CO emission from massive star-forming clumps

(Abridged) Context. Massive star formation is associated with energetic processes, which result in significant gas cooling via far-infrared (IR) lines. Velocity-resolved observations can constrain the kinematics of the gas, allowing the identification of the physical mechanisms responsible for gas heating. Aims. Our aim is to quantify far-infrared CO line emission toward high-mass star-forming regions, identify the high-velocity gas component associated with outflows, and estimate the physical conditions required for the excitation of the observed lines. Methods. Velocity-resolved SOFIA/GREAT spectra of 13 high-mass star forming clumps of various luminosities and evolutionary stages are studied using CO 11-10 and 16-15 lines. Results. All targets show strong high-J CO emission in the far-IR, characterized by broad line wings associated with outflows, thereby significantly increasing the sample of sources with velocity-resolved high-J CO spectra. The contribution of the emission in the line wings does not correlate with the envelope mass or evolutionary stage. Gas rotational temperatures cover a narrow range of 120-220 K for the line wings. The non-LTE radiative transfer models indicate gas densities of 1e5-1e7 cm-3 and N(CO) of 1e17- 1e18 cm-2, similar to physical conditions in deeply-embedded low- and high-mass protostars. The velocity-integrated CO line fluxes correlate with the bolometric luminosity over 7 orders of magnitude including data on the low-mass protostars, suggesting similar processes are responsible for the high-J CO excitation over a significant range of physical scales. Conclusions. Velocity-resolved line profiles allow the detection of outflows toward massive star-forming clumps spanning a broad range of evolutionary stages. The lack of clear evolutionary trends suggest that mass accretion and ejection prevail during the entire lifetime of star-forming clumps.


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).
⋆ tdhoang@mpifr-bonn.mpg.deOutflows, a ubiquitous phenomenon in both low and high mass 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/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) velocityresolved rotational lines of carbon monoxide (CO), which is the second most abundant molecule in the interstellar medium (CO/H 2 = 1.2 × 10 −4 , Frerking et al. 1982).The low-lying rotational levels of CO are easily collisionally-excited even at low densities and can readily be observed at millimeter 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 toward 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(Yang et al. , 2022)).
Observations of high-J CO (J ≳ 10) lines provide an opportunity to study denser and warmer parts of star-forming clumps and the ouflows 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(Karska et al. , 2014(Karska et al. , 2018;;van Dishoeck et al. 2021).The velocity-resolved profiles of high-J CO toward low-mass protostars revealed a significant contribution of high-velocity (v ∼ 20-30 km s −1 ) gas to the total far-IR line emission (San José-García et al. 2013;Yıldız et al. 2013), and similarity to the H 2 O emission likely arising from the same gas (San José-García et al. 2016;Kristensen et al. 2017).Single-pointing observations toward high-mass sources have also revealed broad, outflow wings in high-J CO line profiles, but 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 1 (GREAT, Risacher et al. 2018) onboard 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 highvelocity gas (EHV; up to ∼ 140 km s −1 ) tracing shocks associated with the jet and intermediate-to-high velocity gas ( from 50 to 100 km s −1 ) 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 toward two highmass sources, however, lacked the EHV component and revealed broad line wings extending up to ∼ 50 km s −1 (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. 2013Goicoechea et al. , 2015)).
In this paper, we use SOFIA/GREAT to quantify high-J CO emission toward 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 examine how the high-J CO emission varies as a function of clump properties and evolutionary stages.
The paper is organized as follows: Section 2 describes the source sample, observations with SOFIA, and the complementary CO observations with the APEX telescope and Herschel.
In Section 3, we present line profiles of high-J CO transitions (Section 3.1) and decompose the emission that belongs to the line wings (Section 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 (Section 3.3).Subsequently, we analyse the excitation of high-J CO lines using LTE and non-LTE approaches (Sections 3.4 and 3.5).Section 4 consists of the discussion of our results in the context of previous studies and Section 5 presents a summary and our conclusions.

Sample
All sources have been selected from the ATLASGAL survey covering 420 deg 2 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 catalog contains 5007 clumps spanning a wide range of masses (M clump ) and luminosities (L bol ), and divided into four evolutionary stages -quiescent, protostellar, young stellar objects, and H II regions (H II), see Urquhart et al. (2022).
For this work, we originally selected a representative sample of 20 sources grouped within 4 star-forming regions in the Galactic plane.Among them, 13 sources within 3 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 3 protostellar (24d), 7 young stellar object (IRb), and 3 H II regions (HII), with L bol from 1.6×10 3 to 4.6 × 10 5 L ⊙ and M clump from 1.6 × 10 2 to 2.3 × 10 3 M ⊙ (König et al. 2017;Urquhart et al. 2019Urquhart et al. , 2022)).
GREAT was a high resolution, dual-color spectrometer (R ≥ 10 7 ) 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.52THz and 1.80-1.90THz windows, respectively.In 2016, we combined the GREAT's L1 channel with the upGREAT Low Frequency Array (LFA) which covered the 1.83-2.07THz window in two polarizations.The 7-pixel hexagonal setup of the LFA provided spatial information about the line emission whereas each pixel had an FWHM beam size of 14.8 ′′ on the sky 2 .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 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 are 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.6where 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.0 km s −1 .
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 done in 2019 under  Freq.Notes.Molecular data adopted from the Leiden Atomic and Molecular Database (LAMDA, Schöier et al. 2005).∆v is the original spectral resolution of our data, before they are smoothed to a common size of 1.0 km s −1 .The CO 11-10 and CO 16-15 have two ∆v each, which correspond to observations in 2014 and 2016, respectively.
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.52THz and 2.49-2.69THz windows, were tuned to map the CO 13-12 and 22-21 transitions.The maps were scanned in 5 × 5 grids with centers 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 ∆v of 0.6 km s −1 .Data reduction for the CO 13-12 line was performed with the CLASS program, which is part of the GILDAS 3 software developed by the Institut de Radioastronomie Millimétrique (IRAM).
3 https://www.iram.fr/IRAMFR/GILDAS/A second order baseline was subtracted, and the spectra were also smoothed to an adequate 1.0 km s −1 .For this study, we extracted averaged spectra within a beam of 20 ′′ .

Additional observations and ancillary data
Additional single pointing observations of 13 CO 10-9 and C 18 O 9-8 were conducted with the Herschel-Heterodyne Instrument for the Far-Infrared (HIFI, de Graauw et al. 2010) onboard of the Herschel space telescope.Observations for 10 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 region with Herschel' program (San José-García et al. 2013;van Dishoeck et al. 2021).Data from the H and V polarizations 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 Fig. 1: SOFIA line profiles of CO J=11-10 (black), 16-15 (blue), 13-12 (bottom right) transitions.All spectra are resampled to a common spectral resolution of 1.0 km s −1 .Black vertical lines show values of V lsr (see Table 1).Green horizontal lines show baselines.of 0.96 and a main beam efficiency of 0.64 for the 13 CO 10-9 line and 0.74 for the C 18 O 9-8 lines, respectively.Finally, the spectra were smoothed to 1.0 km s −1 .Angular and original spectral resolutions are listed in Table .2.
We also use high spectral resolution 13 CO 6-5 and C 18 O 6-5 data from Dat et al. (in preparation) 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 12 m submillimeter telescope (APEX) (Güsten et al. 2006).The on-the-fly (OTF) scans resulted in datacubes of 80 ′′ × 80 ′′ with angular resolution of ∼9 ′′ .For comparisons with the higher-J CO observations, averaged spectra with an effective beam size of 20 ′′ around the Fig. 2: SOFIA/GREAT line profiles of the CO 11-10 and 16-15 lines as well as the CO 6-5 lines.Source velocities (V lsr ) are shown with vertical lines.The lines are smoothed to a common bin of 1.0 km s −1 .sources were extracted and then smoothed to 1.0 km s −1 for all three lines.

Line detections
In this section, we examine detection rates of high-J CO lines toward high-mass clumps from ATLASGAL, 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 toward 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 centers toward four sources, see Appendix A.
The CO 11-10 line is detected at 3σ or higher levels toward all sources, which span a broad range of evolutionary stages and have diverse properties.The CO 16-15 line, however, is firmly detected toward 10 out of 13 clumps; G13.66−0.6 and G34.41+0.2show only 2σ peaks and G14.19−0.2shows a nondetection.In addition, the CO 13-12 line was successfully observed and detected toward 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 4 (FWZP) of 45 km s −1 4 The FWZP is calculated following a procedure described in San José-García et al. (2016).We first resample the spectra to 3 km s −1 , and subsequently check the velocity of the channel where the line emission drops below 1σ.
is measured for the CO 11-10 line and 33 km s −1 for the CO 16-15 line (Appendix C).The broadest profile, with FWZP of 165 km s −1 , is seen toward 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 EHV gas component toward 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 toward 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 toward the same sources (Fig. 2).Selfabsorption features are seen in the CO 11-10 line profiles toward 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 shows also a sign of selfabsorption 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 ∼38 km s −1 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 ∼1 km s −1 from the source velocity obtained from the C 18 O 9-8 line (Fig. 1).The uncertainty of the Gaussian fit to the C 18 O line is smaller than 0.25 km s −1 , and thus cannot ac- We calculate CO line luminosities, L CO , as 4πD 2 F CO λ , where D is the distance to the source (Table 1) and F CO λ is the velocity integrated flux in W m −1 .The flux conversion from K km s −1 to W m −1 follows Equation 1 in Indriolo et al. (2017).Figure 3 shows the correlations between L CO and source properties (Table 1).The significance of the correlations is quantified by the Pearson correlation coefficient r, which depends also 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 , tracing primarily 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 Section 4).Noteworthy, clumps at different evolutionary stages do not show any clear trend in Figure .3, suggesting 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 most likely associated with the outflows.

Profile decomposition
We use 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 isolate 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 use the version described in Yang et al. (2018a) which do 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 18 O 9-8 line for eight sources, the 13 CO 10-9 line for three sources, the C 18 O 6-5 for one source, and the 13 CO 6-5 for one source (see Appendix D).
We identify line wing emission in the CO 11-10 line toward all sources except G13.66−0.6 (Table 3).The wings in the CO 16-15 line are seen only toward 8 out of 10 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 toward 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 the low S/N (Fig. 1) or a lack of recent heating of the outflow gas due to shocks (Karska et al. 2013;   3), indicating that additional factors play a role in the excitation of SiO and CO lines.
Detecting outflows toward distant star-forming clumps is often a 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 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 toward 12 and 8 sources in the CO 11-10 and CO 16-15 lines, respectively.

CO line wing emission
Decomposition of the line profiles allows 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 toward 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), The fraction of high-J CO emission has been also 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. 2012Herpin et al. , 2016;;Karska et al. 2014;Jacq et al. 2016).In case of NGC 7538 IRS1, 70-100% of velocity-unresolved CO J = 15-14 emission and 3-22% 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 CO line.
The increase of the relative contribution of the wing emission from J up = 11 to 16 is indeed measured for 4 out of 8 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 thermalized, the higher CO line ratios would correspond to higher gas kinetic temperatures, T kin (see Sections 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 absolute offset from the source velocity.The ratio is calculated in steps of 1.0 km s −1 , avoiding the line centers (±5 km s −1 ), and presented for channels where signal-to-noise 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, e.g., 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 the CO 16-15 than in the CO 11-10.Such trends are consistent with similar studies using CO 3-2, 10-9, and 16-15 toward a sample of low-to high-mass protostars (see Section 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.Yet, 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 velocity in several sources.

Molecular excitation in LTE (full profile + wings)
Detection of at least two CO lines allows 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 are also performed for the velocity-integrated line profiles.
Emission line fluxes of CO 11-10 and CO 16-15 are used to calculate the number of emitting molecules, N u , for each molecular transition as: where L CO refers to the line luminosity of CO line at wavelength λ, A to the Einstein coefficient, c to the speed of light, and h to the Planck's constant.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, N u , is used instead of column densities, because the size of the emitting region is unresolved.The relation between N u and the total number of emitting molecules, N tot , follows the equation: 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 is calculated from the 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, N tot , is determined from the fit intercept b as: (3) Figure 6 shows example Boltzmann diagrams for G351.25+0.7 and G12.81−0.2,constructed using the velocityintegrated emission of CO (full profile).Table 4 shows T rot and N tot for all sources with at least two CO line detections, separately for the integrated-profile emission and the line wings (see Section 3.2).
The two sources with three CO line detections are characterized 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 are 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 source velocity.
For G34.26+0.15,T rot of ∼150 K 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, tracing colder gas component, is used in the calculation, in line with results obtained for G34.26+0.15.Noteworthy, it is essential to have many CO transitions to determine all the underlying physical conditions.
Finally, we note that the ratio of the total number of emitting molecules (N 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 (Section 3.2).In absolute terms, log 10 N 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).

Molecular excitation in non-LTE (wings)
Due to relatively low densities in the regions of 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.Here, we use 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 midand high-J CO lines:  We calculated model grids for a range of kinetic temperatures, T kin , from 150 to 3000 K, H 2 number densities, n H 2 , from 10 3 to 10 7 cm −3 , and CO column densities, N(CO), of 10 16 , 10 17 , Fig. 7: CO excitation conditions from RADEX models versus observations.The plots present models at different N(CO) of 10 17 and 10 18 cm −2 .The models are presented in empty circles, and their colors correspond to different hydrogen volume density, n H 2 , between 10 3 and 10 7 cm −3 .At each volume density level, four temperatures: 150, 250, 500, 1000, and 3000 K are 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 small beam filling factor, are shown in triangles.and 10 18 cm −2 .We assumed H 2 as the only collision partner, and a background temperature of 2.73 K.The linewidths of all lines were fixed at 19 km s −1 , 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 have also converted the observations from T MB to T r through T r = T MB /η.Because we do 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); (ii) the source size of 2 ′′ or η ∼ 8 × 10 −3 -5 × 10 −2 , consistent with 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 lines5 (Section 3.2).The ratio of CO 16-15 and CO 11-10 depends on both T kin and n H 2 , 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 16 and 10 18 cm −2 .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 17 and 10 18 cm −2 (Fig. 7).The solutions for temperature and density are degenerate and can be split into two regimes: (i) lower-density with n H 2 of 10 3 -10 4 cm −3 and T kin of at least 1000 K, and (ii) high-density, moderate-temperature scenario with n H 2 of 10 5 -10 7 cm −3 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 bestmatching 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 n H 2 of 10 3 cm −3 and N(CO) of 10 18 cm −2 , models would match the observations assuming the filling factor of 1, whereas n H 2 of 10 4 cm −3 and N(CO) of 10 17 -10 18 cm −2 , 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 Section 3.4, and consistent with detections of molecular species excited exclusively in high-density environments toward other high-mass clumps (e.g., van der Tak et al. 2013van der Tak et al. , 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 highdensity, moderate-temperature scenario to describe the physical conditions toward 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 summarize, 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 is obtained for gas densities of 10 5 -10 7 cm −3 , T kin between 150 and 500 K, and CO column densities of 10 17 and 10 18 cm −2 .Such conditions are consistent with CO excitation in outflows and will be discussed further in Section 4.2.

Discussion
High spectral resolution observations from SOFIA/GREAT allow us to disentangle dynamical properties of high-J CO emission toward 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 .Here, we discuss our results in the context of previous observations of high-mass protostars with Herschel and SOFIA.
Broad line profiles of high-J CO lines provide a solid evidence of the outflow origin of a part of CO emission in highmass clumps from the ATLASGAL survey (Section 3.1, see also San José-García et al. 2013;Indriolo et al. 2017).Noteworthy, 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 (Section 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%, Section 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 toward ATLASGAL sources spanning all evolutionary stages, including 25% of infrared-quiet clumps (Csengeri et al. 2016).Indeed, molecular outflows are also commonly detected toward 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 (Section 3.1).The correlation extends even to low-and intermediate-mass protostars (Figure 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 ATLAS-GAL 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.

Excitation conditions
Observations of multiple CO lines allow 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 ∼120 K to 219 K, and are similar to the temperatures obtained from the full line profiles (Section 3.4).Five IR-bright clumps show mean T rot of 169 ± 30 K, whereas two H ii clumps are characterized by T rot of 147 ± 16 K.Thus, a possible decrease of gas temperature in the outflows as the clumps evolve might be present, but for the sources in our sample the difference is not significant.Notes.Rotational temperatures obtained from the integrated line profiles of CO transitions with J u of 11, 13, and 16. (a) Calculation using CO lines with J u from 14 to 17. (b) Calculation using CO lines with J u from 14 to 16. (c) Based only on CO 11-10 and CO 16-15 detections, see Section 3.3.
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 modeling 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 the CO 10-9 and 3-2 lines' velocities range covered by an outflow in this region (from 5 to 20 km s −1 ), 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 is also clearly detected in the outflow wing emission from ATLASGAL clumps observed with SOFIA/GREAT (Section 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 toward 10 high-mass protostars provided an average T rot of ∼300(23)±60 K (Karska et al. 2014).Protostars with detections of higher-J lines were generally characterized by higher-T rot , suggesting the possible presence of an additional T rot ≳ 700 K gas component detected toward 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. 2013Karska et al. , 2018)).Clearly, any comparisons of T rot should consider the similar J−levels for their calculation (Section 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 (Section 3.4).All sources with at least 3 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 toward low-, intermediate-, and high-mass protostars (Karska et al. 2014;Matuszak et al. 2015;Karska et al. 2018;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 ∼100 K (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 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(Karska et al. , 2018;;Green et al. 2013;Kristensen et al. 2017).The first scenario requires n H 2 ∼10 3 cm −3 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 2 O emission most evidently seen in low-mass protostars, with high-densities required for H 2 O excitation (Karska et al. 2013;van Dishoeck et al. 2021).
Non-LTE modeling of massive clumps provides strong support for a high-density scenario of CO excitation (Section 3.5).In the regime of moderate gas temperatures (T kin from 150 to 500 K), gas densities of 10 5 -10 7 cm −3 match the data best.Such physical conditions are fully-consistent with the modeling of high-J CO and H 2 O emission toward 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, suggests the origin in dissociative or UV-irradiated shock models with pre-shock densities above 10 5 cm −3 (Gusdorf et al. 2016(Gusdorf et al. , 2017)).Assuming a compression factor of ∼100, typical for dissociative shocks (Karska et al. 2013), such models would be also in agreement with radiative-transfer modeling for high-mass clumps (Section 3.5).However, a fraction of high-J CO emission detected at source velocity could also originate from the central hot core.

Conclusions
We have characterized the SOFIA/GREAT line profiles observed toward 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 enables to quantify the line components and the properties of their emitting sources.We summarise and draw the following conclusions: -CO 11-10 emission is detected toward all the sources, as early as the 24d stage.10 out of 13 clumps also show a clear detection of CO 16-15.Additionally, CO 13-12 is detected toward two sources.The lines exhibit broad line wing emission typical for outflows from YSOs.-We detect wing emission in the CO 11-10 line from 12 clumps and in the CO 16-15 line from 8 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 find no signatures of high-velocity gas (i.e., bullets) in the far-IR spectra.-Under the LTE assumption, we find 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 5 -10 7 cm −3 , moderate temperatures of 150 K and 500 K, and CO column densities of 10 17 and 10 18 cm −2 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 extends also to low-mass protostars, where high-J CO 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.

Fig. 3 :
Fig. 3: Line luminosities of CO 11-10 and 16-15, as a function of M clump and L bol .IR-weak (24d) sources are shown in blue circles, IR-bright (IRb) sources in red triangles, and HII regions (HII) in green squares.Linear regression fit with Markov chain Monte Carlo is shown in dashed black lines and yellow shades.The linear log-log and Pearson correlation coefficients, r, are presented on each plot.Objects with self-absorption are shown with an upward arrow, indicating the lower limit for calculated luminosities.
Fig. 4: The fraction of integrated emission in the line wings of CO 11-10 (left) and CO 16-15 (right), as a function of L bol /M clump (Table1).Dashed horizontal lines show the mean wing fraction at each evolutionary stage.The color-coding is the same as in Fig 3.

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

Fig. 6 :
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 J u of 11, 13, and 16.The natural logarithm of the number of emitting molecules from a level u, N u (dimensionless), divided by the degeneracy of the level, g u , is shown as a function of the upper level energy, E u /k B , in Kelvins.Detections are shown as blue circles.Dashed orange lines show linear regression fits to the data; the resulting rotational temperatures are provided on the plots with the associated errors from the fit.

Fig. 8 :
Fig. 8: Velocity-integrated CO line luminosity of 11-10 and 16-15 transitions versus source bolometric luminosity from lowto 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); 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).

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

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

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

Table 1 :
Catalog of source properties

Table 2 :
Overview of the observations

Table 3 :
Detection of line wing emission in high-J CO lines and their comparison to prior outflow detections ✓Notes.The symbols refer to: '-' non-detection, ✓ detections, '...' lack of data or unable to detect.

Table 4 :
CO rotational excitation for both the integrated line profiles and line wings only assuming LTE

Table 5 :
CO rotational excitation determined from integrated line profiles toward high-mass objects