© The Authors 2025

1 Introduction

In the vicinity of massive young stars (>8–10M_⊙), molecular clouds are exposed to extreme ultraviolet (UV) radiation and ionized gas pressure causing compression between molecular gas and ionized gas regions (Zinnecker & Yorke 2007). The far ultraviolet radiation (FUV, 6 e V < hν < 13.6 eV) of stars, irradiating the nearby molecular gas creates photodissociation regions (PDRs, Tielens & Hollenbach 1985; Tielens 2013).

Strong, energetic injections reshape parental clouds, triggering low (and intermediate) and high mass star formation (e.g., Hester et al. 1996; Sugitani et al. 2002; Schneider et al. 2021; Comerón et al. 2022) and also contributing to the complexity of the PDR chemistry. This makes the immediate neighborhood of OB stars an important laboratory for studying the formation and destruction of molecular species.

The creation of substructures at the interface between an expanding HII region and its native molecular clouds is clear evidence of stellar feedback. Prominent substructures may be comprised of pillars resembling column-like features that are still connected to their natal molecular clouds, with sizes ranging from 0.5 pc up to a few pc. Unlike pillars, globules are isolated and have a cometary-like shape with a head-tail morphology (e.g., Schneider et al. 2016a). These pillars and globules mainly host low-mass star formation (e.g., Hester et al. 1996; White et al. 1999), while only a few among them may contain high-mass or intermediate-mass early B-type stars (e.g., Schneider et al. 2012; Djupvik et al. 2017; Comerón et al. 2022; Schneider et al. 2021). (Magneto)-hydrodynamic simulations have demonstrated that turbulence in UV-irradiated environments plays a crucial role in the formation of pillars and globules (Gritschneder et al. 2010; Tremblin et al. 2012a,b).

Schneider et al. (2016a) characterized and classified pillars, globules, and evaporating gaseous globules (EGGs), which are in the form of small globules, along with proplyd-like objects resembling evaporating circumstellar disks, according to their physical properties determined with the Herschel far-infrared observations at 70–500 µm toward the Cygnus OB2 (Cyg OB2) association, under the Herschel imaging survey of OB Young Stellar objects (HOBYS) program (Motte et al. 2010). The Cyg OB2 association significantly affects the surrounding medium, especially the nearby molecular clouds, by forming pillars and globules (Schneider et al. 2016a). DR18 is one of the largest globules in the Cygnus-X region and is located close to the OB2 association region, which contains 169 primary OB stars (age ~1–7 Myr, Wright et al. 2015) and it is clearly visible in the Spitzer/IRAC 8 µm emission shown in Fig. 1. Schneider et al. (2006) also show that DR18 is associated with the Cygnus-X region according to the CO observations and, thus, for the DR18 distance, we adopted a distance of 1.4 kpc determined by trigonometric parallaxes of 6.7 GHz CH₃ OH masers (Rygl et al. 2012), which is also consistent with the results of Dzib et al. (2013) measuring the VLBI parallax toward Cyg OB2 #5. This globule is exposed to the high-UV radiation field (~1000 G₀¹, Schneider et al. 2016a). Despite such strong UV radiation, this region contains several young stellar objects (YSOs) (Wright et al. 2015; Comerón et al. 2022). In addition, its surroundings are filled by low-density ionized gas having filament-like structures (Piddington & Minnett 1952; Emig et al. 2022; see Reipurth & Schneider 2008 for an overview). This makes the DR18 region a perfect laboratory to investigate how stellar feedback, such as turbulence and UV radiation, influences the chemical complexity of a molecular cloud and its evolution.

Previous observations of Cygnus-X have mostly focused on continuum emission from submillimeter (submm) (e.g., Motte et al. 2007; Schneider et al. 2016a; Cao et al. 2019) to centimeter (cm) wavelengths (e.g., Brunthaler et al. 2021; Emig et al. 2022). Most of the molecular line observations were done with a low angular resolution (>15″) for this environment, mainly from CO surveys with single dish telescopes (Schneider et al. 2006; Gottschalk et al. 2012), except for the [OI] observations with a better resolution of 6″ (Schneider et al. 2021). Several high angular resolution interferometric observations targeting molecular outflows (Duarte-Cabral et al. 2013, 2014; Skretas & Kristensen 2022), 6.7 GHz CH₃ OH masers (Ortiz-León et al. 2021), and H₂CO absorption line (Gong et al. 2023), the latter two studies as part of the GLObal view on STAR formation in the Milky Way (GLOSTAR) survey (Brunthaler et al. 2021). Thus, more research needs to be carried out on the dense, cold gas components of DR18 and Cygnus-X as a whole with high angular resolution to fully understand the influence of the stellar feedback and UV radiation from the OB2 association.

This work makes use of the results of observations made in the framework of the Cygnus Allscale Survey of Chemistry and Dynamical Environments (CASCADE; see Beuther et al. 2022 for an overview), which provides a unique view of molecular gas components toward the Cygnus-X region, with high angular resolution millimeter (mm) observations of different molecular gas tracers. The observations were carried out at wavelengths of 3–4 mm, centered at ≈3.6 mm (or frequencies of approximately 70.5–78.2 GHz and 85.2–93.9 GHz) with unprecedentedly high angular resolutions (3″–4″) over an area of 5.1′ × 5.1′. The spatial resolution is approximately 0.02– 0.03 pc (or ~4000–6000 au) at an adopted distance of 1.4 kpc. In Sect. 2, we describe the observations obtained within the CASCADE program and aspects related to the data reduction process. Supplementary data sets for comparisons with our 3.6 mm data and a description of DR18 are presented in Sect. 3. In Section 4, we compare the millimeter and centimeter (mm and cm) continuum emission datasets and present molecular gas structures traced by HCO⁺, N₂H⁺, SiO, and NH₂D. In that section, we describe compact core identification and the physical properties of the identified cores. In Section 5, we discuss the origin of the free-free radiation contribution to the 3.6 mm continuum emission to understand the ionized gas structures toward DR18 and the stellar feedback to the surrounding dense molecular gas. Furthermore, we investigate the association between SiO emission and the PDRs of DR18. Lastly, we summarize our findings and emphasize the main conclusions in Section 6.

Fig. 1 Spitzer/IRAC 8 µm emission map (in a flux range of 0–98 MJy/sr) of a part of the Cygnus-X complex (Beerer et al. 2010). The yellow rectangle outlines DR18, which is studied here. The white circle has a radius 12.2 pc (at a distance of 1.4 kpc), with its center at the position of the trapezium (Cyg OB2 #8, RA 20h33m 16s and Dec +41°18′45″) of O stars. The cyan and bright green markers are known members of O- and B-type stars, respectively, in the OB2 association (Wright et al. 2015).

2 Observation and data reduction

2.1 IRAM 30 m and NOEMA

The observations of the DR18 region were done with the Northern Extended Millimeter Array (NOEMA) for high angular resolution imaging and the Institut de Radioastronomie Mil- limétrique (IRAM) 30 m telescope for lower angular resolution imaging. The NOEMA observations comprised one mosaic tile consisting of 78 pointings in the C and D configurations between December 2019 and April 2020. The center coordinates of the covered field of view (FoV) of 5.1′ × 5.1′ are RA 20^h35^m08.370^s and Dec +41°13′30.40″ (l, b = 080.3622, +00.4479). With the goal to restore missing fluxes of extended molecular line resolved out by the interferometer, we used the IRAM 30 m telescope to cover the 5.1′ × 5.1′ area, in on-the-fly (OTF) mode, which resulted in typical root mean square (rms) noise values of ~0.1 K on the T_mb scale. A total bandwidth of 16 GHz is covered in dual polarization, the lower sideband covering 70.4 GHz to 78.2 GHz, and the upper sideband from 85.8 GHz to 93.6 GHz. The 3.6 mm continuum emission is obtained over the whole bandwidth with a spectral resolution of 2.0 MHz, equivalent to a velocity resolution of ~7.3 km s^–1 after masking and discarding portions of the frequency range with emission or absorption lines. In the NOEMA observations, a number of individual spectral lines targeted for the CASCADE program were covered with the higher spectral resolution of 62.5 kHz, yielding velocity resolutions between 0.26 and 0.19 km s^–1, while for the IRAM 30 m observations, the spectral resolution is consistently resampled to ~0.8 km s^–1. Since this study uses the NOEMA plus IRAM 30 m combined data, we note that all analyses for spectral lines are limited to a velocity channel width of 0.8 km s^–1. Table 1 lists all the molecular lines detected toward DR18 and the spectroscopic parameters of the lines, along with information on the synthesized beams, as well as noise levels of 1 σ_rms. All the data reductions, including calibration and imaging, for the NOEMA and 30 m telescope were performed with the Grenoble Image and Line Data Analysis Software (GILDAS) package (Pety 2005). The further details of the observations with NOEMA and the IRAM 30 m telescope and data reduction procedures for the CASCADE program are described in the CASCADE overview paper (see Beuther et al. 2022).

2.2 GBT MUSTANG-2

To complement also the mm continuum data with short-spacing information, we used the MUSTANG-2 array on the Green Bank Telescope (GBT). The MUSTANG-2 instrument at GBT has a 4.2' instantaneous FoV comprised of a 215-element array of feedhorn-coupled TES bolometers (Dicker et al. 2014). It has a 30 GHz bandpass with an effective bandpass center near 90 GHz (Ginsburg et al. 2020). Observations for this work were performed under project GBT22A-280 on 2022 February 20 and March 14, with 3.7 hours spent on source (spread over multiple pointings). The data were calibrated and reduced with the MUSTANG IDL Data Analysis System (MIDAS), with improvements in the iterative mapmaking pipeline (e.g., see Appendix C in Romero et al. 2020, for details of the pipeline). The primary improvements are the reduction of ringing, which is inherent in the reduction pipeline, as well as the improved recovery of larger scales.

2.3 Merging millimeter continuum emission data obtained with NOEMA and MUSTANG-2

The bandpass of the MUSTANG-2 camera on the GBT roughly covers the frequency range between 75 and 109 GHz. With 90 GHz being almost exactly at the center of that, we use for the combination of the data only the NOEMA USB data that are centered around 90 GHz as well. We acknowledge that the bandpass widths of MUSTANG-2 and NOEMA are different (30 versus 8 GHz), but that the difference is not accounted for in the combination. Having the same central frequency should make the two datasets compatible in good approximation. To combine NOEMA and MUSTANG-2 data, we tried different approaches employing both the Common Astronomy Software Applications (CASA) and GILDAS packages. The results obtained by the two approaches agreed very well. Therefore, in the following, we describe and use the data combined in GILDAS. The MUSTANG-2 data have a beam size of 9″ and a pixel size of 1″. We first reprojected the MUSTANG-2 data to the phase center of the NOEMA data and updated a few MUSTANG header parameters to the GILDAS format. For example, the rest frequency or beam size have to be explicitly given in the MUSTANG-2 data header for proper data combination. The final data combination was then conducted with the GILDAS task uv_short using a uv truncation radius of the MUSTANG-2 data of 15 m. The synthesized beam of the combined data is 3.″32 × 2.″66 with a position angle of –165° (E of N). The 1σ_ms noise of the combined data is 0.07 mJy beam^–1.

3 Ancillary and archival data

To better understand the environment of the DR18 region, we utilized ancillary continuum data from various surveys carried out at different wavelengths. To trace the ionized gas, we used the high angular resolution radio continuum emission at 6 cm wavelength from the GLOSTAR survey⁴ (Medina et al. 2019; Brunthaler et al. 2021). The 6 cm radio continuum data were generated by combining the Karl G. Jansky Very Large Array (VLA) D- and B-configuration data, which yields a synthetic beam size, θ, of 4″ and a sensitivity of ~0.1 mJy beam^–1.

In addition, to investigate the properties of emission from dust and polycyclic aromatic hydrocarbons (PAHs) tracing PDRs, we compared the Spitzer Legacy Survey of the Cygnus-X complex (Cygnus-X; Hora et al. 2007; Beerer et al. 2010), which mapped a 6° × 6° area with Spitzer/IRAC four bands (3.6, 4.5, 5.8, and 8.0 µm; Papovich et al. 2016). We also used H₂ column density and dust temperature maps obtained from Bonne et al. (2023), adopting the methods described in Palmeirim et al. (2013). The H₂ column density is determined with the difference maps of the convolved maps from 250 µm to 500 µm from the HOBYS program (Motte et al. 2010) and the dust temperature map is generated with the 160 µm/250 µm flux ratio. Lastly, we utilized 2.166 µm Brγ line and 2.122 µm H₂ (1–0) S(1) vibrational line obtained from Comerón et al. (2022) to identify PDRs in DR18 by comparing our NOEMA data and other ancillary datasets. Table 2 summarizes the continuum emission datasets and line ancillary data.

4 Results and analysis

The head of the DR18 globule observed in the CASCADE survey is visually striking in infrared emission, particularly in the Spitzer/IRAC 8 µm wavelength band, which is dominated by radiation from PAHs. The head faces the main Cyg OB2 population (see Fig. 1) and exhibits excellent spatial alignment between the highly saturated 8 µm emission and strong 6 cm radio continuum emission obtained from the GLOSTAR survey (as shown in the larger FoV image in the left panel of Fig. 2). Also, this mid-infrared emission image reveals an apparent hollow feature, such as an ear-like shape (delineated by a yellow curved line). In addition, from the distinct appearance of the infrared emission, several YSOs have been identified and classified by Comerón et al. (2022), indicating active star formation inside of the globule.

A group of Class II and transition-disk YSOs is clearly situated in the HII region traced by the 6 cm radio continuum emission. In the vicinity of the curved illuminated 8 µm emission (yellow curve line), a Class I YSO (yellow square) and a flat spectrum YSO (red upside-down triangle) are found, which are younger than the other YSOs in this region. Further out toward the tail of the globule, only more evolved YSOs have been found (e.g., Class III YSOs, purple circles). A recent nearinfrared spectral study conducted by Comerón et al. (2022) found the spectral type of the transition disk YSO (DR18-05) close to the 6 cm radio continuum emission peak to be aB2 star, which is thought to have created the cavity in the head of DR18 (Comerón & Torra 1999; Comerón et al. 2022). According to Comerón et al. (2022), the B star could have required approximately 3 × 10⁴ years to carve the arc-shape with a radius of 0.22 pc on the tip of DR18.

Fig. 2 Three-color composite images of DR18 (8 µm in red, 4.5 µm in green, and 3.6 µm in blue) with different maximum flux levels for three Spitzer/IRAC bands to show the extended emission in a larger FoV (left panel) and the emission details, which are saturated in the left panel, in a zoom-in view (right panel) covering our 5.1′ × 5.1′ target region (outlined by a white square box in the large FoV image). In both images, cyan contour lines represent the combined VLA B+D array 6 cm radio continuum emission from the GLOSTAR survey by levels of 0.35, 0.70, 1.05, 1.40, 1.75, 2.10, 2.45, 2.80, 3.15, and 3.50 mJy/beam. In the zoom-in image, the different color symbols indicate young stellar objects identified by infrared colors (Comerón et al. 2022): yellow squares mark Class I YSOs, light blue diamonds for Class II YSOs, red upside-down triangle for a flat spectrum YSO, purple circles for Class III YSOs, and bright green upside triangles represent transition-disk YSOs. The yellow half-arc line indicates the open shell or the ear-like feature in the infrared emission image.

4.1 Dust thermal and free-free emission contributions to 3.6 mm continuum emission

Figure 3 shows the 3.6 mm continuum and 6 cm continuum emission maps from this study and the GLOSTAR survey (from left to right). While the 3.6 mm and the 6 cm continuum emission are similarly extended, it is notable that the 6 cm emission extends more to the west compared to the 3.6 mm emission. On the other hand, many of the bright, compact features on the 3.6 mm continuum emission do not have significant 6 cm emission counterparts. The overall emission distributions at these two wavelengths imply that the 3.6 mm continuum emission potentially has significant contributions from free-free emission originating in the ionized gas, apart from the dust thermal emission.

To examine the contribution of free-free radiation emission to the 3.6 mm continuum emission, we analyzed spectral indices between 6 cm and 3.6 mm wavelengths using the 6 cm GLOSTAR combined the D+B VLA array continuum image and the 3.6 mm MUSTANG-2+NOEMA combined continuum image. To compare both datasets, we re-convolved the MUSTANG-2+NOEMA combined data (3732 × 2766 with a position angle of –165°) to the exact angular resolution (4″) of the VLA D+B combined data, and then resampled the reconvolved MUSTANG-2+NOEMA data to match the pixel size of the VLA data. This spectral index measurement consists of fitting a simple slope to the continuum fluxes at the central frequencies of the two datasets (5.8 GHz for the GLOSTAR survey and 90 GHz for the CASCADE survey) in logarithmic scale, given by the formula $${\alpha }_{90\text{\hspace{0.17em}GHz}-5.8\text{\hspace{0.17em}GHz}}=\frac{\log \left(S_{90\text{\hspace{0.17em}GHz}}/S_{5.8\text{\hspace{0.17em}GHz}}\right)}{\log (90\text{\hspace{0.17em}GHz}/5.8\text{\hspace{0.17em}GHz})}.$$(1)

The S_{90 GHz} and S_{5.8 GHz} are the continuum flux densities per pixel. The spectral index uncertainty is the propagating error estimated as $${\sigma }_{90\text{\hspace{0.17em}GHz}-5.8\text{\hspace{0.17em}GHz}}=\frac{\sqrt{{\left(3{\sigma }_{S_{90\text{\hspace{0.17em}GHz}, \text{rms}}}/S_{90\text{\hspace{0.17em}GHz}}\right)}^2+{\left(3{\sigma }_{S_{5.8\text{\hspace{0.17em}GHz}, \text{rms}}}/S_{5.8\text{\hspace{0.17em}GHz}}\right)}^2}}{\log (90\text{\hspace{0.17em}GHz}/5.8\text{\hspace{0.17em}GHz})}.$$(2)

The ${\sigma }_{S{}_{90\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}GHz}}}$ and ${\sigma }_{S_{5.8\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}GHz}}}$ are the continuum flux uncertainties corresponding to the frequencies. Instead of using 1 σ_rms observational uncertainty for the spectral index uncertainty, we applied 3 σ_rms threshold because the GLOSTAR continuum emission data used for the spectral index does not include the total-power emission data, which is more sensitive to extended emission structures. This concerns missing fluxes for faint, extended emission in the outer regions of the cm continuum emission image. To mitigate this, we restricted our analysis to regions with emission above 3 σ_rms at 5.8 and 90 GHz, respectively.

Figure 4 displays the 90 GHz – 5.8 GHz spectral index map in the left panel and its spectral index error map in the right panel. The determined spectral index in large parts is consistent with free-free emission close to a spectral index of α = −0.1 with several compact regions that show a clear dust thermal contribution with positive spectral indices. While the western edge shows slightly more negative values than those values expected from dust thermal emission and free-free emission regions, these parts also have larger uncertainties approaching 0.2. Furthermore, these values must be taken cautiously because of uncertainties introduced by different UV coverages of both observations. In general, the spectral index is between −0.1 and 2.0 for typical HII regions and in a range of 2.0 < α < 4.0 for sources with pure dust thermal emission (e.g., Sánchez-Monge et al. 2017). We note that the compact regions with thermal emission are found exclusively east of the B2 star, while the western part of the emission, which faces the OB2 association, is dominated by free-free emission.

Fig. 3 Millimeter (left) and centimeter (right) continuum emission maps. The mm continuum emission map represents the combination of the GBT MUSTANG2 and NOEMA data sets, whereas the cm continuum emission map is generated by combining the VLA D and B configuration observations. The FWHM synthetic beams of these maps are displayed in the lower left corners. The contour levels of the mm continuum emission are 0.46, 0.96, 1.38, 1.80, 4.14, 5.33, and 6.51 mJy beam–1, whereas for the cm emission, the levels are 0.35, 0.70, 1.05, 1.40, 1.75, 2.10, 2.45, 2.80, 3.15, and 3.50 mJy beam–1. The star symbol indicates the position of the B2 star.

Fig. 4 90 GHz – 5.8 GHz spectral index map in the left panel and its spectral index error map in the right panel. The higher frequency data are smoothed to the beam size of the lower frequency data to achieve a common resolution of 4″. The star symbol indicates the position of the B2 star. The black contours are the 3.6 mm continuum flux, and their levels are the same as in the left panel of Fig. 3, and the blue contours are the 6 cm continuum emission with the same contour levels shown in the right panel.

4.2 Molecular line profiles and emission distributions

The CASCADE observations have provided, for the first time, a high angular resolution view of spectral line emission in the DR18 region, allowing for detailed investigations of its molecular gas content to be undertaken on physical scales of ∼0.02 – 0.03 pc. Figure 5 shows the velocity-integrated maps of HCO⁺, N₂H⁺, SiO, and NH₂D emission lines with the overlaid mm continuum emission in black contours. The intensity maps of the other species detected in this region are presented in Figs. A.1 and A.2. The HCO⁺ intensity integrated over a velocity range from 6km s⁻¹ to 15.6 km s⁻¹ shows an extended distribution and exhibits the arc-shape morphology resembling the bright 8 µm emission feature indicated by a yellow arc line in the right panel of Fig. 2. In addition, weaker HCO⁺ emission extends toward the tail of the DR18 visible in the infrared emission. Such spatially widespread emission distributions are also found in the intensity maps of CCH, H₂CO, HCN, and HNC exhibited in Fig. A.1. The outline of HCO⁺ line emission follows the morphology of the extended 3.6 mm and 6 cm continuum emission regions. Such spatial coincidence features between molecular gas and free-free emission are also found toward the pillars in other OB association regions, such as the elephant trunks in M16 (Sofue 2020), as well as toward an ultracompact HII region Mon R2 (Treviño-Morales et al. 2016). The bright HCO⁺ emission peaks spatially match with the 3.6 mm compact continuum emission features, which also coincide significantly with bright features in the N₂H⁺ and NH₂D maps.

In contrast, N₂H⁺ line emission better coincides with the compact mm continuum emission (see the upper right panel in Fig. 5) rather than with pure free-free emission found at cm wavelengths. This supports that the continuum emission features traced by dense molecular tracers (n_cr > 10⁵ cm⁻³) represent dust thermal emission rather than free-free emission. These dense gas regions also have counterparts in the emission from other molecular lines (as shown in Fig. A.2), including lines from deuterated species; for example, from NH₂D in the right bottom panel in Fig. 5. Interestingly, some of the compact emission components are aligned as a “crown-like” feature, indicated by the magenta- dotted curve opening toward the northern direction in the N₂H⁺ intensity map. The feature is also pronounced in the HNC, HN¹³, HC₃N, and H¹³CO⁺ emission maps. The eastern portion of the crown-like feature does not show any continuum emission counterpart and only appears in those species. NH₂D, DCN, and DNC shown in Figs. 5 and A.2 are also detected toward the dense molecular gas regions. However, DCN shows weaker emission and is spatially less matched with DNC and NH₂D emission.

The left bottom panel of Fig. 5 displays the SiO emission integrated over a velocity range from 8 km s⁻¹ to 10km s⁻¹. The SiO emission has a remarkably narrow linewidth, considering that SiO emission is often used as a tracer of shocked gas and this type of emission has very broad line wings (e.g., Martin-Pintado et al. 1992; Bachiller & Pérez Gutiérrez 1997; Codella et al. 1999; Csengeri et al. 2016). Not only does SiO have a narrow linewidth, but its emission shows clumpy distributions concentrated in spatially compact regions around the HII region of D18-05. Compared to other molecules, the distribution of SiO emission is distinct and spatially separated from that of other dense gas tracers. The SiO emission peaks do not coincide with 3.6 mm continuum emission peaks or the bright N₂H⁺ emission components and appear in regions without dense gas, unlike other dense gas-tracing molecular lines, which are preferably detected toward the mm continuum emission peaks. There are only a few weak SiO detections toward the compact mm features. The north of the HII region shows slightly more extended SiO emission. Figure 6 shows the intensity ratio map of SiO over N₂H⁺ lines. It clearly shows that the intensity ratios increase outwards of the N₂H⁺ emission regions (close to the black contours), indicating the relatively weaker N₂H⁺ and brighter SiO in the regions. On the other hand, in the inward parts of the N₂H⁺ emission regions, the ratios decrease. This pattern appears consistently in all the areas having detections of both species. We explore the origin of this narrow line SiO emission in Section 5.2.

4.3 Gas temperature distribution

Based on the dust temperature map (shown in the upper panel of Fig. 7) taken from Bonne et al. (2023), we know that toward the interior HII region of DR18 the dust temperatures are higher than 23 K and can reach up to 28 K. However, due to the larger angular resolution of 18″ (Bonne et al. 2023) than the NOEMA synthesized beams, it is difficult to investigate the temperature structure toward the N₂H⁺ cores. However, it is still evident that the cores are relatively colder than the interior of the HII region, given the large spatial resolution. As these cores are expected to be denser, colder parts in the DR18 region, we expect that they have lower gas temperatures than those obtained from the dust temperature map. To determine the temperature distributions of extended and dense molecular gas components, we used the line emission ratios of HCN (1−0) and HNC (1−0) and their isotopologues, as suggested by Hacar et al. (2020). The authors show that the ratio of these ground state transitions of the two species for values from 1 to 4 (corresponding to temperatures between 10 and 40 K) provides a good approximation of gas kinetic temperature, as derived from the NH₃ inversion line data toward the Orion molecular cloud. Unlike other molecular transitions, the HCN and HNC lines have similar distributions as the HCO⁺ emission and ridge structures in the N (H₂) map, with N (H₂) > a few ×10²² cm⁻². To generate HCN/HNC 1−0 ratio maps, we smoothed HCN, HNC, and HN¹³C maps with the beam size of H¹³CN, which has the largest beam size among these lines. After the spatial smoothing, we sequentially performed regridding all the smoothed data to the same pixel size and coordinates of H¹³CN. From the HCN/HNC and H¹³CN/HN¹³C maps, we generated gas temperature maps following equations presented in Hacar et al. (2020): $$T_{\text{HCN}/\text{HNC\hspace{0.17em}}}(\text{K})=10\times \left[\frac{I(\text{HCN})}{I(\text{HNC})}\right]\text{ if }\frac{I(\text{HCN})}{I(\text{HNC})}\le 4$$(3) $$T_{\text{HCN}/\text{HNC\hspace{0.17em}}}(\text{K})=3\times \left[\frac{I(\text{HCN})}{I(\text{HNC})}-4\right]+40\text{ if }\frac{I(\text{HCN})}{I(\text{HNC})}>4.$$(4)

The bottom panel of Fig. 7 displays a map of the gas temperatures in high angular resolution (~3″) obtained by combining both gas kinetic temperature maps (see Fig. B.1) determined with the HCN/HNC and H¹³CN/HN¹³C line ratios. The combination was performed by substituting the temperatures of HCN/HNC with those of H¹³CN/HN¹³C. The study by Hacar et al. (2020) mentions that the expected uncertainty with this method is ΔT_HCN/HNC ≈ 5 K for temperatures ranging from 10 < T_HCN/HNC < 40 K and ΔT_HCN/HNC ≈ 10 K for temperatures above 40 K. The effect of optical depth for HCN and HNC can be ignored because the highest column density of hydrogen molecules toward the DR18 region is not above 4.2 × 10²² cm⁻², which corresponds to a visual extinction, A_V, of 42 mag calculated using the A_V = N(H₂)/9.4 × 10²⁰ cm⁻². The study by Hacar et al. (2020) states that for regions with up to A_v ≈ 100, their method provides HCN/HNC line ratios that are still acceptable. Nevertheless, toward higher H₂ column density regions, we estimated the gas kinetic temperatures using H¹³CN and HN¹³C emission lines, whose optical depths are smaller than 0.2–0.3.

Figure 7 shows the kinetic temperature map determined from HCN/HNC and H¹³CN/HN¹³C ratios. In the T_HCN/HNC map, the east of the HII region and the crown-like structure appearing in the N₂H⁺ intensity map (Fig. 5) show colder temperatures 5 K ≲ T_HCN/HNC ≲ 25 K, whereas the T_dust map shows slightly higher temperatures of 25–27 K. In contrast to the dense regions, the extended molecular gas regions show higher temperatures as high as above 40 K in the T_HCN/HNC map compared with those ~22 K–26 K in the T_dust map. This temperature difference occurs because the far-infrared wavelengths between 160 µm and 500 µm used for determining the dust temperature map are not sensitive to warmer gas above 40K. In addition, higher gas temperatures have been observed toward the extended features and between the N₂H⁺ cores, which are considered to be the lower gas density regions and highly UV-illuminated regions as seen in Fig. 2. Thus, in such regions, dust and gas could be decoupled as observed in other HII and photodissociation regions (e.g., Koumpia et al. 2015; Salgado et al. 2016). However, as mentioned above, in these regions with T_HCN/HNC > 40 K, the uncertainties are higher about 10 K. Nevertheless, the high-temperature values in the T_HCN/HNC map could be unlikely to be artifacts and might imply an interaction between the ionized gas and molecular gas, with careful considerations of the absolute gas temperature values.

Fig. 5 Velocity-integrated intensity maps of HCO+ (6km s−1 to 15.6km s−1), N2H+ (−1.2km s−1 to 17.2km s−1 spanning all the hfs transitions) in the upper row and SiO (8km s−1 to 10km s−1) and NH2D (0km s−1 to 15.2km s−1 covering all the hfs transitions) in the lower row. The black contours represent the GBT+NOEMA combined 3.6 mm continuum emission, and the bright green marker indicates the position of the ionizing source, DR18-05. The contour levels of the mm continuum are the same as shown in Fig. 3. The horizontal white scale bar has a size of 0.1 pc at a given distance of 1.4 kpc. The beam of the maps of each species is displayed on the left corner of each image panel.

Fig. 6 Intensity ratio map of SiO/N2H+ in color scale for the overlapping areas of both molecular species. The black and pink contours represent the 5σ level of N2H+ velocity-integrated intensity and a 3σ level of SiO velocity-integrated intensity, respectively. The bright green star indicates DR18-05, the ionizing source in DR18.

Fig. 7 Dust and gas temperature maps. Top: dust temperature map obtained from Bonne et al. (2023) generated with the 160 µm/250µm flux density ratio as described in the work by Palmeirim et al. (2013). Bottom: gas kinetic temperature map combined temperature maps derived from the HCN/HNC and H13CN/HN13C emission line ratios. The black contours represent the 3.6 mm continuum emission. The bright green star indicates the position of the B2-type star D18-05. The beams for the Tdust and THCN/HNC maps are 18″ and 3.43 × 2.76 with a position angle of 13.66° displayed in the lower left corner.

4.4 Compact core identification

We utilized the N₂H⁺ intensity map, integrated over a velocity range encompassing all the hyperfine structure (hfs) transitions, to extract compact cores that trace the dense gas regions. To identify substructures (e.g., cores) in DR18, we used the astrodendro package⁶ which applies the Dendrogram algorithm (Rosolowsky et al. 2008). Among the substructures classified by astrodendro, “leaves” are the smallest substructures referred to here as compact cores. Three input parameters were required for the astrodendro package: min_value, min_delta, and min_npix. We applied a 5σ threshold value to min_value, where the 1σ_rms noise (0.049 Jy beam⁻¹ km s⁻¹) is measured from the velocity-integrated intensity map. Additionally, we chose min_delta for the significance of a leaf compared with its neighboring leaf or branch features to be equal to 1 σ_rms . Lastly, we used min_npix as the minimum number of pixels to identify a leaf, set as 45 pixels, the number of pixels covered by the beam for N₂H⁺ emission. The leaf structures are defined as brighter than the five σ_rms level and larger than the beam size. With these input parameters, we identify 18 leaf substructures that we hereafter consider as cores. In Fig. 8, the black contours indicate the substructures, and the blue ellipses show their schematic structures, which are not extracted source sizes. The central coordinates of the structures are their mean positions in the x and y directions and are listed in the second and third columns of Table 3.

4.5 Derivation of the core physical parameters

Here, we determine the physical quantities of the N₂H⁺ cores identified in the previous section. The effective core radius is derived as $R_{\text{core }}=\sqrt{A/\pi }$, where A is the area of the leaves on the plane of the sky and adopting a distance of 1.4 kpc. The average radius of the cores is 0.025 pc. The obtained core radius in pc and its diameter in arcseconds are listed in Table 3. With the given core radius, we extracted the N₂H⁺ spectra toward the 18 compact cores to derive N₂H⁺ column densities and to determine line parameters (e.g., linewidths).

To deconvolve the hyperfine structure of the N₂H⁺ 1–0 line, we utilized the eXtended CASA Line Analysis Software Suite (XCLASS⁷, Möller et al. 2017) by solving the one-dimensional (1D) radiative transfer equation. Here, we assume an isothermal source under the assumption of local thermodynamic equilibrium (LTE). To model and generate synthetic spectra for fitting spectral lines extracted toward the cores and for creating a velocity field map, the myXCLASSFit and myXCLASSMapFit functions were used, respectively. These functions compute a synthetic spectrum for each pixel or an input spectrum by fitting an observed spectrum with Gaussian profiles, assuring that opacity effects on the line shapes are appropriately taken into account. The optimization package MAGIX (Möller et al. 2013) was used to perform the fit. This software acts as an interface between the input codes and an iterating engine. The aim is to minimize the deviation (χ² values) between the modeled outputs and the observational data. An embedded SQLite database was used to acquire molecular properties, such as Einstein A coefficients and partition functions. This database contains entries from the Cologne Database for Molecular Spectroscopy (CDMS, Müller et al. 2005) in its Virtual Atomic and Molecular Data Center (VAMDC, Endres et al. 2016) implementation. Additionally, an extended set of partition function calculations is also available. The initial parameter set for fitting each emission line component comprises the excitation temperature (T_ex), the total column density (N_tot), the line width (Δυ), and the velocity offset (υ_off) from the systemic velocity. The column density is related to the optical depth, τ, by $$N=\frac{8\pi v^3}{c^3}\frac{Q\left(T_{\text{ex}}\right)}{{𝑔}_u{A}_{ul}}\frac{\exp \left(E_u/k{T}_{\text{ex}}\right)}{\exp \left(hv/k{T}_{\text{ex}}\right)-1}\tau \mathrm{\hspace{0.17em}d}v {\text{cm}}^{-2},$$(5)

where E_u and g_u are the energy of the upper state for the selected transition and its degeneracy, respectively; ν is the rest frequency of a selected molecular transition, c is the speed of light, and h is Planck constant. The quantity A_ul is the Einstein coefficient for spontaneous emission, and Q(T_ex) and k are the rotational partition function and the Boltzmann constant, respectively. For the excitation temperature of the N₂H⁺ components, we set a prior with a minimum T_ex of 5 K and a maximum T_ex of 30 K. To acquire the best-fit output parameters by minimizing the χ² value, we used the Levenberg-Marquardt (LM) algorithm. In the fitting procedure, all the hfs transitions (F₁ = 1 – 1, F₁ = 2 – 1, and F₁ = 0 – 1) of N₂H⁺ are considered, and thus the output parameters of the fitting are the hfs deconvolved values. To improve the signal-to-noise level, we separately fitted spectral lines extracted from the areas of cores. Thus, we fixed the core size to the core diameter to fit the individual spectra. Figures 9 and C.1 display the XCLASS modeled spectra and the input spectra, extracted over the core size, in a main beam temperature scale (T_mb). The average values of the decomposed line widths and the column density of N₂H⁺ are 1 km s⁻¹ and 2.4 × 10¹³ cm⁻², respectively. In addition, we also estimate N(N2H⁺) with a fixed temperature using the THCN/HNC averaged over an individual core size. The derived column densities are consistent with the results obtained with T_ex(N₂H⁺) and, thus, we use the N(N₂H⁺) values obtained with T_ex(N₂H⁺).

With the fitted line parameters and the astrodendro output results, we can examine the physical properties of the identified cores and derive the core mass as: $$M_{\text{core }}={\mu }_{{\text{H}}_2}{m}_{\text{H}}\pi N\left({\text{H}}_2\right)R_{\text{core }}^2{\text{M}}_{\odot },$$(6)

by integrating N(H₂) over the core size, where ${\mu }_{{\text{H}}_2}=2.8$ is the mean molecular weight, and m_H is the mass of a hydrogen atom. Here, N(H₂) is the H₂ column density and it is determined by adopting the relative N₂H⁺ abundance to be 4 × 10⁻¹⁰ to H₂ because this is the average value between the abundances of 3.0 × 10⁻¹⁰ in low-mass molecular cloud cores (Caselli et al. 2002) and 5.2 × 10⁻¹⁰ in high-mass molecular cloud cores (Pirogov et al. 2003). This abundance assumption is reasonable at the given N₂H⁺ core sizes with the average radius of 0.025 pc, based on the 1D prestellar core model presented in Fig. 5 of Redaelli et al. (2019) showing X(N₂H⁺) between 10⁻¹⁰ and 10⁻⁹ for core radius of 0.02–0.03 pc. The median and maximum H₂ column densities are 5.27 × 10²² cm⁻² and 1.48 × 10²³ cm⁻², respectively. The mean and median values of M_core are 4.10 M_⊙ and 1.63 M_⊙. With the determined core masses, to address whether the cores are gravitationally bound or unbound, we calculate the virial mass (M_vir) and the virial parameter (α_vir) using the N₂H⁺ line widths. The virial mass of an equivalent homogeneous density sphere (Caselli et al. 2002) is derived as $$M_{\text{vir}}=210\left(\frac{R_{\text{core }}}{\text{pc}}\right){\left(\frac{\text{Δ}v}{ \text{km} {\text{s}}^{-1}}\right)}^2{\text{M}}_{\odot }.$$(7)

Thus, the mean and median values of M_vir are 5.27 M_⊙ and 3.68M_⊙. The virial parameter (Bertoldi & McKee 1992) is defined as the ratio of the virial mass and the actual core mass, $${\alpha }_{\text{vir}}=\frac{M_{\text{vir}}}{M_{\text{core}}},$$(8)

and can be used to evaluate the stability of a core against gravitational collapse. This simplified measurement assumes the core is isolated and does not interchange mass, momentum, or energy with its surrounding medium. Kauffmann et al. (2013) shows that α_vir ≲ 2 means that a core is gravitationally bound, while if α_vir ≳ 2, a core is unbound as random turbulent gas motions create significant support against collapse. In addition, α_vir = 1 indicates that a core is in virial equilibrium, and for the case with α_vir ~ 2, a core is marginally gravitationally bound. Based on the determined virial parameters, 50 % (9/18) of the cores are unbound, while the other half of the cores are bound. This virial parameter could change with the abundance of N₂H⁺, which can vary in different physical environments. If the N₂H⁺ abundance were higher, ~4 × 10⁻⁹ (as found toward star-forming clumps; e.g., Sanhueza et al. 2012) than the assumption (~4 × 10⁻¹⁰), we would obtain smaller core mass as N(H₂) becomes smaller (~5 × 10²¹ cm⁻²). Consequently the virial parameters of all the cores are derived to be much larger, 6–69, indicating a gravitationally unbound nature. According to the 1D prestellar core model presented by Redaelli et al. (2019), the N₂H⁺ abundances are lower, about 10⁻¹⁰, in the inner core regions where a core radius is smaller than 0.03 pc; while the abundance peaks to 10⁻⁸ at a core radius of 0.05 pc. Thus, the virial parameters obtained with ~4 × 10⁻¹⁰ are reasonable as the N₂H⁺ abundances toward the clumps (e.g., Sanhueza et al. 2012) were derived from larger physical scales compared to the cores analyzed in this study. We will discuss the distribution of these bound and unbound cores in connection with stellar feedback in Sect. 5.

Fig. 8 N2H+ velocity-integrated intensity over all the hfs transitions. The leaf structures identified by the Dendrogram algorithm are shown in black contours, and schematic structures of the extracted cores computed from the intensity-weighted second moment (e.g., variance) of intensities in the plane of the sky are marked as blue ellipses. The white contours are 3σrms (0.146 Jy beam−1 km s−1) of the intensity map integrated over a velocity range of 18 km s−1 . All the identified cores are labeled the same as in Table 3. The synthesized beam size for N2H+ emission is displayed in the lower left corner.

Fig. 9 N2H+ spectral lines (black curves) extracted over the circular regions of the compact cores, overlapped with N2H+ spectral lines (red dotted curves) extracted over the leaf structures from the astrodendro and the XCLASS modeled spectra (green curves) in a Tmb scale.

4.6 Gas kinematics toward DR18

To investigate the structures and kinematics of the molecular gas, we generated velocity-channel maps of the HCO⁺ line for the velocity range from 6.5 km s⁻¹ to 14.5 km s⁻¹, with a channel spacing, Δυ, of 0.8 km s⁻¹. As shown in Fig. 10, we detected both compact and extended emission distributions. The peak velocity of the H¹³CO⁺ spectrum over the observed entire area peaks at the systemic velocity, 8.4 km s⁻¹, of DR18 and this velocity component gas is widespread in the Cygnus-X region (e.g., Schneider et al. 2023, 2016b). At velocities from 7.5–9.5 km s⁻¹, close to the systemic velocity, the compact substructures appear, and they show clear emission counterparts in the mm continuum emission and N₂H⁺ integrated emission maps. The emission morphologies at velocities ranging from 6.5 km s⁻¹ to 8.5 km s⁻¹ spatially overlap with the cometary head of the radio continuum emission, indicated by the magenta solid arc lines. The redshifted emission components, which are visible at velocities between 10.5 and 11.5 km s⁻¹, are consistently aligned with the extended mm and cm continuum emission and more concave against the inner cometary HII region, as indicated by the dotted magenta curves. In addition, there are several bright emission knots that prevail even at higher redshifted velocities, reaching up to 14.5 km s⁻¹.

In the channel maps from 10.5 km s⁻¹ to 14.5 km s⁻¹, we notice spatially elongated features as indicated by blue and red ellipses. Unlike other red-shifted features, the locations of these emission components appear to be devoid of cm radio continuum emission; only the northern one (blue ellipse) shows the weak mm continuum. In particular, the southeastern feature clearly shows bullet-like shapes, as are often found in protostellar jets (e.g., Tafalla et al. 2017). The jet-like red-shifted features seem to be originating from the region where C2 and C3 are located, but not from YSO counterparts, which were identified by Comerón et al. (2022). On the other hand, their possible blue-shifted counterparts might be the features elongated toward the northwest direction marked with blue ellipses at velocities of 10.5 km s⁻¹–12.5 km s⁻¹. Their bent shapes could be due to an interaction with the ambient gas components surrounding the HII region. A detailed discussion of these features is beyond the scope of this work.

Apart from the jet-like features, we also identified knotlike features with red-shifted velocities in a range of 11.5– 14.5 km s⁻¹, which are distributed along the northeast edges of the extended radio continuum emission. These knot-like features are certainly compact but not associated with any mm continuum emission features or N₂H⁺ emission. Also, these features are not bright compared with the N₂H⁺ cores in the HCO⁺ intensity map (Fig. 5). To investigate the knots seen in the HCO⁺ emission, we extracted HCO⁺ and H¹³CO⁺ spectra over a twice larger area than the average beam size at the locations of six conspicuous knot-like features selected by visual inspection on the channel maps. In the 12.5 km s⁻¹ channel map, the blue crosses indicate the locations of the extracted HCO⁺ and H¹³CO⁺ spectra. The spectra shown in Fig. 11 correspond to the positions of the markers, numbered 1 to 6, with decreasing declination. Most of the knots display double-peak profiles, except for Knot 1, which only shows the red-shifted component. The lower velocity component is at the systemic velocity of 8.4 km s⁻¹ (shown by the blue dashed line), while the other component is red-shifted (approximately 12.5 km s⁻¹) relative to the systemic velocity. We also compared the HCO⁺ spectra with H¹³CO⁺ spectra shown in green (scaled up by a factor of 5). The gray filled area indicates the rms noise level (not scaled up). H¹³CO⁺ emission is only detected toward knot 5, at the systemic velocity. On the other hand, the other positions are not detected at our noise level. Thus, the double-peak profiles of the main-isotopologue spectra cannot be attributed to self-absorption as no optical thin emission is seen at their peak velocities. We observe two kinematic components toward these features, likely associated with the extended free- free radiation detected in the mm and cm continuum emission based on the morphologies of the continuum emission as well as the systemic gas.

5 Discussion

5.1 Ionized gas and molecular gas emission in DR18

In the 3.6 mm and 6cm radio continuum emission maps (Fig. 3), the morphology of the continuum emission around the DR18- 05 star shows a cometary head pointing toward the dense gas regions. On the other hand, the extended continuum emission has a cometary shape heading toward the Cyg OB2 association, which is the opposite direction of the arc-like dense molecular gas distribution. To investigate whether the continuum emission at 6 cm and 3.6 mm wavelengths is primarily due to DR18-05 or whether there is an additional contribution from the external OB2 association, we estimated its Lyman continuum photon flux (N_i), namely, the number of ionizing photons emitted by a star, to define the spectral type of an ionizing source, with the following equation (e.g., Carpenter et al. 1990), $$\left(\frac{N_i}{{\text{ photons }}^{-1}}\right)=9\times {10}^{43}\left(\frac{{\displaystyle \int S_v}}{\text{mJy}}\right){\left(\frac{D}{\text{kpc}}\right)}^2{\left(\frac{v}{5\text{\hspace{0.17em}GHz}}\right)}^{0.1}.$$(9)

Here, ∫S_ν is the integrated radio flux density estimated at a frequency of ν, while D is the distance to DR18. The equation assumes that the HII region is optically thin. This is consistent with our spectral index analysis discussed above and is to be expected for the DR18 HII region as it has an extended morphology and is not deeply embedded in molecular clouds, for instance, ultra-compact or hyper-compact HII regions (see, e.g., Kurtz et al. 1994).

Considering the flux over the entire 6 cm radio continuum emission obtained from the GLOSTAR, we obtain a Lyman photon continuum flux of 2.4×10⁴⁹ s⁻¹, corresponding to the spectral type of an 05.5 star by following the relation between Lyman continuum photons and zero-age main sequence (ZAMS) spectral type presented in Panagia (1973). Hence, the estimated Lyman continuum photon rate is nearly five orders of magnitude larger than that of a typical B2-type star (i.e., 4.5×10⁴⁴ s⁻¹).

However, if we only consider the continuum flux within the interior HII region exhibited as blue-green color in the three-color composite image in the left panel of Fig. 12, we obtained ZAMS spectral types of B0 – B0.5 (6.9×10⁴⁶ s⁻¹), which are reasonably consistent with the one (B2 type) from the optical spectral analysis carried out by Comerón et al. (2022). In addition, using the 3.6 mm continuum emission toward the HII region by assuming that the mm continuum emission from the HII region is mainly free-free emission, we obtained a very similar Lyman photon continuum flux (i.e., 7.4 × 10⁴⁶ s⁻¹). The ionized gas along the ridge cannot be formed by DR18-05 because the molecular gas seems to block the UV radiation emitted from the B2 star as the PDR layer traced by the PAH 8 µm and H₂ (1–0) S(1) emission lines are located between the ridge and the DR18-05 HII region. Thus, the ionization of the gas on the ridge must have another origin. This suggests that DR18-05 cannot itself excite the entire volume of ionized gas associated with the globule. Thus, the ionization by the external OB stars of the OB2 association likely contributes to the free-free emission observed at 6 cm and 3.6 mm to a significant extent, as discovered in another globule in the Cygnus-X region (Schneider et al. 2021).

In the left panel of Fig. 12, the three-color composite image of HCO⁺, 8 µm, and 6 cm radio continuum emission shows that strong radio continuum emission is extended toward the west within the eroded cavity delineated by the PAH and molecular emission. This feature is visible in the flux difference map of Brγ and H₂ (1–0) S (1) emission lines, as shown in the right panel of Figure 12. The 2.12 µm H₂ (1–0) S(1) emission line is emitted in a cascade after absorbing UV photons (10 eV < hν < 13.6 eV) and it is thus an excellent tracer of the photodissociation interfaces, where FUV radiation (6 eV < hν < 13.6 eV) is absorbed in the cold molecular cloud (Black & van Dishoeck 1987). In addition, the 2.17 µm Brγ emission line is a recombination line of hydrogen coming from an HII region or ionized plasma region (e.g., Djupvik et al. 2017; Carlsten & Hartigan 2018). In the difference map, negative values indicate strong H₂ emission, implying the FUV absorbing surface of the molecular cloud, while positive values indicate an ionized gas region. In the image, as expected, the bright Brγ emission closely resembles the 6 cm radio continuum emission (blue color in the three-color image).

Except for the ionized gas originating from the HII region of DR18-05, we notice widespread ionized gas emission along the rim-brightened edges of DR18, which stretches toward the north of the HII region cavity. The intensities of 6 cm radio continuum, 8 µm, and the H₂ (1–0) S(1) emission are lower toward the rim-brightened ridge evident in the three-color image, compared to values close to the HII region. However, the distributions of 6 cm, 8 µm, H₂ (1–0) S (1), and HCO⁺ toward the ridge show the well-aligned paralleled geometric structure of an ionization front (traced by the 6 cm continuum), a photodissociation layer (in the 8 µm emission), and a molecular gas layer (indicated by the HCO⁺ emission), from west to east. As shown in Fig. 1, the surface of DR18 is irradiated directly by external O/B star members of the OB2 association, ionizing the low-density gas on its surface. This suggests that the ionized gas in the ridge is likely photoevaporating, which is consistent with the analysis with 3.6 mm and 6 cm continuum emission. Such photoevaporating ionized gas surfaces are found in intense radiation regions in the Cygnus-X complex (Emig et al. 2022; Carlsten & Hartigan 2018) and the M16 region (Sofue 2020). Such cases showing both an internal PDR caused by embedded stars and external PDRs due to OB associations have been found in infrared spectroscopic observations of Brγ and H₂ (1–0) S(1) lines (e.g., Djupvik et al. 2017) and PDR tracers, such as the 158 µm [CII] line (e.g., Schneider et al. 2012).

To investigate whether ionizing radiation from the HII region is eroding the globule and affecting the surrounding ambient molecular gas, we must estimate its physical properties. First, we estimated the propagation velocity, υexp., of the eroding front, using the following equation for its propagation velocity (Comerón et al. 2005), assuming that the ionized gas from the molecular cloud is free to flow away from the cloud: $$v_{\text{exp}}=\frac{N_{\text{i}}}{4\pi \left(n_0\right)r^2} \text{km} {\text{s}}^{-1},$$(10)

where n₀ is $2n_{{\text{H}}_2}$, the initial molecular gas density. This assumption is acceptable as the ionized gas is traced by the 6 cm radio continuum, and the Brγ line emission is extended outward from DR18-05. The derived range of υ_exp is −19.8–24.2 km s⁻¹, by adopting N_i of 6.9×10⁴⁶ s⁻¹ obtained from the 6 cm radio continuum emission from the volume of the shell, a molecular gas density of 2 × 10³ cm⁻³ toward DR18, derived from Herschel far- infrared observations from 70 µm to 500 µm (Schneider et al. 2016a), and a distance from DR18-05 of ~0.24–0.27 pc. As mentioned by Comerón et al. (2005), the velocities derived from the above formula stand as upper limits because the ionizing flux impinging on the molecular gas is diminished by ionization of the region between the molecular gas and the ionizing star and also because the real distance from the ionizing star DR18- 05 can be larger than the projected distance. The propagation velocities determined here are typical velocities found toward expanding HII regions from radio recombination lines (e.g., Kim et al. 2017). The structure of the radio continuum emission displays a cometary head pointing toward the molecular gas regions and a tail extending toward the OB2 association. This could be explained by the motion of the DR18-05 star and the expansion of ionized gas in a champagne flow (Churchwell 1999). This explanation is also consistent with the results of Comerón & Torra (1999). However, we also notice a bow shock-like feature around the cometary head. Massive stars such as O/B type stars develop strong stellar winds, injecting their mechanical energy into the ambient gas. Thus, winds can play a significant role in the dynamical evolution of an ionized gas region along with the uniform expansion and the moving star (e.g., hybrid model; Cyganowski et al. 2003; Veena et al. 2017). The shell radius is described by the following equation (Castor et al. 1975), $$R_{\text{sh}}=0.042{\left(\frac{L_{\text{w}}}{{10}^{36}\text{ergs} {\text{s}}^{-1}}\right)}^{1/5}{\left(\frac{n_0}{{10}^5 {\text{cm}}^{-3}}\right)}^{-1/5}{\left(\frac{t}{{10}^3\text{yr}}\right)}^{3/5}\text{pc},$$(11)

where L_w is the mechanical luminosity of the wind. The shell expansion velocity is $$V_{\text{sh}}=24.7{\left(\frac{L_{\text{w}}}{{10}^{36}\text{ergs} {\text{s}}^{-1}}\right)}^{1/5}{\left(\frac{n_0}{{10}^5 {\text{cm}}^{-3}}\right)}^{-1/5}{\left(\frac{t}{{10}^3\text{yr}}\right)}^{-2/5} \text{km} {\text{s}}^{-1}.$$(12)

We adopted a distance of 0.2 pc between the ionizing star and the shell, a velocity of 5 km s⁻¹ for the shell velocity due to the stellar wind, and an ambient gas density of 2 × 10³ cm⁻³ (Schneider et al. 2016a). The shell velocity is defined using the HCO⁺ emission (Fig. 10). The determined mechanical luminosity and the shell expansion time are 3.9 × 10³³ ergs s⁻¹ and 2.4 × 10⁴ yr, respectively. Since DR18 is nearly in the plane of the sky, we also estimated a case with a shell velocity of 10 km s⁻¹, and obtained L_w of 3.1 × 10³⁴ ergs s⁻¹ and a shell expansion time of 1.2×10⁴ yr. To investigate whether the wind is dynamically more important for forming the shell than the expansion driven by the difference in pressure between the ionized and surrounding gas, we compare the relation between the mechanical luminosity, N_i, and n₀ (Shull 1980; Garay et al. 1994), for the case that the wind is more important than the expansion. For this to apply, the following relation must be satisfied, $$\left(\frac{L_{\text{w}}}{{10}^{36}\text{ergs} {\text{s}}^{-1}}\right)>0.33{\left(\frac{N_{\text{i}}}{{10}^{49} {\text{s}}^{-1}}\right)}^{2/3}{\left(\frac{n_0}{{10}^5 {\text{cm}}^{-3}}\right)}^{-1/3}.$$(13)

The derived wind mechanical luminosity (for a case with a shell velocity of 5 km s⁻¹, 3.9 × 10³³ erg s⁻¹, and, for a case with 10km s⁻¹, 3.1 × 10³⁴ erg s⁻¹) is at best comparable to (and possibly only a small fraction of) the luminosity of the ionized gas, 4.4 × 10³⁴ ergs s⁻¹. This implies that the stellar wind alone cannot produce the shell. Hence, it is likely that the cavity of DR18 is shaped by the ionized gas produced by DR18-05, with a minor contribution from its stellar wind and also by the influence of the OB2 association. The cavity is open toward the OB2 association and is photoionized by its combined UV radiation. In addition, along with DR18-05, other YSOs located in DR18 could contribute to forming the shell as well. Located around the shell, there are four compact N₂H⁺ cores (C1, C4, C7, and C9) close to the HII region and, interestingly, these cores are gravitationally unbound, with αvir ≳ 2.0. This may imply that they are affected by the expansion of the HII region toward ambient gas associated with these unbound cores. The majority of the gravitationally bound cores showing αvir ≲ 2.0 (6 out of 8 cores; C6, C8, C10, C11, C13, and C15) have compact NH₃ core counterparts (Zhang et al. 2024), which were observed with a comparable spatial resolution (~0.02 pc) with to CASCADE observations (~0.02–0.03 pc). In addition, these bound cores partially overlap with the 8 µm and H₂ emission, as depicted in Fig. 12. This suggests that these bound cores are located in cold molecular gas regions behind the PDR and the core fragmentation might be influenced by the HII region.

Fig. 10 Velocity channel maps of HCO+ from 6.5 km s−1 until 14.5 km s−1. The velocity difference between the channel maps is 1 km s−1. Each channel map has different color scales and shows emission stronger than a 10 σrms. The flux scale per velocity channel map is varied to better highlight the individual emission features. The black and white contours are the mm and cm continuum emission, respectively, and the contour levels are the same as in Fig. 3. The pink sold lines and dashed lines at velocities from 6.5 km s−1 to 10.5 km s−1 show the outlines of the cavity created by the HII region of DR18-05. The blue ellipses marked at velocities from 10.5 km s−1 to 12.5 km s−1 and the red ellipses marked at velocities of 11.5 km s−1–14.5 km s−1 are indicated to blue and red-shifted lobes of the potential outflow (or jet) candidate, respectively. In the 12.5 km s−1 channel map, the blue cross symbols mark the positions of the extracted HCO+ and H13CO+ spectral lines (shown in Fig. 11) toward compact features.

Fig. 11 HCO+ (black) and H13CO+ (green) spectral lines extracted over an area with a diameter of 2 × θbeam toward pronounced bright knot features marked on the velocity channel map at 12.5 km s−1 in Fig. 10. The blue dashed-vertical lines indicate the systemic velocity of DR18, and the red dashed lines have a velocity of 12.5 km s−1. The gray filled areas indicate ±1σrms (= 7.15 mJy beam−1) for HCO+.

Fig. 12 Different emission distributions toward DR18. Left: three-color composite image of DR18; HCO+ emission in red, 8 µm emission in green, and the GLOSTAR 6 cm radio continuum emission in blue. The 6 cm emission represents a combination VLA D- and B-configuration data from the GLOSTAR survey. The synthesized beam size of the GLOSTAR data is 4″, comparable to that of NOEMA data (~3–4″). The pixel size, of the Spitzer/IRAC 8 µm data is 0.″6. The black dashed lines indicate the cuts for intensity profiles of different emission tracers, presented in Fig. 15. The white rectangles represent the areas of the East and North shells, as well as the Ridge. Right: flux difference image of Brγ and H2 1–0 S (1) emission lines obtained from Comerón et al. (2022). The black and yellow contours show the 5σ level of N2H+ velocity-integrated intensity and a 3.0σ level of SiO velocity-integrated intensity, respectively. In both images, the bright green star indicates DR18-05, the ionizing source in DR18.

5.2 Association between SiO emission and the internal PDR

As mentioned in Section 4.2, the SiO emission has a distinct spatial distribution compared with that of other species and is concentrated between the cold molecular gas and ionized gas regions (see Fig. 5). Figure 13 shows the SiO spectral lines extracted from the NOEMA+30 m combined data (black curves) and NOEMA-only data (red curves). The NOEMA-only data has a better velocity resolution of 0.25 km s⁻¹ than the NOEMA+30 m combined data used for this study, but we only use this data for presenting the SiO spectra. These spectral lines are extremely narrow and weak (~0.8–2.0 km s⁻¹ and a peak intensity ~0.05 Jy beam⁻¹).

In previous studies, SiO emission has been mainly associated with shocks in outflows of embedded YSOs (e.g., Schilke et al. 1997; Gusdorf et al. 2008; Sánchez-Monge et al. 2013; Duarte-Cabral et al. 2014; Csengeri et al. 2016). Lines from SiO associated with outflows typically exhibit broad velocity components >8 km s⁻¹ and its abundance in outflows relative to H₂, X(SiO) is approximately 10⁻⁶ – 10⁻⁹ (e.g., Martin-Pintado et al. 1992; Sánchez-Monge et al. 2013; Duarte-Cabral et al. 2014; Csengeri et al. 2016). On the other hand, SiO detected in this study shows extremely narrow line widths υ_FWHM ~ 0.8–2.0 km s⁻¹, which are much narrower than the typical line widths originating from outflows. Nevertheless, if the axes of the outflows are on the plane of the sky, the Doppler width of outflow components can be quite modest. Thus, we compared the distributions of SiO emission in DR18 with that of the known YSOs identified by Comerón et al. (2022), as displayed in the SiO velocity field map in the upper panel of Fig. 14. However, there is only one Class I YSO (indicated by a black square marker in Fig. 14), which is likely to drive an outflow (e.g., Bontemps et al. 1996), located at the gap between the east shell and the north shell of the HII region. We can distinguish separate blue- and red-shifted emission regions, none of which appear to be associated with any YSO; this would contradict an outflow-based interpretation. Toward the east shell region, the SiO emission is red-shifted toward the HII region, whereas the SiO emission, toward the N₂H⁺, is blue-shifted. This gradient is the opposite of the velocity gradient seen in N₂H⁺ (bottom panel of Fig. 14), showing that the further components from the HII region show red-shifted velocities. In the north shell region, N₂H⁺ and SiO emission lines show more complicated velocity gradient features. In addition, we determined N (SiO) over areas twice bigger than the SiO beam size, using 1.8 × 10¹² ʃ T_mb dυ cm⁻², with a T_ex of 10 K, taken from Csengeri et al. (2016). SiO column densities toward different positions along the east and north shells are ~1.8 – 2.0 × 10¹² cm⁻². Using H₂ column densities (1.8–3.8 × 10²² cm⁻²) obtained from Bonne et al. (2023), the derived X(SiO) values are ~5 × 10⁻¹¹–1.0 × 10⁻¹⁰. The abundances are much lower than the typical abundances from outflows and closer to the values (~10⁻¹⁰–10⁻¹¹) found toward low-velocity shocked gas regions experiencing global infall, converging flows, or cloud-cloud collisions toward star-forming molecular clouds (e.g., Jiménez-Serra et al. 2010; Nguyen-Lu’o’ng et al. 2013; Duarte-Cabral et al. 2014; Csengeri et al. 2016; Cosentino et al. 2018, 2019, 2020; Kim et al. 2023), along with (~10⁻¹¹) determined toward the Orion bar region and the S 140 PDR (i.e., Schilke et al. 2001). Thus, this SiO emission is unlikely related to any potential outflows associated with the Class I YSO source.

Intriguingly, the SiO emission is significantly confined to narrow regions and spatially overlaps well with strong vibrational H₂ line emission, evident in the negative values in the flux difference map in the right image of Fig. 12. The H₂ (1–0) S (1) emission line is enhanced in the vicinity of the HII region where SiO emission is also detected. The excellent spatial agreement with the enhanced H₂ v = 1 − 0 S (1), narrow line widths, and lower abundances of SiO emission suggest that both trace PDR material around the HII region. H₂ lines can be excited by absorbing FUV radiation, followed by fluorescence (Black & van Dishoeck 1987). Finally, SiO molecules can enter the gas phase via the photo-desorption of the icy mantles of dust grains (Walmsley et al. 1999; Schilke et al. 2001).

However, it is hard to discard a shock origin for SiO and H₂ 1−0 S (1) emission because H₂ can also be excited by collisions due to shock heating of H₂ molecules (Burton et al. 1989). Thus, to investigate whether SiO and H₂ are associated with PDRs or shocked regions, we compared their emission with infrared emission in the 4.5 µm and 3.6 µm Spitzer/IRAC bands. Toward star-forming regions, especially infrared dark clouds, 4.5 µm emission is found to be enhanced and extended, known as green fuzzies⁸ (Chambers et al. 2009) or extended green objects (EGOs) (Cyganowski et al. 2008). This enhancement is likely due to shock excited line emission, either from the H₂ 0–0 S (9) line at 4.69 µm or CO v = 1 - 0 rovibrational bandhead at 4.3–5.2 µm (Noriega-Crespo et al. 2004; Ray et al. 2023). Consequently, green fuzzies mostly mark outflows (e.g., Beuther et al. 2005; Cyganowski et al. 2011). In addition, unlike other IRAC bands, this band does not cover PAH lines and is considered a “PAH-free” band. On the other hand, the 3.6 µm band includes 3.3 µm features from neutral PAHs, while the 5.8 µm and 8µm bands include emission from cationic PAHs (PAH⁺) (Bakes et al. 2001a,b; Benjamin et al. 2003). Comparing the emission in the 4.5 µm band with that from these other bands can help to determine the contributions of H₂ excited by shocks and by UV radiation, which is at the origin of the fluorescent IR emission from PAHs. It may also help elucidate the origin of the enhanced 2.12 µm H₂ emission and the narrow SiO emission. However, the [4.5 µm]/[3.6 µm] ratios toward the regions with the SiO detections and strong H₂ emission are remarkably lower ([4.5 µm]/[3.6 µm] < 0.65) than the color criteria value ([4.5 µm]/[3.6 µm] ≥ 1.8) for EGOs (Chambers et al. 2009) (see Fig. D.1). Such low ratios might be due to either a lack of only a very small contribution from shock-excited H₂ and/or CO lines to the 4.5 µm band or strong PAH features at 3.6 µm band, or both. The study of van den Ancker et al. (2000b,a) reported that H₂ lines of embedded sources are due to shocks, while toward more evolved objects (e.g., the shell of DR18), such lines mainly originate from a PDR as the envelopes disperse; meanwhile, the PDR components become dominant compared to the shock-excited emission.

Figure 15 shows normalized intensity profiles of N₂H⁺, SiO, 8 µm, and 6 cm radio continuum emission along different cuts (black dashed lines in Fig. 12) toward the east and north shells. In many cases, the SiO (indicated by pink curves) is enhanced at offset positions with enhanced 2.12 µm H₂ emission, evident by negative values on the Brγ − 2.12 µm H₂ map (blue curves), toward the east and north shell regions. Overall, these two tracers show substantial correlation over most of the six cuts of those shells. Interestingly, the north and middle cutouts of the east shell show that SiO emission is spatially detected between peaks of N₂H⁺ emission (black curves), a dense gas tracer, and 8 µm emission (bright green curves), which traces PDRs. Toward the middle and south cutouts in the east shell, we recognize that some peaks of SiO and N₂H⁺ are situated close to 6 cm continuum emission peaks tracing the position of DR18-05. These N₂H⁺ peaks correspond to the gravitationally unbound cores (i.e., C9 and C4) and one bound core (i.e., C5). For the SiO peaks in the south cutout, the SiO emission is slightly red-shifted (υ_offset > 8.8 km s⁻¹, while the N₂H⁺ emission is relatively blue- shifted, ~υ_offset < 8.4 km s⁻¹. Thus, these two tracers do not trace the same molecular gas regions. All the spatial distributions and intensity profiles for different tracers indicate that the SiO emission toward DR18 seems conspicuously associated with the PDR surrounding the HII region of DR18-05.

To investigate the origin of SiO, we analyzed irradiated low- velocity shock models (project ID: shockgrid_110_q1m8_2022) using the Paris-Durham shock code (version 1.1.0, revision 115) obtained from the InterStellar Medium DataBase (ISMDB)⁹. The detailed concept of irradiated shock models using the Paris- Durham shock code is described in Godard et al. (2019) and the descriptions for this grid of shock models used here are presented and explained in Kristensen et al. (2023). The grid of shock models covers low-velocity shocks propagating through molecular gas environments with various external UV radiation intensities. Table 4 lists the fixed and free parameters used in the Paris-Durham low-velocity shock models we studied in this work. We examined shock velocities (v_S) ranging from 5 km s⁻¹ to 20 km s⁻¹, initial density (n_H = n(H) + 2n(Hշ)) in the range of 10² cm⁻³ to 10⁵ cm⁻³, and UV radiation field intensities (G₀) from 0.1 to 10³ in Mathis unit (Mathis et al. 1983). However, the magnetic parameter, b is fixed to 1, corresponding to a C-type shock. The initial transverse magnetic field and the magnetic parameter are expressed as $\frac{B_0}{1\text{\hspace{0.17em}μG}}=b\times {\left(\frac{n_{\text{H}}}{1 {\text{cm}}^{-3}}\right)}^{0.5}$.

In addition, the fractional abundance of PAHs relative to n_H is fixed to 10⁻⁸. In particular, the initial fractional elemental abundances (n_x/n_H) of silicon (Si) are 3 × 10⁻⁶ in the gas-phase, as Si⁺, and 3.37 × 10⁻⁵ in grain cores. This grid only considers gasphase chemistry and does not include grain-grain interactions or ice mantles on the grains (see Kristensen et al. 2023 for details). However, grain erosion resulting in the release of elemental Si, Fe, and others into the gas phase is taken into account. Thus, in cases with low-velocity shocks, the final SiO abundance greatly depends on the initial abundance of Si⁺ and formation routes of silicon-contained molecules in different physical conditions.

Figure 16 presents SiO and N₂H⁺ abundances calculated across irradiated shocks propagating in different media, as a function of radiation field intensities. For all the models, X(SiO) and X(N₂H⁺) decrease as radiation field intensities increase. In these models, the maximum values of the fractional abundances relative to H₂, N(X)/N(H₂) are approximately close to the initial fractional abundance (n_x/n_H) of the gas-phase Si (i.e., Si⁺). Overall, the computed SiO abundances for all the input shock velocities in high-density environments are several orders of magnitude greater than those in low-density cases. In the middle and right panels, for shock velocities of 10 km s⁻¹ and 20 km s⁻¹, the abundances of X(SiO) can still be produced at higher levels in the UV-irradiated environments (G₀ ≳ 100) than the observed values if the initial densities are ~10³–10⁵ cm⁻³. According to Schneider et al. (2016a), the measured UV radiation field in the region is about G₀ ~ 10³. However, the initial external UV radiation might be lower than G₀ ~ 10³ determined from the Herschel 70 µm and 160 µm fluxes because the DR18 region may be less exposed to the OB association and the external UV radiation field in the HII region may be smaller. Considering that SiO emission is only detected toward the molecular shells where dense molecular cores are closely located, it is likely that the initial density is higher than 10³ cm⁻³. Therefore, the shock velocity creating the SiO emission in DR18 is likely slower than 10 km s⁻¹.

In the left panels of Fig. 16, the 5 km s⁻¹ shock model irradiated by UV radiation of G₀ = 10³ and with an initial density of 10⁵ cm⁻³ demonstrates a strong agreement with the observed SiO abundances and the physical conditions of the DR18 region. Godard et al. (2019) found that only low-velocity C-type shocks (v_s ≤ 5 km s⁻¹) can exist in typical PDRs (n_H ~ 10⁴ cm⁻³ and G₀ ~ 100). Low-velocity shocks in highly UV-illuminated regions (G₀ = 10⁴) only cause a slight compression of gas, which increases its thermal pressure and changes its interaction with the external UV field (Godard et al. 2019). The DR18 region likely falls between typical PDRs and the extremely irradiated regions. The SiO emission observed in DR18 is evidently associated with gas regions that have been compressed by low-velocity shocks (~5 km s⁻¹ ) propagating within an initial medium density ranging between 10⁴ and 10⁵ cm⁻³ . These shocks are further influenced by external UV radiation (G₀ ~ 100–10³). In addition, the X(N₂H⁺) values for all the irradiated shock models are significantly lower than the expected values, ~10⁻¹⁰–10⁻⁹ and undetectable. This possibly explains spatial and velocity discrepancies between SiO and N₂H⁺ toward the east and north shells. By comparing the different shock models, we suggest that the SiO emission in the DR18 region originates from the compressed gas layers resulting from low-velocity shocks, which are likely produced by the initial expansion of the HII region (e.g., Liu et al. 2024) and stellar winds, irradiated by external UV radiation.

Fig. 13 SiO spectral lines extracted from the NOEMA+30 m combined data (black curves) and the NOEAM-alone data (red curves). The combined data has a channel width of 0.8 km s−1 while the NOEMA-only data provides a narrower channel width of 0.25 km s−1. The coordinates of each extracted position on the SiO intensity map are individually displayed on each panel.

Fig. 14 Velocity field maps of SiO (top) and N2H+ (bottom). The black contours are the 3.6 mm continuum emission. The N2H+ velocity field map is generated by fitting hyperfine structure lines of N2H+ per pixel using XCLASS while the SiO velocity map is the intensity-weighted velocity of the spectral line (i.e., moment 1 map). The black symbols are YSOs, and they are the same as presented in Fig. 2. The velocity ranges in the color bars for both maps are fixed as the same.

Fig. 15 Intensity profiles of different tracers (see legend in the top left panel) along different cuts through the eastern shell (left column) and the northern shell (right column). All the molecular line intensities were extracted from their velocity-integrated intensity maps. The extents and directions of the cuts are shown in the left three-color image of Fig. 12. For improved visibility, the Brγ-H2 intensity profiles are shifted by −0.8, while the SiO profiles are moved by +0.2 along the y-axis. The vertical dashed lines indicate the depths in the intensity profiles of Brγ-H2.

Fig. 16 Fractional abundances, N (X)/N (H2), of SiO and N2H+ relative to H2, calculated across shocks propagating at 5 (left panels), 10 (middle panels), and 20 (right panels) km s−1 in the initial media density nH = 102 (red squares), 103 (green circles), 104 (cyan diamonds), and 105 (blue triangles) cm−3, as function of irradiating radiation field G0 in Mathis unit. In the upper panels for SiO, the filled gray areas indicate the observed SiO abundances (5 × 10−1–1 × 10−10) toward DR18.

6 Summary

We carried out 3.6 mm wavelength observations with the IRAM 30 m telescope and NOEMA toward DR18 to analyze the continuum radiation emission and selected spectral emission lines from various molecular species. This study addresses the complex molecular gas substructures inside DR18 for the first time.

• We detected ground- or near ground state emission lines from HCO⁺, HCN, HNC, C₂H, H₂CO, SiO, N₂H⁺, and HC₃N, along with the rare isotopologues H¹³CO⁺, H¹³CN, and HN¹³C, and ¹³CS. In addition, we also detected three deuterated species, DCN, DNC, and NH₂D. HCO⁺, HCN, HNC, C₂H, and H₂CO show extended molecular gas distributions, whereas the emission from isotopologues and the three deuterated molecules arise from dense, compact regions. Interestingly, SiO emission shows a spatial distribution different from that of any other species.

• The 3.6 mm continuum emission comprises a significant contribution from free-free radiation, as extrapolated from the 6 cm radio continuum emission. The range of spectral types (B0 – B0.5) of a star emitting the free-free radiation from the continuum emission approximately matches that of the B2 type star DR18-05. The cometary HII region located in a cavity of the DR18 globule is likely formed by a champagne flow. Additionally, the shell structure seen in the cm continuum emission surrounding this HII region is probably the result of expansion driven by the HII region, as well as photoionization from DR18-05 and externally coming from stars in the OB2 association. The spatially extended ionized gas encompassing the entire DR18 globule is attributed to photoevaporation from the external OB stars in the Cyg OB2 association.

• Using the velocity-integrated intensity map of N₂H⁺, we identified 18 compact cores that are all located around the HII region of DR18-05. The average size determined via source extraction is 0.03 pc. These cores are relatively colder (T_HCN/HNC < 30 K) than their surrounding regions (T_HCN/HNC > 30 K). Half of the cores (9 out of 18) are gravitationally bound and most of them are located behind PDRs surrounding the HII region of DR18-05. This might suggest that these cores form as a result of stellar feedback that has triggered star formation.

• We found emission from SiO with narrow line widths (~0.8–2.0 km s⁻¹) and low abundance (~5 × 10⁻¹¹ –1 × 10⁻¹⁰). This SiO emission is not spatially coincident with that of other dense gas tracers, and its velocities show that it is not associated with the N₂H⁺ cores that trace the densest gas. On the other hand, the SiO emission shows an excellent spatial agreement with PDRs located between the N₂H⁺ cores and the HII region of DR18-05. The PDRs also show stronger H₂ emission. In addition, the ratios of emission of the [4.5]/[3.6] IRAC bands toward the SiO emitting regions do not show any evidence of shocks, as the values are much smaller ([4.5 µm]/[3.6 µm] < 0.65) than typical ratios found toward shocked regions ([4.5 µm]/[3.6 µm] ≥ 1.8). In comparison to the irradiated shock models, we suggest that the SiO emission encompassing the HII region indicates that the molecular gas regions are slightly compressed by low- velocity shocks (with ~5 km s⁻¹) irradiated by external UV radiation (with intensities ranging from 100 to 1000), as they propagate through the medium with n_н ~ 10⁴ to 10⁵ cm⁻³. These shocks are presumably generated by the initial HII region expansion and possibly by stellar winds as well.

Acknowledgements

This work is based on observations made with the Institut de Radioastronomie Millimétrique (IRAM) 30 m telescope and the Northern Extended Millimeter Array (NOEMA). W.-J. K. was supported by DLR/Verbundforschung Astronomie und Astrophysik Grant 50 OR 2007. A.S.-M. acknowledges support from the RyC2021-032892-I grant funded by MCIN/AEI/10.13039/501100011033 and by the European Union ‘Next Genera- tionEU’/PRTR, as well as the program Unidad de Excelencia María de Maeztu CEX2020-001058-M, and support from the PID2020-117710GB-I00 (MCI-AEI- FEDER, UE). S.A.D. acknowledges the M2FINDERS project from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant No 101018682). D. S. acknowledges support from the European Research Council under the Horizon 2020 Framework Program via the ERC Advanced Grant Origins 83 24 28 (PI: Th. Henning). W.-J. K., N. S, and P. S. acknowledge support by the Deutsche Forschungsgemeinschaft via the collaborative research center SFB 1601 (project ID 500700252), subprojects A2 and B2.

Appendix A Velocity-integrated intensity maps

Fig. A.1 Velocity-integrated intensity maps of CCH (integrated from 6km s−1 to 14km s−1), H2CO (6 km s−1 to 11.6 km s−1), HCN (−0.4 km s−1to 18.8 km s−1, including all the hfs transitions), HNC (7km s−1 to 13.5 km s−1), H13CN (0km s−1 to 16.4 km s−1 including all the hfs transitions), and HN 13C (6.8 km s−1 to 10.8 km s−1) from top to bottom. The star symbol and contours are the same as in Fig. 5.

Fig. A.2 Velocity-integrated intensity maps of H13CO+ (6.8 km s−1 to 11.6 km s−1), HC3N (J = 10 − 9) (7.6 km s−1 to 10.8 km s−1), HC3N (J = 8 − 7) (7.6 km s−1 to 10.8 km s−1), DCN (5 km s−1 to 17 km s−1 spanning all the hfs transitions), DNC (7 km s−1 to 10km s−1), and 13CS (7 km s−1 to 10 km s−1 ) from top to bottom. The star symbol and contours are the same as in Fig. 5.

Appendix B T_HCN/HNC and $T_{H^{13}CN/H{N}_{13}C}$ maps

Fig. B.1 Gas kinetic temperature map separately derived by using the ratios of the HCN/HNC and H13CN/HN13C emission lines. The black contours represent the 3.6 mm continuum emission. The beams for the THCN/HNC and $T_{{\text{H}}^{13}\text{CN}/{\text{HN}}^{13}\text{C}}$ maps are displayed in the lower left corner.

Appendix C N₂H⁺ spectral lines extracted toward the identified compact cores

Fig. C.1 N2H+ spectral lines (black curves) extracted over the circular regions of the compact cores, overlapped with N2H+ spectral lines (red dotted curves) extracted over the leaf structures from the astrodendro and the XCLASS modeled spectra (green curves) in a Tmb scale

Appendix D Flux ratio map of Spitzer/IRAC 4.5 µm and 3.6 µm bands

Fig. D.1 Ratio map of 4.5 µm over 3.6 µm fluxes. The magenta and black contours for SiO and N2H+, respectively, are the same as the right image of Fig. 12.

References

All Tables

All Figures

Fig. 1 Spitzer/IRAC 8 µm emission map (in a flux range of 0–98 MJy/sr) of a part of the Cygnus-X complex (Beerer et al. 2010). The yellow rectangle outlines DR18, which is studied here. The white circle has a radius 12.2 pc (at a distance of 1.4 kpc), with its center at the position of the trapezium (Cyg OB2 #8, RA 20h33m 16s and Dec +41°18′45″) of O stars. The cyan and bright green markers are known members of O- and B-type stars, respectively, in the OB2 association (Wright et al. 2015).

In the text

Fig. 2 Three-color composite images of DR18 (8 µm in red, 4.5 µm in green, and 3.6 µm in blue) with different maximum flux levels for three Spitzer/IRAC bands to show the extended emission in a larger FoV (left panel) and the emission details, which are saturated in the left panel, in a zoom-in view (right panel) covering our 5.1′ × 5.1′ target region (outlined by a white square box in the large FoV image). In both images, cyan contour lines represent the combined VLA B+D array 6 cm radio continuum emission from the GLOSTAR survey by levels of 0.35, 0.70, 1.05, 1.40, 1.75, 2.10, 2.45, 2.80, 3.15, and 3.50 mJy/beam. In the zoom-in image, the different color symbols indicate young stellar objects identified by infrared colors (Comerón et al. 2022): yellow squares mark Class I YSOs, light blue diamonds for Class II YSOs, red upside-down triangle for a flat spectrum YSO, purple circles for Class III YSOs, and bright green upside triangles represent transition-disk YSOs. The yellow half-arc line indicates the open shell or the ear-like feature in the infrared emission image.

In the text

Fig. 3 Millimeter (left) and centimeter (right) continuum emission maps. The mm continuum emission map represents the combination of the GBT MUSTANG2 and NOEMA data sets, whereas the cm continuum emission map is generated by combining the VLA D and B configuration observations. The FWHM synthetic beams of these maps are displayed in the lower left corners. The contour levels of the mm continuum emission are 0.46, 0.96, 1.38, 1.80, 4.14, 5.33, and 6.51 mJy beam–1, whereas for the cm emission, the levels are 0.35, 0.70, 1.05, 1.40, 1.75, 2.10, 2.45, 2.80, 3.15, and 3.50 mJy beam–1. The star symbol indicates the position of the B2 star.

In the text

Fig. 4 90 GHz – 5.8 GHz spectral index map in the left panel and its spectral index error map in the right panel. The higher frequency data are smoothed to the beam size of the lower frequency data to achieve a common resolution of 4″. The star symbol indicates the position of the B2 star. The black contours are the 3.6 mm continuum flux, and their levels are the same as in the left panel of Fig. 3, and the blue contours are the 6 cm continuum emission with the same contour levels shown in the right panel.

In the text

Fig. 5 Velocity-integrated intensity maps of HCO+ (6km s−1 to 15.6km s−1), N2H+ (−1.2km s−1 to 17.2km s−1 spanning all the hfs transitions) in the upper row and SiO (8km s−1 to 10km s−1) and NH2D (0km s−1 to 15.2km s−1 covering all the hfs transitions) in the lower row. The black contours represent the GBT+NOEMA combined 3.6 mm continuum emission, and the bright green marker indicates the position of the ionizing source, DR18-05. The contour levels of the mm continuum are the same as shown in Fig. 3. The horizontal white scale bar has a size of 0.1 pc at a given distance of 1.4 kpc. The beam of the maps of each species is displayed on the left corner of each image panel.

In the text

	Fig. 6 Intensity ratio map of SiO/N2H+ in color scale for the overlapping areas of both molecular species. The black and pink contours represent the 5σ level of N2H+ velocity-integrated intensity and a 3σ level of SiO velocity-integrated intensity, respectively. The bright green star indicates DR18-05, the ionizing source in DR18.
In the text

Fig. 7 Dust and gas temperature maps. Top: dust temperature map obtained from Bonne et al. (2023) generated with the 160 µm/250µm flux density ratio as described in the work by Palmeirim et al. (2013). Bottom: gas kinetic temperature map combined temperature maps derived from the HCN/HNC and H13CN/HN13C emission line ratios. The black contours represent the 3.6 mm continuum emission. The bright green star indicates the position of the B2-type star D18-05. The beams for the Tdust and THCN/HNC maps are 18″ and 3.43 × 2.76 with a position angle of 13.66° displayed in the lower left corner.

In the text

Fig. 8 N2H+ velocity-integrated intensity over all the hfs transitions. The leaf structures identified by the Dendrogram algorithm are shown in black contours, and schematic structures of the extracted cores computed from the intensity-weighted second moment (e.g., variance) of intensities in the plane of the sky are marked as blue ellipses. The white contours are 3σrms (0.146 Jy beam−1 km s−1) of the intensity map integrated over a velocity range of 18 km s−1 . All the identified cores are labeled the same as in Table 3. The synthesized beam size for N2H+ emission is displayed in the lower left corner.

In the text

	Fig. 9 N2H+ spectral lines (black curves) extracted over the circular regions of the compact cores, overlapped with N2H+ spectral lines (red dotted curves) extracted over the leaf structures from the astrodendro and the XCLASS modeled spectra (green curves) in a Tmb scale.
In the text

Fig. 10 Velocity channel maps of HCO+ from 6.5 km s−1 until 14.5 km s−1. The velocity difference between the channel maps is 1 km s−1. Each channel map has different color scales and shows emission stronger than a 10 σrms. The flux scale per velocity channel map is varied to better highlight the individual emission features. The black and white contours are the mm and cm continuum emission, respectively, and the contour levels are the same as in Fig. 3. The pink sold lines and dashed lines at velocities from 6.5 km s−1 to 10.5 km s−1 show the outlines of the cavity created by the HII region of DR18-05. The blue ellipses marked at velocities from 10.5 km s−1 to 12.5 km s−1 and the red ellipses marked at velocities of 11.5 km s−1–14.5 km s−1 are indicated to blue and red-shifted lobes of the potential outflow (or jet) candidate, respectively. In the 12.5 km s−1 channel map, the blue cross symbols mark the positions of the extracted HCO+ and H13CO+ spectral lines (shown in Fig. 11) toward compact features.

In the text

Fig. 11 HCO+ (black) and H13CO+ (green) spectral lines extracted over an area with a diameter of 2 × θbeam toward pronounced bright knot features marked on the velocity channel map at 12.5 km s−1 in Fig. 10. The blue dashed-vertical lines indicate the systemic velocity of DR18, and the red dashed lines have a velocity of 12.5 km s−1. The gray filled areas indicate ±1σrms (= 7.15 mJy beam−1) for HCO+.

In the text

Fig. 12 Different emission distributions toward DR18. Left: three-color composite image of DR18; HCO+ emission in red, 8 µm emission in green, and the GLOSTAR 6 cm radio continuum emission in blue. The 6 cm emission represents a combination VLA D- and B-configuration data from the GLOSTAR survey. The synthesized beam size of the GLOSTAR data is 4″, comparable to that of NOEMA data (~3–4″). The pixel size, of the Spitzer/IRAC 8 µm data is 0.″6. The black dashed lines indicate the cuts for intensity profiles of different emission tracers, presented in Fig. 15. The white rectangles represent the areas of the East and North shells, as well as the Ridge. Right: flux difference image of Brγ and H2 1–0 S (1) emission lines obtained from Comerón et al. (2022). The black and yellow contours show the 5σ level of N2H+ velocity-integrated intensity and a 3.0σ level of SiO velocity-integrated intensity, respectively. In both images, the bright green star indicates DR18-05, the ionizing source in DR18.

In the text

	Fig. 13 SiO spectral lines extracted from the NOEMA+30 m combined data (black curves) and the NOEAM-alone data (red curves). The combined data has a channel width of 0.8 km s−1 while the NOEMA-only data provides a narrower channel width of 0.25 km s−1. The coordinates of each extracted position on the SiO intensity map are individually displayed on each panel.
In the text

Fig. 14 Velocity field maps of SiO (top) and N2H+ (bottom). The black contours are the 3.6 mm continuum emission. The N2H+ velocity field map is generated by fitting hyperfine structure lines of N2H+ per pixel using XCLASS while the SiO velocity map is the intensity-weighted velocity of the spectral line (i.e., moment 1 map). The black symbols are YSOs, and they are the same as presented in Fig. 2. The velocity ranges in the color bars for both maps are fixed as the same.

In the text

Fig. 15 Intensity profiles of different tracers (see legend in the top left panel) along different cuts through the eastern shell (left column) and the northern shell (right column). All the molecular line intensities were extracted from their velocity-integrated intensity maps. The extents and directions of the cuts are shown in the left three-color image of Fig. 12. For improved visibility, the Brγ-H2 intensity profiles are shifted by −0.8, while the SiO profiles are moved by +0.2 along the y-axis. The vertical dashed lines indicate the depths in the intensity profiles of Brγ-H2.

In the text

Fig. 16 Fractional abundances, N (X)/N (H2), of SiO and N2H+ relative to H2, calculated across shocks propagating at 5 (left panels), 10 (middle panels), and 20 (right panels) km s−1 in the initial media density nH = 102 (red squares), 103 (green circles), 104 (cyan diamonds), and 105 (blue triangles) cm−3, as function of irradiating radiation field G0 in Mathis unit. In the upper panels for SiO, the filled gray areas indicate the observed SiO abundances (5 × 10−1–1 × 10−10) toward DR18.

In the text

Fig. A.1 Velocity-integrated intensity maps of CCH (integrated from 6km s−1 to 14km s−1), H2CO (6 km s−1 to 11.6 km s−1), HCN (−0.4 km s−1to 18.8 km s−1, including all the hfs transitions), HNC (7km s−1 to 13.5 km s−1), H13CN (0km s−1 to 16.4 km s−1 including all the hfs transitions), and HN 13C (6.8 km s−1 to 10.8 km s−1) from top to bottom. The star symbol and contours are the same as in Fig. 5.

In the text

	Fig. A.2 Velocity-integrated intensity maps of H13CO+ (6.8 km s−1 to 11.6 km s−1), HC3N (J = 10 − 9) (7.6 km s−1 to 10.8 km s−1), HC3N (J = 8 − 7) (7.6 km s−1 to 10.8 km s−1), DCN (5 km s−1 to 17 km s−1 spanning all the hfs transitions), DNC (7 km s−1 to 10km s−1), and 13CS (7 km s−1 to 10 km s−1 ) from top to bottom. The star symbol and contours are the same as in Fig. 5.
In the text

	Fig. B.1 Gas kinetic temperature map separately derived by using the ratios of the HCN/HNC and H13CN/HN13C emission lines. The black contours represent the 3.6 mm continuum emission. The beams for the THCN/HNC and $T_{{\text{H}}^{13}\text{CN}/{\text{HN}}^{13}\text{C}}$ maps are displayed in the lower left corner.
In the text

	Fig. C.1 N2H+ spectral lines (black curves) extracted over the circular regions of the compact cores, overlapped with N2H+ spectral lines (red dotted curves) extracted over the leaf structures from the astrodendro and the XCLASS modeled spectra (green curves) in a Tmb scale
In the text

	Fig. D.1 Ratio map of 4.5 µm over 3.6 µm fluxes. The magenta and black contours for SiO and N2H+, respectively, are the same as the right image of Fig. 12.
In the text