© The Authors 2025

1 Introduction

The dynamic evolution of star formation leads to the production of a diverse suite of molecules. Molecular emission thus reveals details about the physical and chemical conditions from which these species arise. In particular, wide-band millimetre and sub-millimetre observations are capable of capturing numerous rotational transition lines which are powerful diagnostics of their environment due to the different excitation conditions each line is sensitive to.

In early-stage protostellar sources in particular, when the accreting protostar is still enshrouded in an envelope, two categories of molecules are distinguished for their diagnostic capabilities: complex organic molecules (COMs; carbon-bearing saturated molecules with ≥6 atoms) and carbon-chain molecules (CCMs; unsaturated hydrocarbons). These species are sensitive to their environment and are valuable tracers alongside more common, simple species such as CO, HCN, or HCO⁺.

Interstellar COMs are abundant in the hot (T > 100 K) and dense (n_H2 > 10⁶ cm⁻³) cores of star-forming regions (Blake et al. 1987; Caselli et al. 1993; Cazaux et al. 2003; Ceccarelli 2004; see Herbst & van Dishoeck 2009 for a review of these molecules).

On the other hand, CCMs are known to exist in cold (T ≤ 10 K) molecular clouds and starless cores (Avery et al. 1976; Broten et al. 1978; Kroto et al. 1978; Little et al. 1978; Sakai et al. 2010). They are now also known to form in the warm gas enveloping protostellar systems in a process dubbed warm carbon chain chemistry (WCCC; Sakai et al. 2008; Aikawa et al. 2008; Sakai & Yamamoto 2013).

Due to this chemical differentiation, COMs and CCMs can be used to classify protostellar systems. While these species have been used to study star-forming regions, especially hot cores, high-mass star-forming regions (HMSFRs), which are capable of producing stars with M > 8 M_⊙, are more difficult to study due to their relative rarity, clustered environments, and quick evolutionary timelines. As the evolutionary processes of HMSFRs are not fully understood, observing their molecular emission is a valuable tool for studying these regions. The first HMSFR found with evidence of WCCC was DR21(OH) (Mookerjea et al. 2012), a nearby and prominent region within Cygnus X.

DR21 (OH) lies 2′ north of the DR21 H II region at a distance of 1.5 kpc (Rygl et al. 2012). Its peak, known as DR21(OH)M or N44, has a total mass of about 10³ M_⊙ and a density n_H2 of 10⁶ cm⁻³; it has been resolved into two sources, MM1 and MM2 (Woody et al. 1989; Mangum et al. 1991; Motte et al. 2007). MM1 has the stronger continuum emission at millimetre wavelengths (Liechti & Walmsley 1997). Both MM1 and MM2 are thought to be very young, with no visible H II region, and only very weak continuum emission at centimetre wavelengths (Argon et al. 2000). They are thought to be in a pre-ultracompact H II region phase. Dust continuum observations suggest that MM1 is the brighter source (L = 1.7 × 10⁴ L_⊙; consistent with a B0V star) and shows evidence of star formation, whereas MM2 (L = 1.2 × 10³ L_⊙; early B star), though more massive, is fainter and most likely at an even earlier stage of evolution. These sources have been further differentiated: Zapata et al. (2012) found nine millimetre sources between MM1 (SMA 5–9) and MM2 (SMA 1–4), each with masses in the range 8–24 M_⊙; Minh et al. (2012) identified two hot sub-cores in millimetre observations of MM1 corresponding to SMA 6 and 7. The presence of methanol masers (Class I, pumped collisionally) along a bipolar outflow is further evidence that DR21(OH)M is a region of active star formation (Plambeck & Menten 1990; Kogan & Slysh 1998; Kurtz et al. 2004; Araya et al. 2009).

Near DR21(OH)M are several other dense clumps: DR21(OH)S (or N48), DR21(OH)W (N38), DR21(OH)N1 (N41), DR21(OH)N2 (N40), and N36 (Mangum et al. 1991; Chandler et al. 1993; Motte et al. 2007). These sources exist along the DR21 ridge, which consists of filamentary structures that influence collapse (Schneider et al. 2010; Hennemann et al. 2012). N48, N38, and N40 are massive IR-quiet cores noted as possible precursors to high-mass protostars; the lack of IR emission suggests that they are not yet forming a star, but their SiO emission is indicative that an outflow is being powered by a protostellar object (Motte et al. 2007). Mangum et al. (1991) found N48 and N38 to be high-luminosity objects, consistent with B-type stars, but they did not detect any maser emission, which suggests that these sources are in the early stages of star formation. In total, these sources make up a parsec-scale region with a mass of 4900 M_⊙ (based on 1.2 mm continuum emission; Schneider et al. 2010).

Interferometric observations have catalogued DR21(OH)M as a rich molecular source that comprises hydrocarbons and COMs (Zapata et al. 2012; Minh et al. 2012; Girart et al. 2013; Orozco-Aguilera et al. 2018). Observations of the larger-scale environment detail the continuum emission and that of simple species such as CO, CS, SiO, N₂H⁺, and HCO⁺(Richardson et al. 1994; Lai et al. 2003; Motte et al. 2007; Schneider et al. 2010). We aim to link the molecular environment of DR21(OH)M at core scales to the larger clump and filament it resides in. We used single-dish observations to map CH₃CCH as a CCM tracer as well as CH₃OH and H₂CO as COM tracers, and used a local thermodynamic equilibrium (LTE) model to determine the physical parameters of different gaseous components within the DR21 (OH) region. The rest of the paper is structured as follows: Sect. 2 presents the observations and the dataset used, and Sect. 3 the LTE modelling results and analysis. In Sect. 4 we present a chemical evolution model and discuss the findings, and in Sect. 5 provide a brief summary.

2 Observational data

2.1 Observations

We observed DR21 (OH) in multiple frequency ranges with the Institut de radioastronomie millimétrique (IRAM) 30 m telescope and the 100 m Green Bank Telescope (GBT). A detailed summary of these observations is found in Freeman et al. (2023).

The IRAM 30-m observations were completed November 2020 and April 2021 (project codes 021-20 and 122-20, PIs S. Bottinelli and R. Plume). We mapped a 1^′ × 1.5^′ region in DR 21(OH) in the on-the-fly observational mode with position switching. We observed the frequency ranges 131.2−138.98 GHz and 287.22−295 GHz simultaneously with the Eight MIxer Receiver (EMIR). The offset position was −120^′′ horizontally and −240^′′ vertically. The map centre is αJ(2000)=20^h39^m01.00 and δJ(2000)= +42^∘22^′48^′′.0. Data reduction, and the production of maps, was completed with GILDAS/CLASS¹. At 291 GHz, the spectral channel width is 781 kHz or 0.80 km s⁻¹, the beam size is 9.3^′′, the pixel size is 4.4^′′, and the average noise level is 79.7 mK. At 135 GHz, the spectral channel width is 391 kHz or 0.88 km s⁻¹, the beam size is 20.5^′′, and the pixel size is 9.7^′′. Given the larger beam size, the 135 GHz data will not be discussed in this paper.

The GBT observations were completed in March 2021, in the frequency ranges 84.5−85.75 GHz and 95.55−96.8 GHz using all 16 beams of the Argus focal plane array and the VErsatile GBT Astronomical Spectrometer (VEGAS) spectral line backend (project code 21A-039, PI P. Freeman). We obtained 1′ DAISY on-the-fly maps for both frequency ranges, centred as above, using position-switching for the reference measurements. The spectral sampling with spectrometer mode 2 is 92 kHz in both frequency ranges, or 0.32 km s⁻¹ at 85 GHz and 0.28 km s⁻¹ at 96 GHz. The beam size is 9.2^′′ at 96 GHz and 10.0^′′ at 85 GHz, both with a pixel size of 2.0^′′. The average noise level is 24.0 mK at 85 GHz and 17.7 mK at 96 GHz. The data were reduced and calibrated using GBTIDL². In Sect. 3, the data are presented in units of T_MB. We converted from $T_{\mathrm{A}}^{\ast }$ using the telescope B_eff/F_eff values provided by the respective observatories: GBT at 85 GHz, 0.4545; GBT at 96 GHz, 0.3838; IRAM EMIR at 290 GHz, 0.547.

2.2 Species of study

In these datasets we detect numerous lines of CH₃CCH, a CCM tracer, and of CH₃OH and H₂CO, to be used as COM tracers. The lines were identified in CASSIS³ (Vastel et al. 2015) using the Cologne Database for Molecular Spectroscopy⁴ (CDMS) catalogue (Müller et al. 2005).

Propyne, or methyl acetylene (CH₃CCH), is a symmetric rotor with many transitions closely spaced in frequency. Numerous lines are thus captured in a small bandwidth, and are useful for analysis of the physical conditions, especially temperature, in LTE (Irvine et al. 1981; Askne et al. 1984; Kuiper et al. 1984; Bergin et al. 1994). It has long been detected in low- and highmass star-forming regions and is used as a tracer of chemical complexity (Snyder & Buhl 1973; Lovas et al. 1976; Cazaux et al. 2003; Taniguchi et al. 2018; Giannetti et al. 2017; Santos et al. 2022).

Methanol, CH₃OH, is the simplest alcohol molecule and an asymmetric top molecule (Ball et al. 1970). It has numerous transitions in the millimetre and sub-millimetre ranges, which are abundant and commonly used as probes of the physical environment. H₂CO was the first polyatomic organic molecule detected in space (Snyder et al. 1969), and while not technically a complex species it is chemically associated with larger organic molecules. Both H₂CO and CH₃OH are formed by the successive hydrogenation of CO on dust grain surfaces, and are precursor molecules to larger COMs (Charnley et al. 1997; Garrod & Widicus Weaver 2013). These two species are detected across interstellar space: in cold clouds, hot cores, outflows, shocks, the Galactic centre, and external galaxies (see Herbst & van Dishoeck 2009, and references therein).

Table 1 lists the transition lines seen in our frequency ranges, with certain conditions based on the sensitivity of our observations. We limited the scope of lines: for CH₃CCH an A_ij > 1 × 10⁻⁶ s⁻¹ and E_up < 200 K (to account for the a-CH₃CCH line that has E_up < 150 K, but has a higher tabulated value when the a- and e-types are combined in CDMS); for H₂CO an A_ij > 1 × 10⁻⁶ s⁻¹ and E_up < 150 K; for CH₃OH an A_ij > 2 × 10⁻⁶ s⁻¹ (to remove the 84.5 GHz maser line) and E_up < 150 K. Integrated intensity maps of select lines covering a range of upper energy levels are displayed in Fig. 1.

Class I masers are often found in young star-forming regions, including in DR21(OH) (Ladeyschikov et al. 2019; Slysh et al. 1997; Araya et al. 2009). In particular, the CH₃OH line at 84.5 GHz is a known maser in DR21(OH) (Batrla & Menten 1988). These masers, which are collisionally pumped, are unlikely to be well reproduced by LTE. Thus, to be cautious, we did not include this line in our LTE analysis.

3 LTE modelling: Results and analysis

An LTE model was used to determine the physical parameters of the observed region – size, line width, column density, excitation temperature, and source velocity. In LTE, we assume the excitation temperature represents the gas temperature. We used a Python-based LTE model based on the CASSIS formalism⁵ (previously described in Freeman et al. 2023 for a similar analysis of AFGL 2591 and IRAS 20126+4104). The LTE model identifies all molecular transitions for one or more species and simultaneously fits these lines defined by the physical parameters. The resulting fit is the combination of parameters that best reproduces the spectral line profile as determined by a least-squares minimisation (represented as a reduced χ² value). In the model, we used the CDMS tags⁶ for each species (Müller et al. 2005). The species types – a- and e-CH₃CCH, A- and E-CH₃OH, and o- and p-H₂CO – were not differentiated, as we did not detect enough lines of each.

In the simplest case, an observed spectral line will arise from one physical gas component modelled by a single Gaussian profile with a systemic velocity (v_LSR), line width (Δ v_FWHM), and brightness (T_A^*). However, in the case where the telescope’s beam encompasses multiple physical components, the observed spectra may have to be modelled using multiple Gaussian profiles.

We treated each species, CH₃CCH, CH₃OH, and H₂CO, separately in our analysis. To determine the components of each species, we started with a one-component fit and looped the model over all pixels in the map. When a one-component fit was not adequate, as determined via either a visual inspection of the spectral line profile or the resulting χ² analysis, we added components as needed. This divided the map into regions described by their own, distinct, physical components. The regions, marked in Fig. 2, correspond to the clumps found in Motte et al. (2007). We gained an idea of what components were necessary based on (a) taking an example spectrum and fitting multiple Gaussians in CASSIS to find velocity and line width information, and (b) knowledge from the literature on what velocity and temperature components we expect to see at this scale. The results for these fits are presented in Table 2. To measure the quality of the fits, we present the reduced χ² values in Fig. A.1. This reduced value is the χ² divided by the number of degrees of freedom – the number of data points minus the number of variables – and should be close to 1 for a good fit.

Fig. 1 Integrated intensity maps for select transition lines of CH3CCH, H2CO, and CH3OH. The dense clumps as identified in Motte et al. (2007) are shown as stars, with the main clump N44 in black (see also Fig. 2). The quantum states and upper energy level of the transition are listed in the top left of each plot.

Fig. 2 Regions isolated for LTE modelling near and in DR 21(OH), marked in white boxes over an integrated intensity map of the 12 K CH3CCH 50−40 line. From the top, the regions include: N40, N36, and N41 (labelled N); N44; N38; and N48. The dense clumps, as designated by Motte et al. (2007), are marked by stars.

3.1 N44

All three species required multiple components near DR21 (OH)M, as broad and multi-peaked spectra are seen (Fig. A.2 shows the spectra of the central pixel, marked by the black star in Fig. 1).

We started with two components, N44 1 and N44 2, based on the v_lsr values of the dusty condensations MM1 (about −4.1 km s⁻¹) and MM2 (−0.7 km s⁻¹) provided by Mangum et al. (1992), Mookerjea et al. (2012), and Minh et al. (2012). In all three species, we find components with similar characteristics. Table 2 lists the component parameters and Fig. 3 shows the column density maps (Fig. A.9 shows the error).

We find two cooler and extended components in CH₃OH 3 and 4 and CH₃CCH 1 and 2. We limit the region in Fig. 3 to 18^′′, to encompass the MM1 and MM2 regions where the lines are strong enough to fit well. CH₃OH 3 (22 K) and CH₃CCH 1 (42 K) peak just south-east of MM1 and extend to most of the modelled N44 region. CH₃OH 4 (15 K) and CH₃CCH 2 (38 K) peak in between MM1 and MM2 and are slightly less extended. They all have the same Δ v = 2.8 km s⁻¹. In the interferometric observations of Minh et al. (2012) and single-dish observations of Weaver et al. (2017) there are CH₃OH components of 1820 K. Weaver et al. (2017) also find one component of CH₃CCH at 37 K. At these temperatures, our cooler components likely trace the extended envelope around the MM1 and MM2 dusty condensations.

We find two warmer and more concentrated components in CH₃OH 1 and 2 and H₂CO 1 and 2. CH₃OH 1 (82 K) and H₂CO 1 (81 K) peak on MM1. CH₃OH 2 (75 K) and H₂CO 2 (84 K) are the most compact and are concentrated just north of MM2, showing only a few pixels of prominence. These components are broader than the cool components, with a Δ v = 4.0 km s⁻¹. Comparably, in H₂CO, Zhao et al. (2024) find a kinetic temperature of 103 K for N44. Their spectral line profile clearly shows two velocity components (visually, similar to ours) but in their model they could not distinguish different components. In CH₃OH, Minh et al. (2012) find a 200 K core component and Weaver et al. (2017) find a 93 K component. This suggests that at our spatial scales, on clump scales, we are likely smoothing out the hot gas near the SMA cores, thus the lower temperatures.

For the velocities of each of these components, we find a distinction between those peaking near MM1 or near MM2. Near MM1 the warm component is near −3.6 km s⁻¹ in CH₃OH and −4.1 km s⁻¹ in H₂CO, while the cold component is around −4.1 to −4.2 km s⁻¹. Near MM2, the warm component is at −0.5 km s⁻¹ while the cold component is at −0.9 km s⁻¹. While our v_lsr values (Table 2) do not exactly match previous results, which contain components varying from −5 km s⁻¹ to 1 km s⁻¹ (Mangum et al. 1992; Minh et al. 2012; Girart et al. 2013; Mookerjea et al. 2012), the spatial scales over which we are modelling may provide slightly different average velocities.

A broad feature that is blueshifted relative to the source v_lsr of −4.1 km s⁻¹ is seen in both CH₃OH 5 and H₂CO 3 (Fig. A.3 shows the spectra of a pixel about 6^′′ west of MM1 and MM2, where the outflow is most prominent). We fixed this component to a Δ v of 15 km s⁻¹ and a v_lsr of −8 km s⁻¹. We are not able to reconcile the temperature of these two species in this component. H₂CO is fit well at 80 K, while CH₃OH is at 25 K. This could be due to the energy levels of the lines we detect, as CH₃OH contains much lower energy lines, or, it could be due to nonthermal behaviour.

We first introduced this broad component based on other, similar features seen the H₂CO and CH₃OH maps of Zapata et al. (2012) Fig. 2. Their broad lines emanate from SMA 6 and 7 with the blue- to redshifted features going along an E-W gradient. There is a known maser outflow seen in Vallée & Fiege (2006) or Araya et al. (2009), with a redshifted wing of a low-velocity outflow towards the NW. This also corresponds well to the maser observations of Plambeck & Menten (1990) and Kogan & Slysh (1998); the former suggest that the maser emission could be powered by outflows interacting with dense clumps of gas. However, our blueshifted westward component may not correspond to this particular outflow. Instead, it may correspond to the outflow emission in CO(2−1) and SiO presented in Zapata et al. (2012) that differs from their methanol outflow. These are two classic outflow tracers. The CO(2−1) emission features a high-velocity bipolar flow with its blueshifted wing towards the west, and SiO shows a low-velocity flow with the same spatial pattern. This blueshifted feature is also seen in CO and SiO in Lai et al. (2003), Motte et al. (2007), and Schneider et al. (2010). The broad H₂CO and CH₃OH component we observe, then, is likely tracing the CO outflow and is not associated with the maser outflow.

In the model spectrum (Figs. A.2 and A.3), certain transition lines are underfit in CH₃OH. There is known non-thermal emission in this region, identified mainly through the prominent maser emission (Plambeck & Menten 1990; Kogan & Slysh 1998; Kurtz et al. 2004; Araya et al. 2009). Weaver et al. (2017) found two transition lines at 230 and 250 GHz, which are not known as masers, with non-LTE behaviour and suggest they could be powered by similar mechanisms. Our lines at 290 GHz are similarly above the frequency at which any maser lines have been distinguished, but indicate non-LTE behaviour in our models.

To investigate the possible non-LTE behaviour of CH₃OH, we used a Markov chain Monte Carlo RADEX (van der Tak et al. 2007) code in CASSIS. Unlike the LTE model, RADEX is too computationally intensive to run over all pixels in the map. We chose to model the central pixel of the map as representative of MM1 and MM2 (Fig. A.2) and a pixel 6^′′ (three pixels) west of the central pixel in N44 where the outflow component is seen to be the most prominent (Fig. A.3).

We used the same parameters as in the LTE model for the size, line widths, and velocities. Otherwise, the remaining fitted parameters are in Table 3. We used the VASTEL Database⁷, which uses the A- and E-CH₃OH forms, as there are only collisional coefficients for the separated A- and E- forms. These use p-H₂ as the collision species. The isotopic ratio was set at 1.

The results show CH₃OH is near LTE with most excitation and kinetic temperatures within 1−10 K of each other. In the centre pixel (Fig. A.2) most lines are just as well represented by the LTE model. In the outflow (Fig. A.3), certain lines are much better represented in the RADEX model. There, the GBT triplet at 96.74 GHz (E_up = 20, 7, 13 K), and the IRAM lines at 290.1, 290.3, and 292.7 GHz (E_up = 49, 75, 71, 64 K). The GBT lines, which are dominated by the colder components in LTE, are likely tracing sub-thermally excited gas. These GBT and IRAM lines show visually the most prominent outflow wings, and could have non-thermal excitation due to the shock conditions.

We tried running the LTE model with adjusted temperatures from these RADEX results. In component 2, a higher temperature does not make a clear difference. In components 3 and 4, a higher temperature from the outflow results produces worse results. Similarly, an adjusted temperature in component 5 produces worse results. However, components 3 and 4 peak near the centre pixel, where the LTE and RADEX temperatures align, and component 5 peaks near the outflow pixel, where similarly the LTE and RADEX temperatures align.

For the other clumps, as they show a similar pattern in the CH₃OH models, we discuss the RADEX modelling along with the LTE modelling. The LTE modelling provides an estimate of the gas conditions on a clump scale, while the RADEX modelling provides additional information on the parameters for a single, central pixel.

Fig. 3 Column density maps, as resulting from an LTE model, for all species (rows) in components (columns) of N44. The cores SMA 7 in MM1 and SMA 3 in MM2 are labelled ‘1’ and ‘2’, respectively, and the clump of N44 is marked by a star (as designated in Zapata et al. 2012 and Motte et al. 2007). CH3OH 5 has been moved to align with H2CO 3, the broad components. Values for all LTE model parameters are listed in Table 2.

3.2 N48

N48 is a well-studied clump, with Mangum et al. (1991) showing it also as a high luminosity source with a potentially embedded star and Motte et al. (2007) classifying N48 as a massive IR-quiet protostellar source. In our results, a one-component fit worked well for CH₃CCH and H₂CO for most of the map (Fig. A.4). Figure 4 shows the column density maps for these species, which were found at 34 K and 68 K, respectively (Fig. A.11 shows the error).

In CH₃OH, we identify one component at 28 K and another at 8 K (Table 2). They had equal line widths of 3.2 km s⁻¹ and we fixed the velocity of the second component to that of the first component. The region is known to have structure on smaller scales – with the Plateau de Bure Interferometer Bontemps et al. (2010) found five fragments in 1 mm emission and two main fragments in 3 mm emission, while Csengeri et al. (2011a) found two components in H¹³CN at a v_lsr of −2.39 and −4.77 km s⁻¹. Our two components are at −3.83 and −4.03 km s⁻¹, similar, but again likely offset due to the spatial scales and environments we can trace. In addition, Mangum et al. (1991) see an E-W extension in the source. Our second component of CH₃OH (Fig. 4) has a westward extension, and towards the edge of the modelled region extends towards the north-west, towards N38.

Similar to N44, there are underfit CH₃OH lines in the LTE model. Using RADEX, we fitted two components with very similar T_ex, v_lsr, and full width at half maximum (FWHM) values. There was a larger N_tot in component 1 and a smaller one for component 2, and, overall, this model better represented both the lower energy GBT lines and the higher energy IRAM lines (Fig. A.4, last column). However, with such similar values we conclude that the gas is in near-LTE.

Fig. 4 Column density maps for all species in components of the region N48. The clump location is indicated with a star (Motte et al. 2007).

3.3 N38

N38 is seen as a luminous source with a potentially embedded star and a N–S extension, according to Mangum et al. (1991), with Motte et al. (2007) classifying it as a massive IR-quiet protostellar source. In our data, we identify one component each in CH₃CCH and H₂CO and two components in CH₃OH. Figs. 5, A.10, and A.5 display the column density maps, error maps, and spectra.

Similar to the previous sources, we are seeing warm gas at the scale of the clump. H₂CO peaks at 53 K just west of the Motte et al. (2007) clump peak while the CH₃CCH component peaks just south of the clump at 29 K. This is similar to other temperature derivations of the region, where Mangum et al. (1992) find two components at 26 and 25 K in NH₃ and Zhao et al. (2024) find a H₂CO kinetic temperature of 57 K. We are not necessarily smoothing out a hot core here.

In CH₃OH, we see one warm, 30 K peak that dominates the fit in the IRAM lines and spatially matching the H₂CO peak (Fig. 5). We also find one cold, 8 K, component that dominates the fit in the GBT lines and peaks to the south-east of the warm component. This aligns spatially with the cold CH₃OH component that extends from N48 to the north-west (Fig. 4). We infer this cold gas is coming from the extended filament of the DR21 ridge. Similarly, Csengeri et al. (2011b) found this cold extended gas in between N38 and N48 with single-dish observations of N₂H⁺.

Once again, the LTE model under-fits CH₃OH lines in both the GBT and IRAM ranges, and we are able to reproduce the lines well in RADEX (Fig. A.5). We find in RADEX one warmer, 25 K, component and one colder, 13 K, component. Overall, the parameters are similar to the LTE model, except with the higher N_tot contained in the warmer component 1 compared to the LTE model where it was in the cold component 2.

Fig. 5 Same as Fig. 4 but for the clump N38.

Fig. 6 Same as Fig. 4 but for the N region, containing the clumps N40, N41, and N36.

3.4 N: N40, N41, and N36

Towards the north region (N), we grouped together N40 (Fig. A.6), N41 (Fig. A.7), and N36 (Fig. A.8), as seen in Fig. 2. The column density maps in Fig. 6 (Fig. A.12 show the error), show that the H₂CO is fit well by one component, while CH₃CCH requires two components: one 30 K component with a line width of 2.6 km s⁻¹ in between the three clumps identified by Motte et al. (2007), and a 25 K component with a line width of 1.2 km s⁻¹ towards the outer regions of N41 and N36. In CH3OH we identify three components: one warm, 30 K component and two cold, 6 K components, one moderately broad at 3.7 km s⁻¹ and one narrow at 1.2 km s⁻¹ (see also Table 2). The velocities were fixed to that of the first component.

Motte et al. (2007) classify N40 as a massive IR-quiet protostellar source from 1.2 mm continuum emission, yet Bontemps et al. (2010) note it is weak at 3 mm and may be extended, unlike the other dense clumps they study. Compared to our other clumps, we also see weaker emission in this region (Fig. 1). In N48 and N38 the strength of the CH₃OH lines with E_up = 49 K, < 21 K are comparable (Figs. A.4 and A.5), while in N40 and N41 the 49 K CH₃OH lines are much weaker and the emission is dominated by the cold components.

We fitted each clump separately in RADEX (Figs. A.6, A.7, and A.8). All three clumps were fitted with two warmer components, at ∼30 K and ∼24 K, and one colder component at ∼10 K. As the components are not equally prominent in each clump, we expect to see variation in the RADEX parameters between clumps. In the LTE model, we are capturing the larger-scale trends as we try to find parameters suitable to the whole modelled region. Overall, the RADEX results only varied slightly from the LTE model.

4 Chemical modelling: Results and analysis

The components distinguished by our LTE model show the large-scale environment around many, known core-scale features. DR21(OH), within the larger DR21 ridge, is a young complex with embedded protostars that have not yet disrupted the parent cloud.

In order to determine the possible origin of these molecules in the observed environments, we modelled the species abundances using the NAUTILUS gas-grain code (Ruaud et al. 2016). We did this for all species in all components of DR21(OH), or N44 (we do not use the N44 designation for each component in this section, as all are N44). NAUTILUS uses three phases for chemical evolution – gas, grain surface, and grain mantle and stimulates the abundances over time as a function of different physical conditions – gas density, dust and gas temperature, visual extinction, ultraviolet flux and cosmic ray ionisation. In the code, coupled differential equations describe a network of over 10 000 reactions: 7627 gas-phase reactions and 3799 grain reactions, which involve adsorption onto surfaces, desorption from surfaces, surface reactions, processing within the mantle, and migration between the mantle and the surface. These reactions connect 510 gas-phase species, 225 grain surface species, and 225 grain mantle species.

First, we derived the observed fractional abundances for all components of CH₃CCH, H₂CO, and CH₃OH. We took the ratio of the observed column density for a given component, which we set at the peak value of the LTE model (Sect. 3.1), and the column density of molecular hydrogen, H₂. We found the H₂ column density at the same position by re-gridding and smoothing the Cao et al. (2022) H₂ map in CASA⁸ onto our grid and beam size. Cao et al. (2022) derived their map using archival data described in Cao et al. (2019): Herschel/PACS 160 μm and James Clerk Maxwell Telescope (JCMT)/SCUBA-2 450 and 850 μm images.

Reports of the H₂ to A_V conversion value range from $A_V=\frac{N_{H_2}}{2.1\times {10}^{21}}$ (Rachford et al. 2009; Zhu et al. 2017) to $A_V=\frac{N_{H_2}}{1.9\times {10}^{21}}$ (Bohlin et al. 1978; Whittet 1981). In this study, we assumed an ‘average’ value of $A_V=\frac{N_{{\mathrm{H}}_2}}{2\times {10}^{21}}$. We used a standard cosmic ray ionisation rate of 1.3 × 10⁻¹⁷ s⁻¹. For the UV field, the NAUTILUS code provided a function of the form S × 10⁸ photons cm⁻² s⁻¹ (Ruaud et al. 2016) where S is a scaling factor. We used a scaling factor of 1 to model a general interstellar radiation field. We note that this could be higher in the DR21 environment; however, given the high H₂ column densities and resulting values of A_V, higher UV field strengths are unlikely to have any effect on the chemistry due to shielding.

We used volume densities, n_p−H2, similar to the CH₃OH RADEX models. We matched the other species components as: CH₃OH components 1 and 2 at 1 × 10⁸ to H₂CO components 1 and 2; CH₃OH component 5 at 5 × 10⁶ to H₂CO component 3; CH₃OH components 3 and 4 at 2 × 10⁷ to CH₃CCH components 1 and 2. This comparison is based on the temperature, line width, and velocity parameters found for each species (Table 2). Such high densities are not unusual for a dense core (Mookerjea et al. 2012; Girart et al. 2013) and are representative of (near) LTE conditions.

4.1 Model A: Two-stage warmup

In NAUTILUS, we started with a two-stage model (Model A): a cold quiescent cloud followed by a warmup stage powered by embedded protostars. The cold cloud had initial abundances representing the diffuse interstellar medium (ISM; Table 4), a temperature of 10 K, a density of 10⁴ cm⁻³, a visual extinction of 50 magnitudes, and it evolved for 10⁵ years (Ruaud et al. 2016, as in). Then, in the warmup phase the gas and dust temperature is set to the observed temperature of each pixel found from the LTE model, and the visual extinction is calculated from the average H₂ column density (Table 5). These are constant through the warmup stage. There is a discontinuity from the cold cloud stage that is attributed to a rapid change as a result of the formation of the protostar. This stage evolved for another 10⁶ years.

For all three components of H₂CO, the warmup stage was adequate to reproduce our observed gas-phase abundances (Fig. 7, top). For components 1 and 2, it takes just over 100 years to produce the observed abundances, while for component 3, the assumed outflow, it takes over 10³ years for the abundances to match. The same chemistry applies to all components, on different timescales – at the time the abundances match, the dominant production mechanism changes (Figs. B.1 and B.2). The main production mechanisms at the start are the gas-phase reactions

$$\mathrm{O}+{\mathrm{C}\mathrm{H}}_3\to \mathrm{H}+{\mathrm{H}}_2\mathrm{C}\mathrm{O},$$(1)

and

$$\mathrm{C}\mathrm{H}+{\mathrm{H}}_2\mathrm{O}\to \mathrm{H}+{\mathrm{H}}_2\mathrm{C}\mathrm{O}.$$(2)

Then, these are replaced by desorption of H₂CO off the grain surfaces. For all components, the destruction mechanism follows a similar timeline. H₂CO is, at first, solely destroyed by the gasphase reaction

$${\mathrm{H}}_2\mathrm{C}\mathrm{O}+\mathrm{C}\to {\mathrm{C}\mathrm{H}}_2+\mathrm{C}\mathrm{O},$$(3)

and after 100 or 10³ years, depending on the component, the dominant destruction mechanism becomes adsorption onto the grain surface. At this point, the abundance of gas-phase C has sharply decreased.

Overall, the temperatures of ∼80 K are enough to produce the observed abundances. The experimentally determined binding energy of H₂CO is reported as 3260 K, which corresponds to a desorption temperature of 110 K (as derived in Penteado et al. 2017). At smaller scales, Zhao et al. (2024) report the H₂CO kinetic temperature to be 103 K, so we are likely seeing the hot core smoothed out, and diluting the temperature with the cooler, extended gas. The precursor species of H₂CO on grain surfaces is HCO, which has a lower surface binding energy of 1355 K (Penteado et al. 2017), but is highly reactive and therefore not a barrier to forming H₂CO (Stantcheva et al. 2002).

As CH₃CCH is seen as an envelope molecule, we tried simulating a cloud-edge environment with an A_V of 1 (Fig. 7, bottom). In this scenario, the chemical evolution is influenced by the external interstellar radiation field (e.g. Spezzano et al. 2016; Kalvāns 2021). The abundances for both components are reached a few hundred years after the cold cloud stage. At 10⁴ years, there is a sharp decrease in abundance.

The main production route (Fig. B.3), from about 20 years on, is grain surface hydrogenation:

$$\text{gr}(\mathrm{H})+\text{gr}\left({\mathrm{C}\mathrm{H}}_2\mathrm{C}\mathrm{C}\mathrm{H}\right)\to {\mathrm{C}\mathrm{H}}_3\mathrm{C}\mathrm{C}\mathrm{H}.$$(4)

Prior to this, it is the grain surface reaction:

$$\text{gr}(\mathrm{C}\mathrm{H})+\text{gr}\left({\mathrm{C}}_2{\mathrm{H}}_3\right)\to {\mathrm{C}\mathrm{H}}_3\mathrm{C}\mathrm{C}\mathrm{H}.$$(5)

The grain surface production is largely dominant over gas-phase production (as in Hickson et al. 2016; Andron et al. 2018; Calcutt et al. 2019). CH₃CCH has a binding energy of 4290 K, requiring a temperature of 144 K for thermal desorption (Penteado et al. 2017). The precursor species, CH₂CCH and C₂ H₃, have binding energies of 3840 and 1760 K, respectively (Penteado et al. 2017), which are both lower than for CH₃CCH. Whether or not these species desorb prior to forming CH₃CCH, the production of gas-phase CH₃CCH is dominated by chemical desorption processes. This is not consistent with CH₃CCH being formed in a WCCC scenario. Once the abundance peaks, the main destruction route is photodissociation to H and CH₂CCH – assumed to be influenced by the low A_V environment.

At about 10⁴ years, there is a sharp change in the CH₃CCH formation routes. We determine that, in the low-extinction environment where there is a large H abundance, H atoms are cycling on and off the grain surface. At this abrupt peak, this stops and H is primarily frozen off onto the grains. On the grain surface, H is used up in hydrogenation reactions. Notably, the main reactions during this period are

$$\text{gr}(\mathrm{H})+\text{gr}\left({\mathrm{C}\mathrm{H}}_2\mathrm{C}\mathrm{C}\mathrm{H}\right)\to \text{gr}\left({\mathrm{C}\mathrm{H}}_3\mathrm{C}\mathrm{C}\mathrm{H}\right),$$(6)

where CH₃CCH stays on the grain surface, and:

$$\text{gr}(\mathrm{H})+\text{gr}(\mathrm{C}\mathrm{C}\mathrm{H})\to {\mathrm{C}}_2{\mathrm{H}}_2.$$(7)

Once the grain surface CH₂CCH and CCH populations are diminished, the gas-grain cycling of atomic H returns. This also explains the sharp decrease in gas-phase CH₃CCH at about 10⁴ years in Fig. 7, as the above reaction does not liberate CH₃CCH off the grain surface.

However, DR21(OH) is located within a clustered region, and we do not know if the extent of CH₃CCH is truly on the envelope edge. Values of A_V up to ∼5 still worked in Model A if the UV field in DR21 is larger than unity. However, modelling CH₃CCH with a higher A_V corresponding to cloud core conditions (as in H₂CO, A_V = 50) failed to reproduce the observed abundances at any value of the UV field. In this case, a shock – a sharp spike in temperature and density – was needed to reproduce the CH₃CCH abundances. This three-stage shock model, Model B, was also introduced for CH₃OH, and is described next.

Fig. 7 Model A, the two-stage warmup model, for H2CO (top) and CH3CCH (bottom). The 10 K cold cloud stage (left) is followed by a warmup stage where the LTE components are separated based on the temperatures and densities observed. The time axis represents the age of each component independently. The components of each model are in dotted, dashed, or solid lines; red or green lines are the model abundances in the gas phase, yellow the model abundances in the grain phase, and black the observed abundances calculated from the averages in Table 2.

Fig. 8 Model B, the three-stage shock model, for CH3OH, simulating a shock. The 10 K cold cloud stage (left) is followed by a sharp increase in density and temperature to simulate a shock (middle) and then a post-shock stage (right), when the gas settles to the LTE conditions observed. The time axis represents the age of each component independently. The LTE components are separated based on their temperatures and densities. The components of each model are in dotted, dashed, or solid lines; blue lines are the model abundances in the gas phase, yellow the model abundances in the grain phase, and black the observed abundances calculated from the averages in Table 2. As C1 and C2 are warmer, we suggest they are still in the shock stage, whereas the cooler C3, C4, and C5 are in a post-shock state.

4.2 Model B: Three-stage shock

For CH₃OH, Model A, with any value of A_V or n_p−H2, did not reproduce the abundances. We then used the approach as in Freeman et al. (2023) and simulated a shock passing through the cloud. Shocks, which can be produced by outflows such as the ones in our source, can liberate grain species into the gas phase causing a sharp increase in gas-phase abundances (Herbst & van Dishoeck 2009; Palau et al. 2017). While we cannot model a shock directly, we simulated shock-induced, non-thermal desorption by implementing a three-stage shock model, Model B: a cold quiescent cloud, a ‘shock’ with a sharp increase in temperature and density, and a post-shock stage where the gas settles to the observed conditions.

The cold cloud is the same as in Model A. In the shock stage, the gas temperature and density evolve as in Palau et al. (2017) for IRAS 20126 (see their Fig. 5, top). This assumes a C-type shock as in Jiménez-Serra et al. (2008) with a shock velocity of 40 km s⁻¹ and a pre-shock density of 10⁴ cm⁻³. This lasts for 10⁴ years, during which the temperature sharply increased to 2500 K before decreasing, while the density slowly increased from 10⁴ to 8.2 × 10⁴ cm⁻³, where it stabilised. We fixed the dust temperature to 80 K. Then, in the post-shock phase, we returned the temperature to those observed from the LTE model (as used in the warmup stage above), set the density at 10⁵ cm⁻², and allowed the cloud to evolve for a final 10⁶ years (Table 5). For all five components, this reproduced the abundances (Fig. 8).

In component 2, the warm 100 K region, the main, and for the most part only, production route is desorption from the grain surface (Fig. B.4):

$$\text{gr}\left({\mathrm{C}\mathrm{H}}_3\mathrm{O}\mathrm{H}\right)\to {\mathrm{C}\mathrm{H}}_3\mathrm{O}\mathrm{H}.$$(8)

We assumed component 1 follows the same route as it has similar physical parameters (Table 2). Due to their higher temperature compared to the other components, we assumed components 1 and 2 are still in a shock stage, as shown in Fig. 8, where their abundances match the model after less than 10 years. The binding energy for CH₃OH is 3820 K, corresponding to a desorption temperature of 128 K (Penteado et al. 2017). Similar to H₂CO, the precursor species CH₃O has a lower binding energy (2655 K), but is a reactive intermediate product (Stantcheva et al. 2002). In Model B and with the lower temperatures found in the LTE model, we see that some non-thermal desorption mechanism after grain surface production is key to producing the gas-phase abundances.

For both the cold and narrow component 4 and the cold and broad component 5, chemical desorption reactions dominate the production of gas-phase CH₃OH through the reactions (Figs. B.5 and B.6)

$$\text{gr}(\mathrm{H})+\text{gr}\left({\mathrm{C}\mathrm{H}}_3\mathrm{O}\right)\to {\mathrm{C}\mathrm{H}}_3\mathrm{O}\mathrm{H},$$(9)

and

$$\text{gr}(\mathrm{H})+\text{gr}\left({\mathrm{C}\mathrm{H}}_2\mathrm{O}\mathrm{H}\right)\to {\mathrm{C}\mathrm{H}}_3\mathrm{O}\mathrm{H}.$$(10)

This is dominant over 10⁶ years in component 4. For component 5, after 10⁴ years, gas-phase dissociative recombination becomes important:

$${\mathrm{C}\mathrm{H}}_3{\mathrm{O}\mathrm{C}\mathrm{H}}_4^{+}+{\mathrm{e}}^{-}\to {\mathrm{C}\mathrm{H}}_3+{\mathrm{C}\mathrm{H}}_3\mathrm{O}\mathrm{H}.$$(11)

We assumed these colder components are in the post-shock stage (see Fig. 8). Their model abundances are slightly overproduced until about 10⁵ years, where they sharply decrease. In the cold components, like component 2 above, we see a non-thermal desorption method is necessary. Unlike component 2, there are a number of different reactions that contribute to the formation of CH₃OH. There are sharp peaks in some of these routes, namely for Eq. (10). This reaction varies depending on the availability of grain surface H to react. For all components, the main destruction route is adsorption to the grain surface until about 10⁵ years. After that, in all components, there are contributions from gas-phase ion-molecule reactions.

To see if the timescales are physically realistic, we looked at the outflow. For a flow moving at 7 km s⁻¹ (from component N44 5, Table 2), to span 10′′ (Fig. 3; Zapata et al. 2012, Fig. 2), it would need to be 1.6 × 10⁴ years old. This is less than the time it takes for all five components to match model and observed abundances in the post-shock stage.

We find that in higher-density models the same evolution occurs on a faster timescale. We tried to use the RADEX densities in Model B as we had for Model A to reconcile the results of our different analyses. If the RADEX densities were used in the post-shock stage, the abundances would match at a few 10³ years for components 3, 4, and 5. However, there would be a sharp density jump from the shock stage to the post-shock stage in this scenario, for which we have no explanation. We tried using the RADEX densities in the shock stage, and the CH₃OH was destroyed quickly. However, in addition to the abundances not matching, we have no basis for such high densities in the shock stage as we used a previously parametrised shock. Different type shocks have been explored for example in L1157 (James et al. 2020), but it is out of the scope of this paper. We can only conclude that CH₃OH need some form of non-thermal desorption to produce the observed abundances.

We can speculate on possible sources powering non-thermal desorption in DR21(OH). There are numerous outflows seen emanating from MM1 and MM2, including the previously discussed molecular blueshifted outflow and the E-W CH₃OH maser outflows, which contain other species such as H₂CO and H¹³CO⁺(Zapata et al. 2012; Orozco-Aguilera et al. 2018) as well as thermal CH₃OH emission. A N–S H₂CS outflow is seen from MM1 with rotation temperatures of 40−60 K, where lines of CH₃OH are also detected (Minh et al. 2011, 2012). They suggest CH₃OH could arise from the same mechanism as H₂CS : non-thermal desorption from outflows interacting with the surrounding gas. While we cannot discern the exact sources of the non-thermal desorption in the CH₃OH components, the clustered cores of MM1 and MM2 are evidently producing multiple sites of turbulent gas. Further modelling of higher spatial resolution data is needed to determine the structure of DR21(OH) and connect it to the potential molecular formation routes.

5 Summary

We have presented IRAM and GBT observations mapping DR21(OH) at 10^′′ resolution. We used CH₃CCH, CH₃OH and H₂CO emission lines to determine the physical parameters of the different gaseous environments in the region, connecting the core scales to the larger clump and filament scales. Our main findings are:

• DR21(OH) main, or N44, is found to have ∼80 K and ∼30 K components near both MM1 and MM2, interpreted as the warmer inner core and the cooler envelope. The gaseous components of MM1 and MM2 are differentiated spatially on our maps and in velocity – MM1 is near −3.6 to −4.2 km s⁻¹, while MM2 is near −0.5 to −0.9 km s⁻¹. H₂CO and CH₃OH also show a broad, blueshifted outflow component that is not associated with the maser outflows of this region but with a molecular outflow seen previously in CO and SiO;

• Chemical modelling indicates that H₂CO can be produced in the warmup stage through thermal mechanisms, CH₃CCH can be produced in the warmup stage in an environment with low A_V, and CH₃OH needs a non-thermal desorption mechanism to produce the observed abundances. For CH₃CCH and CH₃OH especially, grain surface production is important for reproducing their gas-phase abundances;

• Other dusty condensations in the region show warm, ∼30 K, gas with similar line widths and velocities as those of N44. The cold components seen are likely non-LTE and represent sub-thermally excited gas.

With a sophisticated LTE model we are able to use large-scale single-dish data to disentangle different temperature and velocity components of star-forming clumps and cores. This is useful even for complicated species such as CH₃OH, as we are able to apply the model to the whole map along with targeted use of the more computationally expensive RADEX. The LTE model allows us to confidently determine the abundances of each of these species, which is essential for determining the chemical evolution of the region.

Data availability

The supplemental figures in the Appendix can be found at https://zenodo.org/records/15258578.

Acknowledgements

This work is based on observations carried out under project numbers 021-20 and 122-20 with the 30-m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). The Green Bank Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. We would like to thank the anonymous referee for their thoughtful and constructive comments. PF would like to thank L. Morgan and D. Frayer at the GBT/NRAO for their help developing observing scripts as well as GBTIDL scripts for ARGUS calibration and reduction, as well as P. Torne and M. Rodriguez at IRAM for their help observing with the IRAM telescope. PF and RP would like to thank K. Qiu and Y. Cao for providing the H₂ column density maps. PF and RP acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), through the Canada Graduate Scholarships – Doctoral, Michael Smith Foreign Study Supplement, and Discovery Grant programs. As researchers at the University of Calgary, PF and RP acknowledge and pay tribute to the traditional territories of the peoples of Treaty 7, which include the Blackfoot Confederacy (comprised of the Siksika, the Piikani, and the Kainai First Nations), the Tsuut’ ina First Nation, and the Stoney Nakoda (including Chiniki, Bearspaw, and Goodstoney First Nations). The City of Calgary is also home to the Métis Nation of Alberta Districts 5 and 6.

References

All Tables

All Figures

	Fig. 1 Integrated intensity maps for select transition lines of CH3CCH, H2CO, and CH3OH. The dense clumps as identified in Motte et al. (2007) are shown as stars, with the main clump N44 in black (see also Fig. 2). The quantum states and upper energy level of the transition are listed in the top left of each plot.
In the text

	Fig. 2 Regions isolated for LTE modelling near and in DR 21(OH), marked in white boxes over an integrated intensity map of the 12 K CH3CCH 50−40 line. From the top, the regions include: N40, N36, and N41 (labelled N); N44; N38; and N48. The dense clumps, as designated by Motte et al. (2007), are marked by stars.
In the text

Fig. 3 Column density maps, as resulting from an LTE model, for all species (rows) in components (columns) of N44. The cores SMA 7 in MM1 and SMA 3 in MM2 are labelled ‘1’ and ‘2’, respectively, and the clump of N44 is marked by a star (as designated in Zapata et al. 2012 and Motte et al. 2007). CH3OH 5 has been moved to align with H2CO 3, the broad components. Values for all LTE model parameters are listed in Table 2.

In the text

	Fig. 4 Column density maps for all species in components of the region N48. The clump location is indicated with a star (Motte et al. 2007).
In the text

	Fig. 5 Same as Fig. 4 but for the clump N38.
In the text

	Fig. 6 Same as Fig. 4 but for the N region, containing the clumps N40, N41, and N36.
In the text

Fig. 7 Model A, the two-stage warmup model, for H2CO (top) and CH3CCH (bottom). The 10 K cold cloud stage (left) is followed by a warmup stage where the LTE components are separated based on the temperatures and densities observed. The time axis represents the age of each component independently. The components of each model are in dotted, dashed, or solid lines; red or green lines are the model abundances in the gas phase, yellow the model abundances in the grain phase, and black the observed abundances calculated from the averages in Table 2.

In the text

Fig. 8 Model B, the three-stage shock model, for CH3OH, simulating a shock. The 10 K cold cloud stage (left) is followed by a sharp increase in density and temperature to simulate a shock (middle) and then a post-shock stage (right), when the gas settles to the LTE conditions observed. The time axis represents the age of each component independently. The LTE components are separated based on their temperatures and densities. The components of each model are in dotted, dashed, or solid lines; blue lines are the model abundances in the gas phase, yellow the model abundances in the grain phase, and black the observed abundances calculated from the averages in Table 2. As C1 and C2 are warmer, we suggest they are still in the shock stage, whereas the cooler C3, C4, and C5 are in a post-shock state.

In the text