© H. Beuther et al. 2022

1 Introduction

Stellar feedback in the form of radiation, winds, and supernova explosions can have positive as well as negative impacts on their environment. Positive feedback here is defined as inducing new star formation processes, whereas negative feedback is meant as the destruction process of the natal cloud and its limiting of further star formation. Feedback is important for the determination of the physical processes within individual clouds but also crucial in an integral sense for whole galaxies because it likely determines the general star formation efficiency (e.g., Matzner 2002; Elmegreen 2011; Hopkins et al. 2014; Geen et al. 2016, 2018; Kim et al. 2018).

With the goal to understand the impact of feedback processes onto the environment for individual resolved cases, the SOFIA (Stratospheric Observatory for Infrared Astronomy) legacy program FEEDBACK¹ targets 11 galactic high-mass star-forming regions and H II regions in the ionized atomic carbon fine-structure line [CII] at 158 μm (1.9 THz) and the atomic oxygen line [OI] at 63 μm (4.7 THz). The ionized carbon line [CII] is of particular importance because it allows us to study the kinematics of the ionized gas as well as the interfaces with the molecular clouds. The [CII] line is a direct tracer of bubble or stellar feedback kinematics (e.g., Pabst et al. 2020). Furthermore, the [CII] line is one of the dominant cooling lines in the interstellar medium at low to intermediate densities and UV fields (e.g., Hollenbach & Tielens 1997; Röllig & Ossenkopf 2013). The FEEDBACK program started in spring 2019 and is still ongoing. More details about this SOFIA legacy program are provided in Schneider et al. (2020). One outcome of the FEEDBACK project so far is the detection of expanding shells in [CII] emission in various bubble-shaped or bipolar H II regions, namely, in RCW120 (Luisi et al. 2021; Kabanovic et al. 2022), RCW49 (Tiwari et al. 2021), RCW36 (Bonne et al. 2022), and RCW79 (Zavagno et al., in prep.). These shells were first detected in Orion A (Pabst et al. 2019)and appear to be a ubiquitous phenomena. While the regions RCW49 and RCW120 reveal second -generation star formation within these shells (Tiwari et al. 2021; Luisi et al. 2021), the expanding veil nebula in Orion shows little evidence of dense clumps that could lead to star formation (Pabst et al. 2019; Goicoechea et al. 2020).

While an analysis of the whole sample will follow when all data have been collected, in the following, we concentrate on the prototypicalH II region NGC 7538, shown at optical and cm continuum wavelengths in Fig. 1. The region is located at a distance of 2.65 kpc (Moscadelli et al. 2009) and has a physical extent of several parsec. The two brightest sources exciting the H II region are IRS5 and IRS6, with spectral types of O9 and O3, respectively (Puga et al. 2010, marked in Fig. 1). The H II region is associated with several sites of active star formation, particularly in the south with the regions IRS1, S, and IRS9 (marked in Fig. 1). Toward these southern star-forming regions, the H II region appears to be sharply bounded, whereas in the northeast, the region shows diffuse emission extending well beyond the photon-dominated region (PDR) shell (Luisi et al. 2016). In a recent study combining CO(3–2) and [CII] data, Sandell et al. (2020) characterized a large north-south outflow emanating from the southern star-forming core IRS1 with blue-shifted ionized carbon [CII] and molecular CO emission north of IRS1 in the vicinity of the H II region. Furthermore, Townsley et al. (2018) have shown that diffuse X-ray emission is found toward the H II region, spreading even beyond the extend of the H II region as determined by, for instance, the cm continuum emission.

More generally speaking, the NGC 7538 H II region is embedded in a larger molecular cloud complex (e.g., Ungerechts et al. 2000; Fallscheer et al. 2013). Figure 2 gives an overview of various observed tracers. The cm continuum and 8μm emission appears to stem dominantly from the H II region. We should keep in mind that 8 μm emission is usually emitted in the cavity walls of the H II region, where the poly-cyclic aromatic hydrocarbons (PAH) survive. In contrast to these features, the dense gas as traced by the Herschel dust continuum or CO(3–2) emission is found more in the outskirts of the H II region, in particular toward the south, where the well-know star-forming regions IRS1, S, and IRS9 are located. As discussed in more detail below, the new SOFIA [CII] data combine different worlds: tracing the PDR and also showing strong emission features toward the active star-forming region IRS1 and a bar-like region at the southeastern interface of the H II region with the dense molecular cloud. Hence, the velocity-resolved [CII] data are the ideal probe for studying the feedback processes and the kinematic gas properties from the evolving H II region on its neighboring molecular cloud.

The specific topics we address in this study consider the impact from the expanding H II region onto the environmental dense gas. We seek to identify expanding [CII] shells or other dynamic features, either triggering or preventing new star formation processes, as well as identifying layered structures from a photon-dominated region (PDR). We also consider the carbon budget between ionized C⁺, atomic C⁰, and molecular CO gas specifically at the H II region and molecular cloud interface.

Fig. 1 Optical and cm continuum data of the NGC 7538 region, encompassing the field of view of our observations. The color scale shows the optical DSS2 (Digital Sky Survey v2) obtained with the red filter. The contours outline the 1.4 GHz cm continuum emission at 1′ angular resolution from the Canadian Galactic Plane Survey (CGPS, Taylor et al. 2003). The two stars in the center of the H II region mark the positions of IRS5 and IRS6 from Puga et al. (2010). The three five-pointed stars directly below the H II regions show the positions of the active sites of star formation IRS1, S, and IRS9.

2 Observations and data

2.1 SOFIA observations

The region was mapped with the dual-frequency array heterodyne receiver upGREAT² (Risacher et al. 2018) in the velocity-resolved [CII] fine-structure line at 158 μm/1.9 THz³. The [OI] 63 μm line was observed in parallel. Fast Fourier Transform Spectrometers (FFTS) with 4 GHz instantaneous bandwidth (Klein et al. 2012) provided a nominal spectral resolution for the [CII] line of 0.244 MHz or 0.04 km s⁻¹. The observations were done on two consecutive days during Cycle 7, on December 12 and 13, 2019, with an average T_sys single side-band (SSB) of 4486 and 4580 K and an average precipitable water vapor of 6.2 and 3.6 μm respectively per day. The data were calibrated to the main beam brightness temperature intensity scale, T_mb, with the kalibrate task (Guan et al. 2012) that is part of the standard GREAT pipeline, with a forward efficiency of 0.97 and an average beam efficiency of 0.66.

Mapping was conducted in an on-the-fly (OTF) mode where the entire NGC 7538 map was split into 4 tiles of roughly 436″ on each side. During scanning the data were better than Nyquist-sampled with a dump at every 5.2″. The final map size is 14.5′× 14.5′ = 11.2 pc × 11.2 pc. Each tile should typically be observed four times with different scanning angles to reduce potential striping effects. In the case of NGC 7538, three tiles had been mapped to a deeper sensitivity, whereas for the fourth northwestern tile, full coverage has not been achieved thus far. This results in slightly increased noise in the top-right tile of the map and residual striping apparent in Fig. 2. The data cubes used here have been resampled at a velocity resolution of 0.5 km s⁻¹ to increase the signal-to-noise ratio. Our final velocity-resampled and grided data-cube has an angular resolution of 15.8″ (corresponding to ~0.2 pc linear spatial resolution). The 1σ rms over the well-covered 3/4 of the map are ~0.9 K per 0.5 km s⁻¹ channel on a T_mb scale. The higher 1σ rms in the northwestern tile is ~1.9 K per 0.5 km s⁻¹ channel.

Furthermore, slightly smaller maps of the region were observed in the atomic carbon fine-structure line [CI] at 492 GHz and three CO transitions with J = 8−7, J = 11−10, and J = 22−21. Here, we focus on the lower-density tracer [CI] as well as the CO(8–7) transition at 921.7997 GHz. The [OI] 63 μm and higher-J CO lines tracing higher densities and temperatures will be analyzed in forthcoming studies. The angular resolution of the final data cubes for the [CI] and CO(8–7) lines are 63″ and 42″, respectively. The 1σ rms values for the [CI] and CO(8–7) lines at 0.5 km s⁻¹ are 0.7 and 1.4 K, respectively.

Fig. 2 Integrated intensity molecular and atomic line maps as well as continuum emission are presented in color scale toward NGC 7538. The specific lines, integration ranges, and continuum bands are labeled in each panel. The corresponding angular resolution is shown at the bottom-left of the panels. The contours in most panels as labeled show the [CII] emission at only a few selected contour levels to highlight similarities and differences. The contours in the bottom-middle panel outline the 1.4 GHz cm continuum emission from the CGPS. The three five-pointed stars mark the positions of the embedded sources NGC 7538IRS1, NGC 7538S, and NGC 7538IRS9. The two stars in the middle show the positions of the two main exciting sources of the H II region NGC 7538IRS5 and NGC 7538IRS6. The bar-like PDR as well as a a linear scale bar are presented in the top-left panel.

2.2 Complementary archival data

2.2.1 Optical, mid-infrared, far-infrared, and radio data

We used the Digital Sky Survey version 2 (DSS2) to obtain an optical image of the region in the red filter. The DSS2 digital images are based on photographic images obtained using the Oschin Schmidt Telescope of Palomar Mountain and the UK Schmidt Telescope. The plates were processed into the present compressed digital form with the permission of the institutions.

The 1.4 GHz radio continuum and the 21 cm HI data come from the Canadian Galactic Plane Survey (CGPS, Taylor et al. 2003). The angular resolution of these data products is ~ 1′. The 1σ rms of the continuum and line data are~0.5 and 3 K in 0.82 km s⁻¹ channels, respectively.The 8 μm photometric emission mid-infrared data were obtained with the IRAC camera (Fazio et al. 2004) on board the Spitzer satellite (Werner et al. 2004) at an angular resolution of ~2″. The 70 μm far-infrared data are observed with the PACS camera (Poglitsch et al. 2010) on board the Herschel satellite mission (Pilbratt et al. 2010) in the framework of the HOBYS key program (Motte et al. 2010). The angular resolution is ~ 5.6″. The H_α data were observed in the framework of “The INT Photometric H_α Survey of the Northern Galactic Plane” (IPHAS) at an angular resolution of ~ 1.1″ (Drew et al. 2005).

2.2.2 Herschel column density and CO(3–2) data

The Herschel gas column density map we use in this paper is a higher angular-resolution version (18″) of the one presented in Fallscheer et al. (2013). The map was produced by the HOBYS (Herschel imaging survey of OB Young Stellar objects, Motte et al. 2010) consortium and will soon be publicly provided. The column density map was determined following the HOBYS procedure, namely, from a pixel-to-pixel spectral energy distribution (SED) greybody fit to the 160, 250, 350, and 500 μm wavelength observations. There were no 100 μm maps available and the 70 μm data are not included because the interest lies in the distribution of the cold molecular gas. For the SED fit, we fixed the specific dust opacity per unit mass (dust+gas) approximated by the power law with β_d = 2 and left the dust temperature and column density (atomic and molecular gas, in the following we refer to that as gas column density) as free parameters (see, e.g., Russeil et al. 2013 for more details). The procedure explaining how high angular resolution maps were obtained is described in detail in Appendix A of Palmeirim et al. (2013). The concept is to employ a multi-scale decomposition of the flux maps and assume a constant line-of-sight temperature. The final map at 18″ angular resolution is constructed from the difference maps of the convolved SPIRE maps (at 500, 350, and 250 μm) and the temperature information from the color temperature derived from the 160 to 250 μm ratio.

The CO(3–2) data were obtained with the James Clark Maxwell Telescope (JCMT) as part of a large map of the NGC 7538 region first presented in Fallscheer et al. (2013) and later in Sandell et al. (2020). Following Fallscheer et al. (2013), the typical 1σ rms in 0.42 km s⁻¹ channels is 0.6 K on a scale. Using a beam efficiency of 0.64⁴, the data areconverted to main beam brightness temperature for the quantitative analysis. The beam size of these data is ~ 15″.

3 Results

3.1 Spatial structures

As indicated in the introduction, the [CII] emission in the NGC 7538 complex traces the large-scale PDR as well as the PDR associated with denser star-forming molecular gas at the edge of the H II region. To be more specific, Fig. 2 shows the integrated emission of several tracers and, in particular, the morphology of the [CII], the 8 and 70 μm emission is very similar. The ionized gas as traced by the H_α and 1.4 GHz emission is confined to the real H II region, and the [CII], 8 and 70 μm emission are partly wrapped around it (bottom-right panel of Fig. 2). While the [CII] line can form in PDRs as well as molecular media, diffuse 8 μm emission typically stems from UV-pumped infrared fluorescence from polycyclic aromatic hydrocarbon (PAH) molecules within photon-dominated regions (PDRs). The dense gas tracers (gas column densities from dust continuum and CO(3–2), bottom panels in Fig. 2) are emitted dominantly in the south and south-west of the region. We note that toward the largest part of the [CII] emission, the CO(3–2) line does not show significant self-absorption features (see spectra discussed below). The associated emission toward the young embedded star-forming region IRS1 as well as the bar-like PDR features are visible in all tracers shown in Fig. 2. In addition to this, some dense gas structures can also be found in the northeastern part of NGC 7538.

To show spatial similarities and differences in a bit more detail, Fig. 3 presents a compilation of different pairs of tracers always as contour maps overlaid on color-maps. We added here the higher excited CO(8–7) line (E_u ∕k = 199 K compared to E_u∕k = 33.2 K for the (3–2) transition) as well as the atomic carbon fine structure line [CI] at 492 GHz. The top-left and top-middle panels confirm that the [CII] and 8 μm emission agree very well and appear to at least partially wrap around the ionized gas (see also bottom-right panel of Fig. 2). A comparison of the two CO transitions in the top-right panel shows that in the southwestern bar-like emission structure, the higher excited CO(8–7) transition emits closer to the two exciting sources at the center of the H II region than the lower excited CO(3–2) transition. Similar layered structures within the bar can be identified between atomic and ionized carbon ([CI] and [CII], bottom-left panel of Fig. 3) as well as between theCO(8–7) and [CI] maps. In comparison to those layered structures, the spatial morphologies of the higher excited CO(8–7) and ionized carbon [CII] lines agree much better in the bar. This is a typical PDR layered structure in [CII]/[CI]/CO as well as in CO(8–7) and (3–2). We come back to this in Sect. 4.2.

3.2 Kinematic structures

To get a first proxy of the dynamics, Fig. 4 presents the integrated [CII] emission as well as the first- and second-moment maps of the [CII] data (intensity-weighted peak velocities and velocity dispersion). The first-moment map (middle panel of Fig. 4) reveals a general velocity gradient approximately from east to west over the entire extend of the H II region. On top of the overall velocity gradient, one can identify several blue-shifted emission regions in the middle of the map. In comparison, the second-moment or velocity dispersion map shows over large parts of the map rather uniformly low values below 5 km s⁻¹. However, ontop of that, we can clearly identify an almost ring-like structure with larger velocity dispersion on the order of 10 km s⁻¹. As discussed below, these regions of apparent large velocity dispersion are largely caused by multiple velocity components in the [CII] spectra.

For a more detailed analysis of the kinematics, it is crucial to look at individual spectra and we selected several positions in the region for further analysis. We labeled them as follows and they are all marked in Fig. 4:

• IRS5: one of the main exciting source of the H II region

• IRS1: the densest star-forming core

• MOB: a position in the middle of the bar-like PDR

• BlPN: blue peak in the north in the first-moment map

• BlPS: blue peak in the south of the first-moment map

• RPC: red peak in the center of the first-moment map

• RPE: red peak in the east of the first-moment map

• BrPN: broad emission in the north in the first-moment map

• BrPS: broad emission in the south in the first-moment map.

Figure 5 presents the corresponding [CII] spectra as well as the molecular and atomic counterparts in the CO(3–2) and [CI] lines⁵. Where possible, we fitted Gaussians to the [CII] and CO(3–2) spectra, as well as two components at several positions. While self-absorption may explain smaller dips in some spectra, the two velocity components fitted here are typically that far separated in velocity space that they should indeed be separate components. The corresponding full width half maximum values (FWHM) in km s⁻¹ are also shown in Fig. 5. Interestingly, the only nearly simple, single-Gaussian spectrum is identified toward the position in the middle of the bar (MOB, top-right panel in Fig. 5). We discuss that structure further in Sect. 4.2. The spectra toward the infrared sources IRS5 and IRS1 as well as the red-shifted peak in the east of the region (RPE) show some velocity structure, but no particularly prominent features.

This is very different in the remaining five spectra toward the other blue- and red-shifted positions in the center of the maps (BlPN, BlPS, RPC) as well as toward the positions with a particularly broad velocity dispersion (BrPN, BrPS, and Figs. 4 and 5). While the main emission from the NGC 7538 region is typically in the velocity range with respect to the local standard of rest (v_LSR) between −60 and −50 km s⁻¹, these positions show additional strong emission at even more negative velocities beyond −60 km s⁻¹. In particular, the spectra toward BlPN, RPC, BrPN, and BRPS exhibit a second, clearly separated velocity component in the [CII] emission that has barely any detectable counterpart in the molecular emission of the CO(3–2) line or the atomic [CI] emission. The rms values for these two lines (Sect. 2) correspond at 30 K to 3σ column density sensitivities (see Sect. 4.3 for details on the calculations) of N(CO) ≈ 2.1 × 10¹⁶ cm⁻² and N(C^0) ≈ 2.7 × 10¹⁶ cm⁻² per channel (0.42 and 0.5 km s⁻¹ for CO(3–2) and [CI], respectively).

This second component can also be identified toward IRS5 and BlPS. Hence, the ring-like high velocity dispersion structure seen in the second-moment map of [CII] (Fig. 4, right panel) is in fact not a region of particularly high velocity dispersion, but there we have two velocity components mimicking high values in the second-moment map. The peculiar aspect of this additional velocity component is that it is mainly detected in the ionized carbon [CII] line without a strong molecular or atomic counterpart being detectable in the individual spectra. However, Sandell et al. (2020) detected a north-south outflow in CO(3–2) and [CII] emission (their [CII] map is centered on IRS1 and smaller than the one presented here) emanating from IRS1, where the blue-shifted emission is located north of IRS1. The positions BlPN, BlPS, RPC, BrPN, and BrPS are all in the general vicinity of that outflow. In particular, the morphology of the blue-shifted [CII] and CO(3–2) gas emission near IRS1 (corresponding roughly to our position BlPS) is similar. Hence, the blue-shifted gas we find may at least partly be associated with that large-scale outflow. Nevertheless, in the spectra extracted at individual positions as shown in Fig. 5, the blue-shifted gas is much stronger in the [CII] spectra than in the more commonly as outflow tracer observed CO emission. We discuss this component further in Sect. 4.4.

To look at the velocity structures in a different fashion, Fig. 6 presents a channel map of the [CII] emission in NGC 7538 (a version without any suggested bubble structures is presented in Fig. A.2). While in this representation the main emission of the PDR associated with the H II regions is also dominant approximately between −60 and −50 km s⁻¹, there are considerable extended emission structures at blue-shifted velocities < − 60 km s⁻¹, and to a lesser degree, also at red-shifted velocities > − 50 km s⁻¹. For comparison, Fig. A.3 presents the corresponding CO(3–2) channel map, where these more blue- and more red-shifted features are barely recognizable. The [CII] channel map gives the impression of several bubble- and ring-like structures that may stem even from several star formation events. We will discuss the bubble-features and their possible interpretation in more detail below in Sect. 4.1.

In addition to the molecular CO, atomic [CI], and ionized carbon [CII] emission, we can also investigate the atomic hydrogen by means of the 21 cm line observed with the Canadian Galactic Plane Survey (CGPS, Taylor et al. 2003). As expected, toward the main H II region around the exciting sources IRS5 and IRS6, the HI is seen only in absorption. We can clearly identify a strong HI absorption component associated with the high-velocity blue-shifted gas < − 60 km s⁻¹, but the main velocity component of the region around − 55 km s⁻¹ exhibits only weak absorption. The blue-shifted HI absorption should be related to the expanding shell accelerating also the atomic gas envelope in the direction of the observer. The red-shifted absorption at ~ − 45 km s⁻¹ cannot be the other side of the expanding shell because in absorption spectroscopy, it has to lie in front of the ionized gas. Hence, the red-shifted HI absorption must be related to some foreground gas. Since there is barely any [CII] nor CO emission at these velocities, that component may even be unrelated to the NGC 7538 H II region.

Looking a bit outside the actual H II region toward the east (position RPE, middle panel of Fig. 7), the HI emission spectrum is broader than what we observe in [CII] and CO. This can be understood in a way that HI is more easily excited than [CII] and CO, hence, it can pick up more tenuous gas at higher velocities. Furthermore, at the velocity of the [CII] and CO(3–2) peak emission, the HI shows a dip in the spectrum which may be caused by HI self-absorption (HISA, e.g., Gibson et al. 2005; Syed et al. 2020).

Furthermore, if we look toward positions where we can clearly identify two components in the [CII] emission, e.g, the position BlPN in the north of the H II region (bottom panel of Fig. 7), it may at first sight be surprising that the strong blue-shifted [CII] component at velocities < − 60 km s⁻¹ is inconspicuously weak in the HI emission. One potential explanation for that could be that the radiation field may be that strong that most of the gas is ionized hydrogen and carbon with just a thin layer of neutral HI and molecular gas in the surrounding molecular cloud.

Fig. 3 Overlays of selected pairs of emission tracers as described above each panel. The integration ranges for the atomic and molecular lines are the same as in Fig. 2. The contour levels are chosen individually for each figure to highlight the emission best. The angular resolution elements for both tracers are presented in the top-right of each panel. The five-pointed stars always mark the position of the embedded source NGC 7538IRS1 and the two open stars show the positions of the two main exciting sources of the H II region NGC 7538IRS5 and NGC 7538IRS6. A linear scale bar is presented in the bottom-left panel.

Fig. 4 Moment maps of the [CII] emission. Left, middle, and right panels: present (in color scale) the integrated emission and the first- and second-moment maps (intensity-weighted peak velocities and intensity-weighted velocity dispersion) in the velocity range from −75 to −40 km s−1, respectively.The noisier top-right quarters of the first- and second-moment maps are blanked (see Sect. 2). The contours in the middle and right panels show the integrated [CII] emission from the left panel at only a few selected contour levels to highlight similarities and differences. The three five-pointed stars in the left panel mark the positions of the embedded sources NGC 7538IRS1, NGC 7538S, and NGC 7538IRS9. The two stars show the positions of the two main exciting sources of the H II region NGC 7538IRS5 and NGC 7538IRS6. A linear scale bar is presented in the left panel. The markers in the middle and right panels highlight additional positions (see main text) around the areas where we extracted spectra in Fig. 5. The line in the middle panel outlines the pv-cut from IRS6 through the bar-like structure presented in Fig. 13.

Fig. 5 [CII] (black), CO(3–2) (red), and [CI] (teal, multiplied by 2) spectra toward the positions marked in Fig. 4. The green lines show Gaussian fits to the [CII] and CO(3–2) data. The numbers present the corresponding FWHM values (“*/*” labels two components). The positions are labeled in each panel. The RPE position to the east is not covered by the [CI] map.

Fig. 6 [CII] channel map for the velocities labeled in each panel. The data are binned for that plot in 2 km s−1 channels. The two purple stars mark the positions of IRS1 and IRS5, and a scale bar is shown in the top-left panel. The intensities are always plotted logarithmically from 1 to 40 K. The purple circles outline potential bubbles. These bubbles are labeled in the top-left panel. A version without any suggested bubble structures is presented in Fig. A.2.

Fig. 7 [CII] (black), CO(3–2) (red), and HI (blue, divided by 10) spectra toward three selected positions marked in Fig. 4. These are IRS5, the red-shifted peak toward the east (RPE), and the blue-shifted peak in the north (BlPN). The positions are labeled in each panel.

4 Discussion

4.1 Identifying expanding bubbles or an inhomogeneous medium

Bubble/ring identification

As indicated in Sect. 3.2, the channel map of the [CII] emission shows several ring-like structures by eye-inspection (four rings are outlined in Fig. 6 as magenta circles). These features can be limb-brightened edges of expanding shells or bubbles, but they could also be a ring or torus in the plane of the sky. Before we discuss their velocity structure in more detail, we apply an unbiased automized ring identification by computing the covariance between the map and annuli structures of varying radii. The details of the identification analysis are presented in Appendix B. We stress that the bubbles (or rings) discussed below are only candidates for real physical structures that may stem from the expansion of the H II region. The bubble- or ring-like structures we identify this way have central positions (in J2000.0) and approximate radii, r, of: (A) RA = 23h13m38s, Dec = 61°32′20″; r ~ 100″ ~1.3 pc; (B) RA = 23h13m38s, Dec = 61°30′12″; r ~ 100″ ~1.3 pc; (C) RA = 23h13m46s, Dec = 61°31′00″; r ~ 100″ ~1.3 pc; (D) RA = 23h13m55s, Dec = 61°31′56″; r ~ 80″ ~1.0 pc. As outlined in Appendix B, the approach identifies even a fifth ring-like structure (labeled E), however, that appears rather as an overlap of mainly structures A to C, and we do not consider that further as a separate physical entity. Uncertainties for the central positions and radii are also discussed in Appendix B. For comparison, the ring and bubble candidates are also plotted on the 8 μm emission that mainly stems from PAHs emitted in PDRs (Fig. 8). Several of the ring- and bubble-like structures can clearly be seen in this PAH emission.

The suggested bubble labeled A is best visible as limb-brightened rings in the channels between −65 and −59 km s⁻¹. An additional indicator that this bubble may indeed be a real physical structure is that at the most blue-shifted velocities (channel at −69 km s⁻¹), the emissionis mainly at the center of this proposed bubble A. This is exactly what we would expect when this bubble is expanding along the line of sight (see also Appendix B) and has also been shown for RCW120 (Luisi et al. 2021).

Fig. 8 Potential bubble structures and 8 μm emission. The color scale shows the 8 μm data and the purple circles again outline the tentative bubble structures.

Position-velocity analysis

A different way to visualize the velocity structure of bubble A is via a position-velocity cut. Figure 9 presents three vertical pv-diagrams in the [CII] line through the center of bubble A and shifted 1′ to the east and west, respectively. Corresponding horizontal pv-diagrams are presented in Fig. A.4. While most velocities along these cuts are between − 50 and − 60 km s⁻¹, there is clearly high-velocity blue-shifted gas around ~ − 70 km s⁻¹ close to the center of the bubble, indicating an expanding structure along the line of sight. Figures 9 and A.4 also show the shape a spherical shell with a radius if 100″ and an expansion velocity of 14 km s⁻¹ would have in such a pv-diagram. This is not a fit but just outlining how such a shell would appear as an ellipse in the pv-diagram. The left and right panels with pv-cuts at offsets of 1′ correspond in a shell geometry to angles of ~36 deg (see a sketch of such a geometry in Appendix A.2 from Butterfield et al. 2018). The velocities at such angle are then reduced by cos (36°). While the data and toy model at negative blue-shifted velocities show some resemblance of each other, consistent with a spherical half-shell, the red-shifted side at positive relative velocities shows no indication of such a shell-like geometry. This could either indicate that such a shell model is insufficient or the center of the shell could be so close to the rear side of the cloud that only blue-shifted gas moving toward the observer can be detected, while the red-shifted expansion leaves the cloud almost immediately that no gas gets accelerated to a detectable level. In that context, Luisi et al. (2016) estimated that approximately 15% of the ionized gas is leaking outside of the H II region.

Structure B is centered close to the main excitation sources in the region, IRS5 and IRS6 (Fig. 1). However, we note that the automated bubble-identification approach does not center it exactly on the excitation sources but it is a bit shifted. Bubble B’s border at the southwestern side roughly replicates the bar-like PDR discussed in Sect. 4.2. In the north, structures visible in particular at velocities of − 55/ − 53 km s⁻¹ delineate the potential sphere of that bubble. Carrying out position-velocity cuts in the horizontal as well as vertical direction through this bubble B and source IRS5 (Fig. 10), we find that the highest-velocity gas on the blue- as well as red-shifted side is found close to IRS5. This is again indicative of expanding motions along the line of side, potentially caused by IRS5 and IRS6. We again draw the corresponding shapes for expanding shells into the pv-diagram (Fig. 10). While the radius of 100″ is given by our shell-identification approach above, the velocities on the blue- and red-shifted side differ significantly. To outline those differences, we used an expansion velocity of 10 km s⁻¹ for the blue-shifted side, and 3 km s⁻¹ for the red-shifted side. If the observed structure really belongs to an expanding shell, such velocity differences would indicate a highly inhomogeneous medium where the red-shifted expansion would enter into much denser material. Another option is that the central velocity is not − 55 km s⁻¹ but a few km s⁻¹ higher. In that case, we would just have a blue-shifted bubble-structure without any red-shifted counterpart, similar to structure A, discussed before. That would put the center of structure B also on the rear side of the cloud.

The identified structure C is in-between the other bubble or ring structures. While it may be a real physical entity, it could also arise from an overlap of the other bubble or ring structures. Therefore, we refrain from a further discussion of that structure.

Furthermore, we identify structure D, where the limb-brightened ring-like features are best delineated in the channels between − 57 and − 53 km s⁻¹. In addition to the bubble edges, structure D again shows high-velocity gas at its center, typically signposting expanding structures. These centered emission structures in bubble D are best visible in the blue-shifted channels at − 67/ − 65 km s⁻¹ and in the most red-shifted channels at − 49/ − 47 km s⁻¹. For comparison, Fig. 11 presents the corresponding vertical and horizontal pv-cuts through the central position of structure D.In contrast to the pv-diagrams of bubbles A and B (Figs. 9 and 10) that show the high-velocity gas mainly at negative blue-shifted velocities, the pv-diagram of structure D, in particular the one in vertical direction (right panel in Fig. 11), shows an almost ring-like structure with high velocities on the blue- and red-shifted side, in good agreement with the channel-map discussed before. For comparison, we again show the shape an expanding shell with radius 80″ and an expansion velocity of 6 km s⁻¹.

Fig. 9 Vertical position-velocity cuts through bubble A from south to north derived from the [CII] line emission. The central panel goes through the estimated central position of A whereas the left and right panels are shifted 1′ to the east and west, respectively. The white dashed line marks the approximate center of the bubble. The black dotted lines in the middle panel outline the elliptical shape that a spherical shell with a radius of 100″ and expansion velocity of 14 km s−1 would have in the pv diagram. The dotted lines in the left and right panel correspond to the same expanding shell at the 1′ shift to the east and west where the 1′ corresponds to ~36 deg offset in the sphere (e.g., Butterfield et al. 2018, sketch in their Appendix A.2); hence, the velocities are then reduced by cos (36°).

Driving sources

It is interesting to note that the ring-like structure of high velocity dispersion in the [CII] moment map (Fig. 4) roughly coincides with the combined outline of the four bubble or ring structures discussed here (Fig. 6). Since the high velocity dispersion in Fig. 4 largely stems from several velocity components (Fig. 5), these velocity signatures may stem from expanding motions from the H II region. While for bubble B, the driving source for bubble expansion may be the main exciting sources of the H II region (IRS5 and IRS6, Puga et al. 2010), this is less obvious for the other structures A, C, and D. For these three structures, we have not identified a clear central source that could be associated with bubble expansion. While Puga et al. (2010) identified more young stellar objects in the environment of the H II region, they are not at the bubble centers (and, in fact, IRS5 and IRS6 are not exactly in the center either), and, in addition, they are typically of lower luminosity. Also, the Spitzer data do not allow us to identify clear driving sources.

While the O3 and O9 stars IRS6 and IRS5 (Puga et al. 2010) are most likely the main exciting sources of the H II region, the identification of several more ring- or bubble-like structures does not necessarily mean that there are more bubble-driving sources in the field. While proper motion of the exciting stars may cause the dislocation of IRS5 and IRS6 from the center of bubble-structure B, such displacement may also be caused by density and pressure gradients in the region. A similar displacement is also observed between the Trapezium stars and the main bubble in Orion, namely, the Orion Veil (Pabst et al. 2019, 2020).

Fig. 10 Horizontal (along RA, left) and vertical (along Dec, right) position-velocity diagrams through the position of one of the exciting sources IRS5 (marked in each panel) derived from the [CII] line emission. corresponding roughly to the estimated center of bubble B. The contour levels are in 3σ steps. The black dotted lines outline the elliptical shape a spherical shell would have in the pv diagram – given a radius of 100″ and expansion velocities of 10 and 3 km s−1 at negative and positive relative velocities, respectively.

Fig. 11 Horizontal (left) and vertical (right) position-velocity diagrams through the center of bubble D (marked in each panel) derived from the [CII] line emission. The contour levels are in 3σ steps. The black dotted lines outline the elliptical shape a spherical shell would have in the pv diagram – given a radius of 80″ and expansion velocities of 6 km s−1.

Multiple bubbles versus expansion into an inhomogeneous medium

In addition, the strong sub-structure of the H II region visible in the different data sets, in particular in the [CII] channel map, indicates an almost porous structure of the H II region where radiation and wind energy can leak out and also influence parts of the H II region that are less close to the main exciting sources IRS5 and IRS6. Hence, while bubble-like structures are possible, the NGC 7538 H II region may also expand into a very inhomogeneous environment. If an H II region expands in such an inhomogeneous cloud, changes in the cloud morphology are expected. Such impact may also cause bubble-like morphologies as observed here. Recent analytic and numerical work by Lancaster et al. (2021a,b,c) also showed that the wind-driven expansion of the gas around H II regions can cause different morphological sub-structures. For the NGC 7538 H II region, a clear discrimination between several bubble-structures each driven by individual internal sources or an expansion of mainly one wind-driven shell into an inhomogeneous medium cannot be drawn. Further hydrodynamic modeling of expanding H II regions into an inhomogeneous medium may show whether such multiple ring-like structures of expanding gas can occur or whether separate bubble-driving centers are required.

Thermal expansion or winds of the high-mass stars

Independent of whether all the ring-like structures are real bubbles or whether they are rather caused by an expansion of the gas into an inhomogeneous medium, the high-velocity components visible in the pv-diagrams(Figs. 9–11) all show high-velocity gas likely driven by the expansion of the H II region. In particular the gas toward structures A and B exhibit expansion velocities ≥ 10 km s⁻¹. Taking an expansion velocity of 10 km s⁻¹ and a radius of 100″ at face-value, this would correspond to an expansion time-scale of roughly 0.126 Myrs. This is comparably short with respect to the observationally estimated ages of the cluster. For example, Puga et al. (2010) derived an age range between 0.5 and 2.2 Myrs, and Sharma et al. (2017) estimated a mean age of the young stellar objects of 1.4 Myrs with a range between 0.1 and 2.5 Myrs. Hence, a constant expansion at that velocity is unlikely the real scenario for that region. However, we consider whether the measured expansion velocities might be caused by thermal expansion of the gas or by the winds of the high-mass stars.

To estimate expected expansion velocities, for the thermal expansion, we followed the Spitzer solution (Spitzer 1998), whereas for the wind-driven expansion, the approach by Weaver et al. (1977) and its adaption to include radiative cooling by Mac Low & McCray (1988) was adopted. A comprehensive summary of these two approaches is given in Henshaw et al. (2022). The expansion velocities depending on time for the thermal expansion v_th (t) and the wind-driven expansion v_w(t) are:

Here, c_s is the sound speed in the ionized gas (8 km s⁻¹ at 5000 K), is the mechanical wind luminosity, ρ the ambient density, and t_cool the cooling time following Mac Low & McCray (1988) and Henshaw et al. (2022). The wind mass flow rate Ṁ and the escape velocity v_∞ for an O3 star (IRS6, Puga et al. 2010) are taken from Muijres et al. (2012). The wind-driven velocities, v_w, strongly depend on the density of the ambient gas ρ (also t_cool depends on ρ, Mac Low & McCray 1988).

Following Puga et al. (2010) and Sharma et al. (2017), the age of the NGC 7538 cluster and associated young stellar objects is between roughly 0.1 and 2.5 Myr. Even at the lower boundary of 0.1 Myr, the highest expansion velocities that we can get with a Spitzer-type thermal expansion is ~5.7 km s⁻¹, far below what is measured in NGC 7538. So, purely thermal expansion cannot properly explain the measured high-velocity gas. In contrast, the wind solution gives significant higher velocities. For the given parameters, and assuming densities of 10³, 2 × 10³ and 3 × 10³ cm⁻³, expansion velocities around 10 km s⁻¹ are estimated at times of roughly 1.1, 0.36 and 0.17 Myr, respectively. While the lower end would be more consistent with the constant expansion estimated above, the higher timescales are more consistent with the estimated age limits by Puga et al. (2010) and the mean ages of the young stellar objects derived by Sharma et al. (2017). Hence, wind-driving seems a plausible way to explain the observed high-velocity gas. The wind-driving is also consistent with the diffuse X-ray emission toward that region that is typically attributed to hot plasma caused by the wind shocks of the massive stars (e.g., Güdel et al. 2008; Townsley et al. 2018; Pabst et al. 2020).

Regarding the potential discrepancy between the lower end of the estimated cluster age at 0.5 Myrs and the shorter timescales estimated for either constant expansion or the wind expansion at higher densities, one way to reconcile these is that initially the bubble may have been confined inside the dense molecular core in which IRS5 and IRS6 have formed. It is only once the H II region broke out of this core, expansion into the surrounding lower density material gained speed. In this picture, the expansion timescale estimated here would then refer to the time since the breakout from this dense core.

4.2 Bar-like photon dominated region (PDR)

Figure 12 presents a zoom into the bar-like structure at the south-western edge of the expanding H II region. While the color scale shows the ionized gas as traced by the 1.4 GHz continuum emission, the contours outline various other tracers: 8 μm continuum, [CII], CO(8–7) and (3–2), and [CI] emission. The white, black, blue, and red lines mark the approximate locations of the emission crests in the 8 μm, [CII], CO(3–2) and [CI] emission, respectively. The crests were identified by eye via connecting the emission peaks. The crest of the CO(8–7) emission approximately coincides with that of the [CII] emission.

In this picture, the 8 μm emission, that should stem largely from UV-pumped infrared fluorescence PAHs, peaks closest to the exciting sources of the H II region. Nextin this layered structure, we have the ionized carbon [CII] emission that spatially approximately coincides with the highly excited CO(8–7) emission. Furthest away from the exciting sources, we have the emission crests of the lower excited CO(3–2) and [CI] emission that, given the spatial resolution of the data and when projected on the plane of the sky, approximately coincide spatially. This is different around the young high-mass star-forming region IRS1 south-east of the bar where CO(3–2) peaks closer to IRS1 than the atomic carbon [CI] line (Fig. 12). Coming back to the bar-like structure, in this framework, we can separate the photon dominated region (PDR) into roughly four layers:

• (i)

  The exciting sources and 1.4 GHz continuum emission;

• (ii)

  the 8 μm emission, likely PAHs;

• (iii)

  ionized carbon [CII] and highly excited CO(8–7);

• (iv)

  atomic carbon [CI] and lower excited CO(3–2).

This layeredstructure resembles prototypical edge-on PDR structures observed also for example in the Orion bar (e.g., Tielens et al. 1993; Tauber et al. 1994, 1995; Simon et al. 1997) or M17SW (e.g., Pérez-Beaupuits et al. 2015b,a). We return to the carbon budget within this PDR in Sect. 4.3.

In addition to the layered structure, the first-moment map in Fig. 4 indicates that there may be a velocity gradient across the bar. To get a better view of that, Fig. 13 presents a position-velocity diagram across the bar starting from the O3 star IRS6. We see a clear velocity gradient from approximately − 50 km s⁻¹ at the beginning of the cut near the main exciting source IRS6 to blue-shifted velocities < − 60 km s⁻¹ on the otherside of the bar. The fact that the gas within the bar gets blue-shifted indicates that the bar lies along the line of sight to the observer in front of the exciting source IRS6.

Fig. 12 Zoom on the bar-region south of the exciting sources IRS5 and IRS6. The color scale and the colored contours show the emission tracers labeled in the two panels. Only selected contours are shown for clarity. The four colored lines mark the approximate positions of the emission crests in the corresponding tracers. The sources IRS5, IRS6, and IRS1 are marked. A linear scale bar is shown as well.

Fig. 13 Position-velocity cut from the [CII] line from IRS6 through the bar as outlined in the middle panel of Fig. 4. The white dashed line marks the approximate center of the bar.

4.3 Carbon budget

While carbon is mainly in its ionized C⁺ state in large fractions of the PDR, it is mainly in its molecular CO form in the dense molecular cloud (Fig. 2). But we ought to consider the relative carbon budget within the transition zone from the H II region toward the molecular clouds. This transition is best studied in the southwestern bar-like structure. Figure 14 presents the ionized carbon [CII], the atomic carbon [CI] and the molecular carbon CO(3–2) emission toward that region. The integration regime for all three species is the same from − 70 to − 40 km s⁻¹. In the following, we derive approximate estimates for the column densities and masses of the different carbon components in the bar region outlined by Fig. 14. With the assumptions outlined below, these should only be considered as rough estimates; nevertheless, they allow us to assess roughly how much mass is within each carbon component in and around the bar-like region.

The column densities of the ionized carbon are estimated from the [CII] line following Goldsmith et al. (2012) (Eq. (26)), adding an additional optical depth term of to correct for the optical depth τ (Rohlfs & Wilson 2006):

with the kinetic temperature T_kin (see further below), the collision rate C_ul, the intensity T_mb in K, and the velocity dv in km s⁻¹. The collision rate C_ul = R_uln depends on the temperature, where R_ul is collision rate coefficient with H₂ and n the density, assumed to be around ~10³ cm⁻³. The collision rate coefficients R_ul with H₂ are taken from the Leiden database for molecular spectroscopy (Schöier et al. 2005⁶ ; Lique et al. 2013) where the para- and ortho-rates are weighted following Le Bourlot (1991) and Gerlich (1990). The comparably good spatial correspondence of the [CII] with the molecular line emission around the bar-shaped PDR indicates that H₂ should be the most dominant collisional partner. If, for comparison, all collisions involved atomic hydrogen, with the given assumptions, we would get masses roughly a factor of 1.57 lower (see also Goldsmith et al. 2012). In reality, it is likely a mix of collisional partners dominated by H₂ and hence the difference should be much smaller. In the following we use the collisional rates with H₂. Evaluating the [¹³CII] emission in our data, we detect the line after averaging over areas to decrease the noise. In the case of the bar-like region in NGC 7538, the [¹³CII] data indicate a mean optical depth τ for the main [CII] line of ~2.5. Compared to the optically thin case, this corresponds to a column density and mass correction factor of ~2.7.

Furthermore, the atomic carbon column densities are estimated from the [CI] emission following Frerking et al. (1989):

with theexcitation temperature, T_ex (the temperature quantification is discussed further below). The CO column density can then approximately be estimated following standard equations (e.g., Cabrit et al. 1988):

Here, ν, μ, T_ex, E_u ∕k, τ and T_mb are the frequency, dipole moment (0.112 Debye), the excitation temperature, the upper energy level of 33.19 K, the optical depth, and the main beam brightness temperature of the line, respectively.

Following Cabrit et al. (1988), assuming ¹³CO is optically thin, the optical depth term can be approximated as:

The term corresponds to the ¹²CO/¹³CO abundance ratio. According to Moscadelli et al. (2009), NGC 7538 is in the Perseus arm which results at the given longitude of the region (111.5 deg) in a Galactocentric distance of ~9 kpc (e.g., Reid et al. 2019). Following Wilson & Rood (1994), the ¹²CO/¹³CO abundance ratio is then ~75. Furthermore, the term correspondsto the mean ratio of brightness intensities between the ¹³CO and the ¹²CO transition. Within the IRAM NOEMA large program CORE (Beuther et al. 2018), a small area of 1.5′ × 1.5′ around IRS1 had been mapped in the J = 2−1 transitions of ¹²CO and ¹³CO with the IRAM 30 m telescope. Using these small-area data, we estimate a mean ratio between those two lines of ~ 0.2 in the velocity regime between − 70 and − 40 km s⁻¹. This value is used in the following also for the here observed J = 3−2 transition. This results in the following column density equation:

In the following, we will calculate the mean column densities and the corresponding masses overthe area presented in Fig. 14 that encompassed roughly 11.25 pc². The main missing ingredient in the above equations are the temperatures to be chosen for the column density estimates.

Regarding the ionized carbon C⁺, Langer et al. (2010) typically assume temperatures between 100 and 150 K, which seems a reasonable temperature regime in this PDR as well. Using these two values as brackets for the possible temperatures of the ionized carbon C⁺, we find an ionized carbon C⁺ mass of M_[CII] ≈ 9.7 ± 1.6 M_⊙, assuming the above discussed mean optical depth τ ~ 2.5. This is an upper limit to the mass (Rohlfs & Wilson 2006), whereas we can consider an estimate assuming optically thin emission as a lower limit. In the optically thin limit, we get the lower limit of M_{[CII]_lower_limit} ≈ 3.6 ± 0.6 M_⊙, a factor of ~2.7 below the optical depth corrected mass estimate.

For the atomic and molecular gas, the excitation temperatures are not well determined in that region. However, variations in the excitation temperatures between 30 and 50 K do not vary the estimated masses significantly. Therefore, we estimate the atomic M_[CI] and molecular masses M_CO in that excitation temperature regime. In that framework, we find M_[CI] ≈ 0.45 ± 0.1 M_⊙ and M_CO ≈ 1.2 ± 0.1 M_⊙. The molecular mass can also be compared to the gas mass, which we can derive from the Herschel gas column density map. While a H₂ -to-CO ratio of 10⁴ would result in a gas mass of ~ 1.2 × 10⁴ M_⊙, the gas mass derived from the Herschel map in the same area is ~4600 M_⊙. This difference of around a factor of 2.6 can already be explained by the uncertainties of temperatures and line-of-sight averaging, optical depth estimates, and the assumed dust properties for the Herschel-based mass estimate. Nevertheless, this comparison shows that our molecular mass estimate has an approximate uncertainty of a factor of 2.

While atomic C⁰ and molecular CO carbon have comparably low masses between 0.45 and 1.2 M_⊙, the mass of the ionized carbon C⁺ is in the optically thin lower limit a factor of ~3–8 larger, and for the optical depth corrected mass estimate even a factor of ~8–22 larger in this bar-like PDR (Fig. 14). Hence, the ionized carbon C⁺ is clearly the dominant carbon form in this PDR. As expected, this is very different to the previous carbon studies in infrared dark clouds, where the carbon was found to be dominantly in the molecular CO form (Beuther et al. 2014).

Fig. 14 Integrated [CII] (color), [CI] (black contours), and CO(3–2) (green contours) emission in the bar region where the [CII]/[CI]/CO mass estimates are conducted (Sect. 4.3). All three maps were integrated between −70 and −40 km s−1. The [CI] contours start at 6σ and continuein 3σ steps (1σ = 1.6 K km s−1). The CO(3–2) contours start at 10σ and continuein 10σ steps (1σ = 3.2 K km s−1). The white line marks the approximate location of the bar-like feature.

4.4 Decoupled [CII] component

The [CII] component at velocities below − 60 km s⁻¹ deserves a separate discussion because it is identified mainly in the [CII] emission without strong counterparts in the molecular and atomic gas (Figs. 5 and 6). A spectral feature dominantly visible in [CII] emission has rarely been identified in the past in regions of expanding H II regions and feedback. For comparison, the veil in Orion exhibits barely any extended CO emission, however, compact globules are identified within the [CII] emitting gas structures (Goicoechea et al. 2020).

Inspecting the channel map with proposed bubble-like structures (Fig. 6) in velocity channels at ≤ − 65 km s⁻¹, this blue-shifted high-velocity gas appears to be associated with all four identified bubble-structures as well as with the young star-formingregion IRS1 in the south of the region. Since IRS1 is known to be an active driver of molecular outflows (e.g., Scoville et al. 1986; Kraus et al. 2006; Qiu et al. 2011; Beuther et al. 2012; Sandell et al. 2020), the discovery of high-velocity emission in its vicinity does not come as a surprise. Furthermore, the blue-shifted outflow lobeemanating from IRS1 toward the north is aligned right with the central area of the H II region and may explainpart of the blue-shifted emission (Sandell et al. 2020). Nevertheless, the [CII] emission below − 60 km s⁻¹ is significantly stronger than the CO(3–2) emission at the same velocities (Fig. 5). Hence, a sole explanation by the outflow appears unlikely, and a significant fraction of the highly blue-shifted gas appears to be associated with the suggested bubble-structures.

Using the same assumptions described in the previous section, we can estimate the mass of ionized carbon associated with this blue-shifted [CII] component. Integrating the blue-shifted [CII] line between − 75 and − 65 km s⁻¹ over the entire emission area (~22.5 pc²), we find in the temperature regime between 100 and 150 K a mass of M_⊙. Although the area is roughly a factor of 2 larger than that used for the mass estimates for the bar in the previous section, the mass in this blue-shifted component is more than a factor of 10 smaller than that of the bar-like structure. Nevertheless, these data show that a significant fraction of ionized carbon is still found in that mainly in the [CII] line emitting velocity component.

Assuming that at least part of the blue-shifted high-velocity gas associated with the proposed bubble-like structures is indeed associated with expanding gas along the line of sight, the fact that we see barely any atomic or molecular carbon counterpart implies that there is very low (or even no) interaction of the expanding gas with a denser environmental cloud.

4.5 Shell structures associated with HII regions and PDRs in [CII] emission

As mentioned in the introduction, the existence of expanding shell-like structures around H II regions and PDRs is a recurring result of the SOFIA FEEDBACK program, where NGC 7538 is part of, as well as other [CII] studies. The best studied examples are Orion (Pabst et al. 2019), RCW120 (Luisi et al. 2021), and RCW49 (Tiwari et al. 2021). From a geometry point of view, the Orion Veil and RCW120 exhibit more complete shell-like structures in the [CII] emission. The NGC 7538 data presented here also show bubble-like features, but with much more sub-structure and no real large-scale shell. RCW49 again exhibits a large-scale shell-structure, however, depending on the velocities of the gas, the shell can be broken open, allowing radiation to escape into the environment.

Regarding the expansion speeds measured in the [CII] emission, all four regions discussed here have similar expansion velocities between roughly 10 and 16 km s⁻¹. Assuming constant expansion, the given velocities result in estimated expansion timescales between 0.1 and 0.27 Myr, all within a similar range. The expansion of the four regions can in neither case be explained by thermal expansion like in the classical Spitzer picture (Spitzer 1998), but wind-driving from the central O-stars is needed.

Furthermore, it is interesting to note that in three of the regions – RCW120, RCW49, and NGC 7538 – star formation at the edges of the expanding regions is observed. For NGC 7538, triggering as a potential scenario is discussed in detail by Sharma et al. (2017). Hence, triggered star formation seems a regular outcome in these kinds of regions (see also general bubble and protostar correlation analysis by, e.g., Kendrew et al. 2016 or Palmeirim et al. 2017). In contrast, no obvious star formation is identified at the edges of the Orion Veil. Goicoechea et al. (2020) recently found small globules of molecular gas within the Veil structure, however, these are rather transient objects and unlikely to form stars. Hence, while the star formation activity around RCW120, RCW49, and NGC 7538 hint toward at least partially positive feedback triggering new star formation events – in the Orion Veil, the feedback processes expel and reprocess the gas, limiting further star formation.

5 Conclusions and summary

Using SOFIA, we were able to map roughly 210 square arcmin (125 pc²) around the H II region NGC 7538 in the fine-structure line of ionized carbon [CII] as well as fine-structure [CI] and CO(8–7) emission. These SOFIA data are complemented with archival near-/far-infrared/cm continuum and CO(3–2) line observations.

The [CII] emission reveals a wealth of information about the structure and composition of this region of high-mass star formation. While the overall [CII] morphology may at first sight appear similar to that of the ionized gas traced by the cm continuum emission, the velocity-resolved [CII] emission reveals a considerable substructure resembling the PDR that is otherwise well traced in the 8 and 70 μm emission. Using a automated bubble or ring-identification algorithm, we identify four bubble-like structures. Two of them show kinematic signatures of expanding bubbles. One of these structures is centered near the main exciting sources IRS6 and IRS5. For the other three bubble-like structures, no clear central source can be identified. This may indicate that at least some of the bubble-like structures could be mimicked by the expansion of the gas into an intrinsically porous and inhomogeneous medium within the H II region.

An analysis of the expansion velocities based on high-velocity [CII] gas components reveals that purely thermal expansion from a Spitzer-like H II region is not sufficient, but that wind-driving by the central O-stars is needed. This is supported by Chandra observations that reveal the presence of diffuse X-ray emission due to hot plasma, resulting from the shocked stellar winds from these O stars.

The most blue-shifted [CII] velocity component barely shows any counterpart in atomic or molecular carbon emission. At the interface of the H II region and the molecular cloud, we find a bar-like photon-dominated region where ionized C⁺, atomic C⁰, and molecular carbon CO reveal a layered structured. The carbon in the PDR is dominantly in its ionized C⁺ form. A velocity-gradient across that bar-like PDR is identified, indicating that the PDR is pushed along the line of sight in the direction of the observer.

The [CII] study of this prototypical H II region NGC 7538 reveals a highly sub-structured H II region that interacts with the adjacent molecular cloud, potentially triggering new star formation events in some of the adjacent dense cores. Setting the region into context with other recently investigated [CII] emission studies from the FEEDBACK program (Orion Veil, RCW120, RCW49), we find that bubble-like expanding morphologies are recurring structures in PDRs. While in some of these regions, star formation occurs at the edges of the expanding regions, indicating positive feedback – the data from the Orion Veil reveal the destructive forces of feedback, expelling and reshaping the gas in that region. Hence, both positive and negative feedback can occur as a result of the expansion of H II regions. The SOFIA legacy program FEEDBACK will dissect more H II regions and reveal common processes in H II region-molecular cloud interaction zones.

Acknowledgements

Financial support for the SOFIA Legacy Program, FEEDBACK, at the University of Maryland was provided by NASA through award SOF070077 issued by USRA. The FEEDBACK project is supported by the BMWI via DLR, Projekt Number 50 OR 1916 (FEEDBACK) and Projekt Number 50 OR 1714 (MOBS - MOdellierung von Beobachtungsdaten SOFIA). The Digitized Sky Surveys were produced at the Space Telescope Science Institute under US Government grant NAG W-2166. The images of these surveys are based on photographic data obtained using the Oschin Schmidt Telescope on Palomar Mountain and the UK Schmidt Telescope. The plates were processed into the present compressed digital form with the permission of these institutions. This work is based in part on observations made with the Spitzer Space Telescope, which was operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. We like to thank the referee for a constructive report improving the paper. H.B. acknowledges support from the European Research Council under the Horizon 2020 Framework Program via the ERC Consolidator Grant CSF-648505. H.B. also acknowledges support from the Deutsche Forschungsgemeinschaft in the Collaborative Research Center SFB 881 - Project-ID 138713538 - “The Milky Way System” (subproject B1). N.S. acknowledges support by the Agence National de Recherche (ANR/France) and the Deutsche Forschungsgemeinschaft (DFG/Germany) through the project “GENESIS” (ANR-16-CE92-0035-01/DFG1591/2-1). This work was supported by the German Deutsche Forschungsgemeinschaft, DFG project number SFB 956. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 851435).

Appendix A Complementary figures

Fig. A.1 [CII] (black) and CO(3–2) spectra toward the grid marked in the top-left [CII] first-moment map. The contours show the [CII] emission at only a few selected contour levels to highlight the main structures. The sources IRS5, IRS6, and IRS1 are also shown as black-and-white stars.

Fig. A.2 [CII] channel map for the velocities labeled in each panel. The data are binned for that plot in 2 km s−1 channels. The two purple stars mark the positions of IRS1 and IRS5, and a scale bar is shown in the top-left panel. The intensities are always plotted logarithmically from 1 to 40 K.

Fig. A.3 CO(3–2) channel map for the velocities labeled in each panel. The data are plotted in steps of 2.1 km s−1 (5 channels). The two purple stars mark the positions of IRS1 and IRS5, and a scale bar is shown in the top-left panel. The intensities are always plotted linearly from 0.1 to 10 K. The purple circles outline potential bubbles. These bubbles are labeled in the top-left panel.

Fig. A.4 Horizontal position-velocity cuts through bubble A from east to west derived from the [CII] line emission. The central panel goes through the estimated central position of A whereas the left and right panels are shifted 1′ to the north and south, respectively. The dashed line marks the approximate center of the bubble. The black dotted lines in the middle panel outline the elliptical shape a spherical shell with radius 100″ and expansion velocity of 14 km s−1 would have in the pv diagram. The dotted lines in the left and right panel correspond to the same expanding shell at the 1′ shift to the east and west where the 1′ corresponds to ~36 deg offset in the sphere (e.g., Butterfield et al. 2018, sketch in their appendix A.2). Hence, the velocities are then reduced by cos(36°).

Appendix B Finding bubbles

The assignment of the individual structures visible in the channel maps to bubbles and rings with a definite origin can be highly subjective. To avoid biases of an by-eye-identification, we implemented a method to quantify the amount of ring-like structures in individual maps. The goal of the method was to measure structures that one would identify by eye as rings and quantify their significance. This approach is complementary to ring-finding algorithms like CASI (Van Oort et al. 2019) that allow for huge amounts of data to be scanned, but that only provide a binary detection criterion.

To quantify the match between the map and a ring, we compute the covariance between the map and an annulus structure that is radially symmetric and positive between an inner and an outer radius and zero everywhere else. We normalize the positive value to provide an integral of unity. The covariance map peaks at the points with the strongest match between theactual structure and the ring. A search for the peaks therefore allows for ring finding, but a simple inspection of the covariance map is insufficient for this because individual peaks in the intensity distribution translate linearly intothe covariance coefficients. A reliable quantification therefore requires comparison of the match with different structures. We perform this through a parameter scan comparing structures with different ring radii and widths.

To maintain consistency with the typical by-eye-identification, we limited the ring width by a ratio between outer and inner radius between 1.1 and 1.6 and scanned this range. However, when applied to the data here, we found only a weak dependence of the results on this ratio and usually a good match for ratios around 1.4, so we limit our further discussion to rings of this width to reduce the dimensionality of the problem. Furthermore, we scanned possible ring radii between roughly 38 and 163″.

Figure B.1 demonstrates this effect for the [CII] channel map at − 61 km s⁻¹ from Fig. 6. The first panel shows the channel map emphasizing one of the prominent arcs with a radius of ~ 100″. The other panels show the covariance maps from the convolution with rings of 38, 57, 76, 96, and 115″. They are indicated as thin circles centered at the same position. We can see that for small radii, the covariance maps show a dip in the center of the circles, while only for the 96″ radius, the map shows a peak in the center. Considering only the location at the center of the circle and scanning the ring radii we find the maximum covariance coefficient for the “correct” ring radius. However, it is interesting to see that the covariance map for 76″ rings also shows a peak in the vicinity, shifted to the north by about 20″. This indicates that the arc in the channel map could also be fitted with a smaller ring that is spatially shifted. From these covariance maps only, neither of the two solutions can be excluded.

Fig. B.1 Demonstration of the covariance maps for the convolution of the −61 km s−1 channel map (see Fig. 6) with rings of different radii. The first panel shows the channel map. The magenta circle emphasizes one of the obvious ring-like structures with a radius of 100″. The other five panels show the covariance maps from the convolution of rings with increasing radii. The radii are indicated at the same position in the map.

Additional information can be obtained when including the velocity information if we do not only look for rings but the physical scenario of an expanding bubble as well. In that case, it is only at the systemic zero velocity that the bubble appears as a ring with itsfull diameter. At lower and higher velocities, the channel maps perform a cut through the bubble in front or behind the central plane, so that a smaller broadened ring appears. A different way to visualize this is by plotting the covariance coefficients against the velocity and the ring radius. For an ideal bubble, we can then find an inverse “V” shape where the largest ring radius is detected at zero systemic velocities. At intermediate velocities, we find a smaller ring.

To look for bubbles in the NGC7538 [CII] data, we can look for this inverse “V” shape in the covariance coefficients computed for all channel maps and ring radii. We performed a pre-selection of pixels by looking for rings in the individual channel maps. The principle is straightforward. We searched for peaks on the covariance maps that are small, isotropic, and sufficiently bright (for details see Kong et al. in prep.). Here, we performed a GAUCSCLUMPS (Stutzki & Guesten 1990) decomposition of the covariance maps and locate clumps that have an axis ratio of less than two, a radius below the inner ring radius of the convolution and an intensity of more than 5 K in the [CII] channel maps.

Fig. B.2 [CII] channel map with an overlay of all ring candidates identified in the individual channels. We plotted the rings not with their full width used in the convolution (outer radius 1.4 times the inner radius) but only as very thin rings to allow for a distinction of the different radii by eye.

The result is shown in Fig. B.2 for all [CII] channel maps. The number of ring candidates varies strongly between the different channels. At the extreme velocities, no rings are found, but the largest number of rings is not found at the velocities of the brightest intensity – but, rather, in panel 5, for a velocity of − 61 km s⁻¹. In many cases the radius of the ring is not accurately confined. We find overlapping rings with nearby centers. This especially the case for panel 5 (− 61 km s⁻¹), where five different radii can fit the central ring structure. A typical effect is best visible in panel 10, at − 51 km s⁻¹, where an arc like structure can be fitted by rings of different radii and simultaneously different centers. This demonstrates the fundamental uncertainty of the approach. Based on the plots here, it is only possible to measure the location of the center of each ring if we a priori know its radius. Fitting both leaves some uncertainty.

The locations of all these rings gives a guidance where to look for the inverse “V” structure in the four-dimensional covariance cube. If we look at the distribution of the different ring centers of the different channel maps (Fig. B.2), we find three dominant clusters that involve rings of varying radii:

• (A)

  R.A.(J2000.0)=23h13m38s, Dec=61°32′20″,

• (B)

  R.A.(J2000.0)=23h13m38s, Dec=61°30′12″,

• (C)

  R.A.(J2000.0)=23h13m46s, Dec=61°31′00″; a cluster with only small ring radii is located around:

• (D)

  R.A.(J2000.0)=23h13m55s, Dec=61°31′56″ and a cluster with only large radii around

• (E)

  R.A.(J2000.0)=23h13m39s, Dec=61°31′16″.

From the inspection of Fig. B.2, we can exclude most other positions because they only detected some outer arcs, not driven from inside of the region. In addition to that, the largest ring E appears as an overlap of the other smaller rings, in particular A, B and C (Figs. B.2 & B.4). Therefore, although identified by the algorithm, we do not consider E as a separate physical entity.

Fig. B.3 shows the covariance coefficients for these five locations as a function of channel velocity and ring radius. The relevant information is the location and width of the peak in the coefficients as a function of the ring radius, in particular the shape of the line formed by those peaks as a function of the velocity channel. All panels show somewhat different characteristics but also show at least some indication of the inverse “V” shape. Figure B.3 shows as a black line the maxima of the covariance coefficients that can be considered as a proxy for the ring or bubble radius. The width of these maxima can be considered as an estimator on the uncertainty on the ring radius. The “V” shape is most pronounced for structure D (bottom-left panel in Fig. B.3). There the maximum ring radius is ~80″ at velocities of − 57 and − 55 km s⁻¹ and we find dropping radii to higher and lower velocities. For the other bubble candidates, the transition to smaller radii and lower velocity is less well determined. In the following, we use as approximate radii for structures A, B and C 100″, and for structure D 80″.

Fig. B.3 Covariance coefficients at the positions of the five prominent clusters from Fig. B.2. The coefficients are shown as a function of the channel velocity and the ring radius used in the convolution. The black lines mark the maxima of the covariance coefficients in terms of the ring radius. In those channels without black line, no ring could be identified.

Figure B.4 shows the same result in a different way. Here, we plot the channel map from Fig. 6 again, but now with the over-plotted circles giving the mean ring diameter for the ring that is the local peak in the covariance map for the five points discussed in Fig. B.3. This means that for each channel, we selected the corresponding velocity column in the covariance plots in Fig. B.3, it is necessary to look for the maximum and draw the corresponding circle if the maximum is not at the minimum radius and above 4 % of the peak intensity. The five positions are color coded. The systematic behaviour of an increasing ring diameter with velocity is well visible for the grey and cyan rings representing position C and E. They behave like the ideal bubbleat least for blue-shifted velocities. In contrast the magenta and white rings, positions A and D, are relatively constant in diameter. They rather represent ring-structures in ppv-space. The matches at a velocity of − 51 km s⁻¹ rather look, at least partially, like a coincidence, not necessarily connected to the rings or bubbles at lower velocities. This all indicates that a bubble expansion in the region was hindered with regard to red-shifted velocities and it happens mainly at blue-shifted velocities, that is, toward the observer.

Fig. B.4 [CII] channel map with an overlay of the identified rings around the five candidate positions from Fig. B.3. The ring radii vary with velocity. The plotted rings are narrower than the rings used in the convolution. Different ring positions are encoded in a different ring color. The positions in panels 1–5 from Fig. B.3 correspond to the colors magenta (A), orange (B), grey (C), white (D), and cyan (E).

To estimate the accuracy of the determination of the ring positions and radii, we investigated the covariance plots for varying central position of the identified bubbles. This analysis indicates that for the conditions in NGC7538, with ring or bubble radii on the order of 60–160″, we can pin down their origin with an accuracy of 20–30″. The bubble diameter should be the largest one found in the central velocity channel. If we consider the rings from Fig. B.4 and compare the corresponding ring candidates in Fig. B.2, we see that in those velocity channels, there are typically two or three large rings that appear, corresponding to an uncertainty of the ring radius of 20–30″ (also see thelines in Fig. B.3).

References

All Figures

Fig. 1 Optical and cm continuum data of the NGC 7538 region, encompassing the field of view of our observations. The color scale shows the optical DSS2 (Digital Sky Survey v2) obtained with the red filter. The contours outline the 1.4 GHz cm continuum emission at 1′ angular resolution from the Canadian Galactic Plane Survey (CGPS, Taylor et al. 2003). The two stars in the center of the H II region mark the positions of IRS5 and IRS6 from Puga et al. (2010). The three five-pointed stars directly below the H II regions show the positions of the active sites of star formation IRS1, S, and IRS9.

In the text

Fig. 2 Integrated intensity molecular and atomic line maps as well as continuum emission are presented in color scale toward NGC 7538. The specific lines, integration ranges, and continuum bands are labeled in each panel. The corresponding angular resolution is shown at the bottom-left of the panels. The contours in most panels as labeled show the [CII] emission at only a few selected contour levels to highlight similarities and differences. The contours in the bottom-middle panel outline the 1.4 GHz cm continuum emission from the CGPS. The three five-pointed stars mark the positions of the embedded sources NGC 7538IRS1, NGC 7538S, and NGC 7538IRS9. The two stars in the middle show the positions of the two main exciting sources of the H II region NGC 7538IRS5 and NGC 7538IRS6. The bar-like PDR as well as a a linear scale bar are presented in the top-left panel.

In the text

Fig. 3 Overlays of selected pairs of emission tracers as described above each panel. The integration ranges for the atomic and molecular lines are the same as in Fig. 2. The contour levels are chosen individually for each figure to highlight the emission best. The angular resolution elements for both tracers are presented in the top-right of each panel. The five-pointed stars always mark the position of the embedded source NGC 7538IRS1 and the two open stars show the positions of the two main exciting sources of the H II region NGC 7538IRS5 and NGC 7538IRS6. A linear scale bar is presented in the bottom-left panel.

In the text

Fig. 4 Moment maps of the [CII] emission. Left, middle, and right panels: present (in color scale) the integrated emission and the first- and second-moment maps (intensity-weighted peak velocities and intensity-weighted velocity dispersion) in the velocity range from −75 to −40 km s−1, respectively.The noisier top-right quarters of the first- and second-moment maps are blanked (see Sect. 2). The contours in the middle and right panels show the integrated [CII] emission from the left panel at only a few selected contour levels to highlight similarities and differences. The three five-pointed stars in the left panel mark the positions of the embedded sources NGC 7538IRS1, NGC 7538S, and NGC 7538IRS9. The two stars show the positions of the two main exciting sources of the H II region NGC 7538IRS5 and NGC 7538IRS6. A linear scale bar is presented in the left panel. The markers in the middle and right panels highlight additional positions (see main text) around the areas where we extracted spectra in Fig. 5. The line in the middle panel outlines the pv-cut from IRS6 through the bar-like structure presented in Fig. 13.

In the text

	Fig. 5 [CII] (black), CO(3–2) (red), and [CI] (teal, multiplied by 2) spectra toward the positions marked in Fig. 4. The green lines show Gaussian fits to the [CII] and CO(3–2) data. The numbers present the corresponding FWHM values (“*/*” labels two components). The positions are labeled in each panel. The RPE position to the east is not covered by the [CI] map.
In the text

Fig. 6 [CII] channel map for the velocities labeled in each panel. The data are binned for that plot in 2 km s−1 channels. The two purple stars mark the positions of IRS1 and IRS5, and a scale bar is shown in the top-left panel. The intensities are always plotted logarithmically from 1 to 40 K. The purple circles outline potential bubbles. These bubbles are labeled in the top-left panel. A version without any suggested bubble structures is presented in Fig. A.2.

In the text

	Fig. 7 [CII] (black), CO(3–2) (red), and HI (blue, divided by 10) spectra toward three selected positions marked in Fig. 4. These are IRS5, the red-shifted peak toward the east (RPE), and the blue-shifted peak in the north (BlPN). The positions are labeled in each panel.
In the text

	Fig. 8 Potential bubble structures and 8 μm emission. The color scale shows the 8 μm data and the purple circles again outline the tentative bubble structures.
In the text

Fig. 9 Vertical position-velocity cuts through bubble A from south to north derived from the [CII] line emission. The central panel goes through the estimated central position of A whereas the left and right panels are shifted 1′ to the east and west, respectively. The white dashed line marks the approximate center of the bubble. The black dotted lines in the middle panel outline the elliptical shape that a spherical shell with a radius of 100″ and expansion velocity of 14 km s−1 would have in the pv diagram. The dotted lines in the left and right panel correspond to the same expanding shell at the 1′ shift to the east and west where the 1′ corresponds to ~36 deg offset in the sphere (e.g., Butterfield et al. 2018, sketch in their Appendix A.2); hence, the velocities are then reduced by cos (36°).

In the text

Fig. 10 Horizontal (along RA, left) and vertical (along Dec, right) position-velocity diagrams through the position of one of the exciting sources IRS5 (marked in each panel) derived from the [CII] line emission. corresponding roughly to the estimated center of bubble B. The contour levels are in 3σ steps. The black dotted lines outline the elliptical shape a spherical shell would have in the pv diagram – given a radius of 100″ and expansion velocities of 10 and 3 km s−1 at negative and positive relative velocities, respectively.

In the text

	Fig. 11 Horizontal (left) and vertical (right) position-velocity diagrams through the center of bubble D (marked in each panel) derived from the [CII] line emission. The contour levels are in 3σ steps. The black dotted lines outline the elliptical shape a spherical shell would have in the pv diagram – given a radius of 80″ and expansion velocities of 6 km s−1.
In the text

Fig. 12 Zoom on the bar-region south of the exciting sources IRS5 and IRS6. The color scale and the colored contours show the emission tracers labeled in the two panels. Only selected contours are shown for clarity. The four colored lines mark the approximate positions of the emission crests in the corresponding tracers. The sources IRS5, IRS6, and IRS1 are marked. A linear scale bar is shown as well.

In the text

	Fig. 13 Position-velocity cut from the [CII] line from IRS6 through the bar as outlined in the middle panel of Fig. 4. The white dashed line marks the approximate center of the bar.
In the text

Fig. 14 Integrated [CII] (color), [CI] (black contours), and CO(3–2) (green contours) emission in the bar region where the [CII]/[CI]/CO mass estimates are conducted (Sect. 4.3). All three maps were integrated between −70 and −40 km s−1. The [CI] contours start at 6σ and continuein 3σ steps (1σ = 1.6 K km s−1). The CO(3–2) contours start at 10σ and continuein 10σ steps (1σ = 3.2 K km s−1). The white line marks the approximate location of the bar-like feature.

In the text

	Fig. A.1 [CII] (black) and CO(3–2) spectra toward the grid marked in the top-left [CII] first-moment map. The contours show the [CII] emission at only a few selected contour levels to highlight the main structures. The sources IRS5, IRS6, and IRS1 are also shown as black-and-white stars.
In the text

	Fig. A.2 [CII] channel map for the velocities labeled in each panel. The data are binned for that plot in 2 km s−1 channels. The two purple stars mark the positions of IRS1 and IRS5, and a scale bar is shown in the top-left panel. The intensities are always plotted logarithmically from 1 to 40 K.
In the text

	Fig. A.3 CO(3–2) channel map for the velocities labeled in each panel. The data are plotted in steps of 2.1 km s−1 (5 channels). The two purple stars mark the positions of IRS1 and IRS5, and a scale bar is shown in the top-left panel. The intensities are always plotted linearly from 0.1 to 10 K. The purple circles outline potential bubbles. These bubbles are labeled in the top-left panel.
In the text

Fig. A.4 Horizontal position-velocity cuts through bubble A from east to west derived from the [CII] line emission. The central panel goes through the estimated central position of A whereas the left and right panels are shifted 1′ to the north and south, respectively. The dashed line marks the approximate center of the bubble. The black dotted lines in the middle panel outline the elliptical shape a spherical shell with radius 100″ and expansion velocity of 14 km s−1 would have in the pv diagram. The dotted lines in the left and right panel correspond to the same expanding shell at the 1′ shift to the east and west where the 1′ corresponds to ~36 deg offset in the sphere (e.g., Butterfield et al. 2018, sketch in their appendix A.2). Hence, the velocities are then reduced by cos(36°).

In the text

Fig. B.1 Demonstration of the covariance maps for the convolution of the −61 km s−1 channel map (see Fig. 6) with rings of different radii. The first panel shows the channel map. The magenta circle emphasizes one of the obvious ring-like structures with a radius of 100″. The other five panels show the covariance maps from the convolution of rings with increasing radii. The radii are indicated at the same position in the map.

In the text

	Fig. B.2 [CII] channel map with an overlay of all ring candidates identified in the individual channels. We plotted the rings not with their full width used in the convolution (outer radius 1.4 times the inner radius) but only as very thin rings to allow for a distinction of the different radii by eye.
In the text

	Fig. B.3 Covariance coefficients at the positions of the five prominent clusters from Fig. B.2. The coefficients are shown as a function of the channel velocity and the ring radius used in the convolution. The black lines mark the maxima of the covariance coefficients in terms of the ring radius. In those channels without black line, no ring could be identified.
In the text

Fig. B.4 [CII] channel map with an overlay of the identified rings around the five candidate positions from Fig. B.3. The ring radii vary with velocity. The plotted rings are narrower than the rings used in the convolution. Different ring positions are encoded in a different ring color. The positions in panels 1–5 from Fig. B.3 correspond to the colors magenta (A), orange (B), grey (C), white (D), and cyan (E).

In the text