© The Authors 2023

1 Introduction

Molecular clouds are the sites of star formation, where gas collapses into dense cores before forming protostars. The chemical makeup of molecular clouds therefore leads to that of protostellar systems. Observing the chemical complexity of star formation and modelling its evolution provides an invaluable link between the interstellar medium and planetary composition. In the last few decades, the diverse chemical makeup of star-forming regions has been revealed (Herbst & van Dishoeck 2009; Sakai & Yamamoto 2013; Jørgensen et al. 2020).

There are two major carbon-based groups commonly studied in these regions: interstellar complex organic molecules (iCOMs) and carbon chain molecules (CCMs). Complex, in interstellar terms, means any molecule with six or more atoms. COMs and CCMs are important to observe as their formation is sensitive to environmental conditions, thus can supplement existing information about interstellar objects traced by simpler, often studied, molecules such as CO, HCO+, or HCN.

COMs, which are carbon-bearing saturated complex molecules, have been observed in abundance in the hot (T > 100 K) and dense () cores of star-forming regions (see Herbst & van Dishoeck 2009 for a review of these molecules). The hot cores of high-mass star formation are well-known (high mass stars have M > 8 M_⊙; Blake et al. 1987; Caselli et al. 1993; Macdonald et al. 1996; Helmich & van Dishoeck 1997) as are their compact (≤ 100 AU in radius) equivalent in low-mass regions, dubbed hot corinos (Cazaux et al. 2003; Ceccarelli 2004; Bottinelli et al. 2004, 2007). In starforming regions complex molecules are formed on and released from interstellar grain surfaces. This chemistry is initiated by hydrogenation of CO, depleted onto grain surfaces from the gas phase, to make CH₃OH. COM formation in the gas phase is dominated by ion-neutral reactions, though the formation of saturated molecules in cold phases is inefficient (Herbst & van Dishoeck 2009).

CCMs, on the other hand, are unsaturated hydrocarbons (e.g. C_nH, HC_nN) that are well-known in cold (T ≤ 10 K) molecular clouds and starless cores such as Sagitarrius B2, Taurus Molecular Cloud, and Lupus-1A (Avery et al. 1976; Broten et al. 1978; Kroto et al. 1978; Little et al. 1978; Sakai et al. 2010). The cold starless core phase forms CCMs in the gas phase efficiently prior to the uptake of C in CO, a stable molecule. CCMs, then, were thought to have a short lifetime - present in the early stages of cloud core evolution but becoming deficient in later stage star-forming regions.

However, another formation route for CCMs is ‘Warm Carbon Chain Chemistry’ (WCCC; Sakai et al. 2008; Aikawa et al. 2008; Sakai & Yamamoto 2013). If atomic carbon depletes onto interstellar grains early in the star formation process, successive hydrogenation of C on the grain surface forms CH₄. As the gas warms in the star formation process and CH₄ is liberated from dust grains, it acts as the precursor for CCM formation. The condition keeping carbon in its atomic form prior to depleting onto grain surfaces is thought to be a short timescale for collapse (Sakai et al. 2008) or UV radiation (Spezzano et al. 2016; Higuchi et al. 2018). A WCCC source is characterised by two conditions: that there is an abundance of various carbon-chain molecules, and these molecules are concentrated in the warm and dense part around a protostar (Sakai & Yamamoto 2013). Following the detection of the low-mass star-forming region L1527 as a WCCC source (Sakai et al. 2008), WCCC characteristics have also been seen in a high-mass star-forming region, a giant H II region, and a starless core (Mookerjea et al. 2012; Saul et al. 2015; Wu et al. 2019).

Hot corinos and warm carbon chain chemistry (WCCC) sources now classify two extreme types of protostellar systems, due to this chemical differentiation. However, many sources have characteristics of both types – a hybrid source. The second WCCC source reported by Sakai et al. (2009), IRAS 15398–3359 in Lupus, hosts a hot corino (Okoda et al. 2023). iCOMs are emitted on a scale 10–100× smaller than the CCM emission. Okoda et al. (2023) suggest the hybrid nature supports both the grain mantle composition and the physical environment influencing the chemistry. L483, a Class 0 protostar, shows WCCC characteristics on a 1000 AU scale and also hot corino characteristics on (<) 100 AU scales (Sakai et al. 2009; Oya et al. 2017). The Bok globules B335 and CB 68, which are useful environments to study isolated protostellar sources, show hot corino emission on tens of AU scales and also carbon chain emission on hundreds to a thousand AU scale (Imai et al. 2016, 2022). Higuchi et al. (2018) find numerous intermediate sources in the Perseus Molecular Cloud, though find high CCH/CH₃OH ratios – characteristics of WCCC – towards isolated or edge-of-the-cloud sources owing to both differences in timescales and effects of UV radiation.

Carbon chain chemistry, traced by C₄H (Lindberg et al. 2016; Graninger et al. 2016), CCH (Higuchi et al. 2018; Bouvier et al. 2020), and cyanopolyynes (Taniguchi et al. 2018b, 2021), may be compared to complex organic chemistry, usually traced by CH₃OH, through parameters such as line density, column density, spatial distribution, or temperature. These studies identify that CCMs or their tracers are often in cooler, more extended gas, such as the envelope around hot cores, produced from CH₄ after it is sublimated from dust grain surfaces. There are dissimilarities between CCMs and COMs in spatial distribution, line profiles, and column densities. Lindberg et al. (2016), referencing Rodgers & Charnley (2001), are careful to note that certain classes of CCMs may have different correlations to methanol.

In this paper, we examine tracers of warm carbon chain and complex organic chemistry in the high-mass star-forming regions AFGL 2591 and IRAS 20126. Section 2 describes the molecules and sources. Section 3 describes the observations and Sect. 4 the detected lines and local thermodynamic equilibrium (LTE) model results. Section 5 discusses the results in wider context, and Sect. 6 provides a summary.

2 Objects of study

2.1 Warm carbon chain molecules -propyne and cyclopropenylidene

Propyne, or methyl acetylene, CH₃CCH, is part of the methylpolyyne family (Irvine et al. 1981). It 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 (Askne et al. 1984; Kuiper et al. 1984; Bergin et al. 1994). It was first detected in the giant molecular clouds Sgr B2 and Orion A, two sites with active star formation (Snyder & Buhl 1973; Lovas et al. 1976). It continues to be detected in low- and highmass star-forming regions and is used as a tracer of chemical complexity (Cazaux et al. 2003; Taniguchi et al. 2018b; Giannetti et al. 2017; Santos et al. 2022).

CH₃CCH can form both on grain surfaces and in the gas phase, with the former route directing observed abundances (Calcutt et al. 2019). On grain surfaces, CH₃CCH develops from the successive hydrogenation of C₃ (Hickson et al. 2016; Andron et al. 2018). In the gas phase, the primary route is through ion-neutral reactions (Schiff & Bohme 1979; Taniguchi et al. 2019) with contributions from neutral-neutral reactions (Turner et al. 1999). Recently, studies note the importance of CH₄, forming on then desorbing from grain surfaces in the warmup phase, and producing larger hydrocarbons through gas-phase ion-neutral reactions (Sakai & Yamamoto 2013; Calcutt et al. 2019; Taniguchi et al. 2019). These hydrocarbons dissociatively recombine to form CH₃CCH.

Thaddeus et al. (1985) identified cyclopropenylidene, c-C₃H₂, as the first astronomically observed organic ring molecule. The follow up in Vrtilek et al. (1987), reports on detections in Sgr B2, Orion, and TMC-1. c-C₃H₂ is formed in dense and ‘chemically young’ gas, where carbon is in the atomic form (Spezzano et al. 2016). Gas-phase ion molecule, H-atom transfer, and electron recombination reactions lead to its production (Park et al. 2006; Sakai & Yamamoto 2013). As this process starts with CH₄, desorbed from dust grains at 25 K, this molecule can be representative of WCCC. Aikawa et al. (2020) showed that with a new multi-layered ice mantle model abundances of c-C₃H₂ correlate with CH₄.

2.2 Complex organic molecules – formaldehyde and methanol

Formaldehyde, H₂CO, is not technically a complex species; however, it is chemically associated with larger organic molecules and we include it as a COM precursor in this paper. H₂CO was the first polyatomic organic molecule detected in space, Snyder et al. (1969) observed absorption spectra towards both galactic and extragalactic sources. Methanol, CH₃OH, is the simplest alcohol molecule. It has many low energy levels due to the torsional motion it undergoes as an asymmetric top molecule (Ball et al. 1970). Methanol is an abundant complex molecule with numerous transitions in the mm and sub-mm range which are useful, and commonly used, as a probe of the physical conditions of an interstellar region.

H₂CO and CH₃OH are both formed by the successive hydrogenation of CO after it freezes out onto the surface of grains in the cold cloud phase (Charnley et al. 1997; Garrod & Widicus Weaver 2013). The hydrogenation of CO and H₂CO require activation energy, suggesting not all molecules of these species will be converted into larger molecules. Still, this is a competitive process even at low temperatures; these species are discovered in a variety of sources such as cold clouds, hot cores, outflows, shocks, the Galactic centre, and external galaxies (see Herbst & van Dishoeck 2009 and references therein). These species are observed in both gas and solid phases and are key precursor molecules to larger COMs.

Certain CH₃OH transitions are masers, pumped into an excited state either collisionally (Class I) or by far-infrared radiation (Class II). CH₃OH masers are detected in IRAS 20126 (Class I and II) and AFGL 2591 (Class II) (Batrla & Menten 1988; Plambeck & Menten 1990; Harvey-Smith et al. 2008; Moscadelli et al. 2011; Rodríguez-Garza et al. 2017; Rygl et al. 2012).

Fig. 1 Herschel PACS 160 μπι image of AFGL 2591 (proposal ID KPGT_fmotte_1, PI F. Motte). The colour scale is logarithmic from 0.8 to 80 Jy pixel−1, while the contours are 1, 5, 10, 30, 40, and 60 Jy pixel−1. The black ‘x’ represents the dominant protostellar source VLA 3.

2.3 Sources

AFGL 2591 is the prototypical object in which to study physical and chemical processes during high-mass star formation. Figure 1 displays a Herschel PACS 160 μm continuum image of the source. At a distance of 3.33±0.11 kpc (Rygl et al. 2012), AFGL 2591 has a total mass of about 2×10⁴ M_⊙ and a measured total IR luminosity of 2×10⁵ L_⊙ (Sanna et al. 2012). Within this clump there are several radio continuum sources. VLA 1 and 2 are optically thin H II regions (Trinidad et al. 2003). VLA 3 is believed to be the youngest (dynamical age of about 2×10⁴ yr; Doty et al. 2002; Stäuber et al. 2005) and most massive source in the cluster, with an estimated mass of 38 M_⊙ (Sanna et al. 2012). VLA 4 and 5 are faint sources detected in cm continuum observations (Johnston et al. 2013). A large-scale east-west bipolar outflow powered by VLA 3 is seen in molecular observations which extends past the warm and dense envelope of the source (Lada et al. 1984; Hasegawa & Mitchell 1995; Sanna et al. 2012; Johnston et al. 2013). The simplest COMs and CCMs (CH₃OH, H₂CO, and CCH) are reported in AFGL 2591 via the James Clerk Maxwell Telescope (JCMT) Spectral Legacy and the CHESS surveys (van der Wiel et al. 2011; Kaźmierczak-Barthel et al. 2014). However, these surveys are at high frequencies (>300 GHz) where COMs and CCMs do not necessarily emit their strongest lines. Thus, little is known about the relative abundance of COMs and CCMs in this source.

IRAS 20126+4104 (hereafter IRAS 20126) is a well-studied massive protostar that allows us to observe the early stages of massive star formation in a relatively simple system. Figure 2 displays a Herschel PACS 160 μm continuum image of the source. In the Cygnus X region at a distance of 1.6 kpc, IRAS 20126 has a luminosity of 1.3 × 10⁴ L_⊙ (Moscadelli et al. 2011). The natal dense clump, which is a site of massive star formation, is isolated with a temperature of 40 K (Shepherd et al. 2000). Within this clump observations of C³⁴S by Cesaroni et al. (2005) confirm the presence of a large Keplerian disk with a radius of 5000 AU. Near-infrared K, L′, and M′ band observations by Sridharan et al. (2005) with the United Kingdom Infrared Telescope, however, resolve the central source and suggest the presence of a more compact disk (R ~ 1000 AU). This was confirmed by Chen et al. (2016) whose modelling of SMA observations of CH₃CN confirm a 1.5 M_⊙ disk with a radius of 850 AU rotating about a 12 M_⊙ protostar. Numerous COMs have been detected in JCMT and Plateau de Bure Interferometer (PdBI) observations (Isokoski et al. 2013; Palau et al. 2017) which have also included the detection of some simple CCMs (CCH, CH₃CCH). However no surveys targeting CCMs have been reported. Thus, little is known about the relative abundance of COMs versus CCMs in this source as well.

Fig. 2 Berschel PACS 160 μm image of IRAS 20126+4104 (proposal ID OT1_rcesaron_1, PI R. Cesaroni). The colour scale is logarithmic from 0.4 to 40 Jy pixel−1, while the contours are 1, 5, 10, 20, and 30 Jy pixel−1. The black ‘x’ represents the protostellar source.

3 Observations

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). Using the on-the-fly observational mode with position switching, we mapped a 1′ × 1.5′ region in both AFGL 2591 and IRAS 20126 in the frequency range 287.22–295 GHz. The offset position was −120″ horizontally and – 240″ vertically. The phase centre was αJ(2000) = 20^h29^m24^s.9 and δJ(2000) = +40°11′20.″0 for AFGL 2591 and αJ(2000) = 20^h14^m26.04 and δJ(2000) = +41°13′32.″5 for IRAS 20126. The EMIR receiver, in bands E1 and E3, was connected to the FTS200 backend (Fast Fourier Transform Spectrometer at 200 kHz resolution). The range 131.2–138.98 GHz was observed simultaneously in this set-up, but will not be discussed in this paper due to the much larger beam size. The pointing, checked every 90 min or less, was done using either K3-50A or NGC 7027 and resulted in typical corrections of 1–2″. The focus, completed every 3 h or after a sunset or sunrise, was done with the same source as for pointing and resulted in typical corrections of <0.3″. The atmospheric opacity, as measured by the 225 GHz taumeter, aside from part of the day November 28, 2020, was <0.3, and often <0.1. The system temperatures observing IRAS 20126 were 350–500 K, and observing AFGL 2591 were 300–350 K. The beam size is 9.3″ at 290 GHz, with a pixel size of 4.4″. Data reduction, and the production of maps, was completed with GILDAS/CLASS¹.

The 100 m Green Bank Telescope 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 in AFGL 2591 and IRAS 20126, centred as above, utilising position-switching to a location +3′ in azimuth for the reference measurements. Each petal, with a map radius of 0.5′, took 54.7 s with a scan rate of 2.3″ per second and an integration time of 1 second. The pointing and focus was completed every 30–40 min with the X-band receiver using the GBT’s automatic scan AutoPeakFocus on source 2015+3710. At the ARGUS frequencies, the automatic scan AutoOOF was used to correct for errors in the reflecting surface using source 1229+0203. The system temperatures ranged from 100 to 200 K. With spectrometer mode 2 we obtained a spectral sampling of 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″ or 10.0″ at 96 or 85 GHz respectively, with a pixel size of 2″. The data were reduced and calibrated using GBTIDL².

The spectral axis of the IRAM data was smoothed from the native spectrometer resolution to reduce the high levels of noise. The final spectral sampling is 781 kHz, or 0.80 km s⁻¹. Given the widths of the spectral lines, this is still adequate to confidently detect and fit these lines. The GBT 85 GHz data was spatially smoothed by a factor of two in pixel size to improve the noise levels and to closer match the IRAM pixel size. With this, and the higher quality GBT observations, the GBT data was not spectrally smoothed as the noise levels are lower. The GBT 96 GHz and IRAM 290 GHz data were regridded, smoothed, and rebinned in CASA³ to match the 85 GHz beam size and spatial grid. The final restored beam is 10.0″ with a pixel size of 4″. Average rms noise levels are shown in Table 1. In the Results, Sect. 4, the data are in units of T_A*, adjusted during calibration. In the LTE model, the data are converted to T_mb 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.

4 Results

4.1 Detected lines

Species were identified through the CASSIS⁴ LTE+RADEX⁵ module. First, the CDMS and JPL catalogues under the ‘Species’ tool in CASSIS allowed us to determine possible detected species across the observed spectrum. These two catalogues are cross-referenced for confirmation. Next, we used the LTE Model in the module to model the spectrum. We input an estimated column density, excitation temperature, FWHM line width, source size and υ_LSR. The spectrum was computed across the frequency range using telescope parameters specific to IRAM and the GBT as documented within CASSIS. Line identification was then done visually, matching the model to the data. It is important to note that this model is only for line identification and not to accurately determine the physical conditions.

The combined GBT and IRAM data set contain several simple and complex molecules. The complete spectral surveys are to be described in Freeman et al. (in prep.). Table 2 lists all possible CH₃CCH (CDMS tag 40502 for the vibrational ground state), c-C₃H₂ (CDMS tag 38508), CH₃OH (CDMS tag 32504), and H₂CO (CDMS tag 30501) lines with A_ij > 1×10⁻⁶ s⁻¹ and E_up < 150 K. For CH₃OH the A_ij limit was set to 2×10⁻⁶ s⁻¹ to remove the 40 K maser line. For CH₃CCH the E_up limit was increased to 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. This range was selected to restrict the catalogue transitions to the ones that could be detected with the sensitivity of our observations. The species types - a- and e-CH₃CCH, o- and p-c-C₃H₂, A- and E-CH₃OH, o- and p-H₂CO – were not differentiated as we do not detect enough lines of each. The VASTEL database information for each type is similar to that of CDMS, aside from including higher E_up a-CH₃CCH lines, and we expect the results to not be affected drastically by using the combination of the types.

Figures 3 and 4 show examples of the spectra at the brightest pixel for the CH₃CCH 5₀−4₀ (E_up = 12 K) line and the CH₃OH 2₀−1₀, 2₋₁−1₋₁, 2₀−1₀ triplet (E_up = 7, 13, 20 K). Figures A.1, A.2, A.4, and A.5 display all transition lines of each molecule. The propyne detections with E_up from 10–30 K are easily detected and well-defined; rarely showing blended or complicated structure. The methanol lines from six to several tens of K are clearly detected. The higher energy lines are often blended – in part due to methanol’s complicated line structure – and often show structure that is not represented by a single Gaussian. The IRAM data also suffer from higher noise due to unstable weather conditions and higher observed frequencies. These trends are consistent across both IRAS 20126 and AFGL 2591, with the latter displaying more complicated methanol line profiles.

Figures 5 and 6 show integrated intensity maps for select lines covering a range of upper energy levels. These are notably bright in the lower energy lines, and decrease in strength as the upper state energy of the lines increase. Upwards of a few tens of K, there is little structure to be seen in the carbon chain molecules. In AFGL 2591 (Fig. 5), there are differences in spatial distribution among the different molecules as well as between the different energy level lines. Especially in CH₃OH, the warm lines of E_up = 20–50 K trace the cloud near the known sources VLA 1–5, while the colder 6–20 K lines trace a region, known to be abundant in methanol (referred to as the ‘methanol plume’ by van der Wiel et al. 2011), to the north and north-east. The similarity between molecules is in the quite extended emission around the hot cores region.

In IRAS 20126, all maps are quite featureless. The COM distributions, across energy levels, have a slight south-east-north-west tilt, and are all similarly concentrated around a single region, coinciding with the hot core. The CCM lines show a similar spatial concentration however with a slight south-north extension in the lines detected, as opposed to the south-east-north-west extension seen in CH₃OH. The integrated intensities of the carbon chain lines drop off in strength considerably from the lowest energy lines of 12 K and 6 K, for CH₃CCH and c-C₃H₂, respectively.

The beam size, 10.0″, as displayed in the integrated intensity maps, corresponds to 3.3 × 10⁴ AU (0.16 pc) in AFGL 2591 and 1.6 × 10⁴ AU (0.08 pc) in IRAS 20126 which is sufficient to resolve the extended molecular emission of tens of thousands of AU in each source (Cesaroni et al. 1999; van der Wiel et al. 2011).

Fig. 3 Select transition lines of the carbon chain molecules, CH3CCH and c-C3H2, the complex organic molecules CH3OH, and the precursor complex organic molecule H2CO, at the central pixel of the AFGL 2591 GBT and IRAM maps. The chosen lines represent a range of different upper energy levels, which are indicated in black on the plot. The dashed grey line represents the observed transition frequency at this location. The observed data is shown in black. The model spectra, to be discussed in Sect. 4.2 are shown in colour. The red solid line represents the total model, and the sole component if only one component is fit. For the multi-component fit in CH3OH the separate components are shown in blue, yellow, and brown dashed lines.

Fig. 4 Same as Fig. 3 for the central pixel of the IRAS 20126 GBT and IRAM maps. For the multi-component fits in c-C3H2 and CH3OH the separate components are shown in blue and yellow dotted lines.

4.2 LTE modelling

4.2.1 Model description and implementation

We used a more comprehensive LTE model to determine the physical conditions of AFGL 2591 and IRAS 20126. The model, developed in Python, was based on the CASSIS LTE model scripts in Jython. This formalism is fully documented in Vastel’s Formalism for the CASSIS software⁶. The LTE model simultaneously fit numerous transitions of one or more molecules across a spectrum for one or more physical components – defined by the physical parameters source size, excitation temperature, column density, spectral line width, and source Doppler velocity – of a source.

The model identified the molecular transitions producing spectral lines in a given frequency range. Then, it generated spectral line profiles for any desired molecule iterating over the physical parameters of the gas until it produced a model spectrum that best matches the data. The resulting best fit, therefore, provided these five physical parameters for the gas emitting the observed spectra. As we produced maps with our observations, the code was modified to loop over multiple pixels modelling all the observed spectra, ultimately producing spatial maps of the physical conditions.

The best fit was determined via the Levenberg–Marquardt algorithm in LMFIT (Newville et al. 2014), a curve-fitting method for non-linear least-squares problems. The LevenbergMarquardt algorithm is sufficiently quick for finding the minimum value, and we trust that with appropriately restricted parameter ranges (initial estimations are described below and adjusted in successive iterations of the model) we found the global minimum. The goodness-of-fit was reported with the reduced χ², which describes the least-squares minimisation and takes into account the number of data points and variables included in the model. For an LTE model, there may be degeneracy between parameters such as column density and temperature. In some cases, we have few spectral lines compared to the number of free parameters and cannot exclude the possibility of degeneracy. However, we present a best effort to produce reasonable results given this limitation. In other cases, we have several strong lines or have fixed certain parameters in order to reduce the discrepancy between data points and free parameters. This is discussed in the LTE model results below.

In the analysis, we excluded the possible CH₃OH maser line at 84.5 GHz (5_1,5,2 − 4_0,4,1), seen in Kalenskii et al. (2002), for example. While determining maser sources is outside the scope of this paper, our initial LTE models for identifying lines greatly under-represented the 84.5 GHz line strengths in several pixels with strong CH₃OH emission. Thus, we removed it in case it was a maser. Within our frequency ranges, no other methanol masers are reported in these sources.

We restricted the number of pixels the model loops through with a signal-to-noise ratio (S/N) mask for each molecule of interest. We fit a Gaussian to the spectrum using Astropy (Astropy Collaboration 2022), calculated the peak of the observed spectrum using the Gaussian1D fitter and calculated the noise from a nearby, line-free frequency range with the SigmaClip function. We fit specific, strong lines of each molecule to create this mask. Any pixel that had a S/N > 2.5 for the CH₃CCH 12 K 5₀−4₀ (E_up = 12 K) and 5₁−4₁ (E_up = 19 K) lines at 85 GHz are included in the CH₃CCH and c-C₃H₂ LTE models, and any pixel that had the same threshold for the CH₃OH 2_0,2,0−1_0,1,0 (E_up =7 K) and 2_1,2,2−1_1,1,2 (E_up =12 K) CH₃OH lines at 96 GHz or the 6_0,6,0−5_0,5,0 (E_up =48 K) and 6_1,6,2−5_1,5,2 (E_up = 54 K) lines at 290 GHz (to capture lines in both GBT and IRAM data) are used for the CH₃OH and H₂CO models. The S/N threshold of 2.5 is lower than the standard 3 and was used to ensure we capture the cloud edges where, from a visual inspection of the spectral lines, we believed the emission is still real. As well, the corresponding c-C₃H₂ and H₂CO masks produced more extended maps. To allow comparison between molecules we continued with the more restrictive CH₃CCH and CH₃OH maps with a loosened S/N threshold.

To begin the analysis, we extracted the spectrum from the brightest pixel of certain molecular transition lines - CH₃CCH 5₀−4₀ (E_up = 12 K); c-C₃H₂ 2_1,2−1_0,1 (E_up =6 K); CH₃OH 2_0,2,0−1_0,1,0, 2_1,2,2−1_1,1,2, 2_0,2,1−1_0,1,1 triplet (E_up = 7, 12, 20 K), CH₃OH 6_0,6,0−5_0,5,0 (E_up =48 K); H₂CO 4_0,4−3_0,3 (E_up = 35 K) – and fit one or more Gaussians to these lines with the CASSIS Line Analysis module. This module, with a chosen molecule, frequency range, and transition thresholds (such as E_up,max and A_ij,min), selected and displayed all possible transition lines in the data. By fitting Gaussians by hand, we determined the number of possible components to use in the LTE model and used the line width and velocity as guides for the input parameters. The temperature and column density, on the other hand, are physical properties of the source and input ranges are estimated based on previous knowledge of these sources as well as our integrated intensity maps, Figs. 5, and 6. The size of the source was set at 20″, assumed to be two times the size of the beam. This was to account for the extended emission seen by single dish telescopes, and it was appropriate as we are not resolving the compact protostellar sources but the larger scale envelope. Assuming a source size larger than the beam also removed any beam dilution effects when modelling the molecular spectral lines. With smaller (10″) beam sizes, the resulting column density changed but remained within a factor of 3. With larger (30″) beam sizes, the column density changed by less than a factor of 2. Otherwise, the results did not change and we maintained this constant beam size.

Fig. 5 AFGL 2591 integrated intensity, , maps of rotational transition lines spanning a range of energy levels. The upper state energy level is indicated in the top right-hand corner and the beam size in the bottom left-hand corner. The source VLA 3 is noted by a white or black ‘x’. CH3CCH, in the first column, c-C3H2, in the second, H2CO, in the third, have contour levels of 10, 20, 30, and 40 times the noise level of the respective frequency range. CH3OH, in the fourth column, has contour levels of 10, 20, 40, 80, and 120 times the noise level of the respective frequency range.

Fig. 6 Same as Fig. 5 for IRAS 20126, the source noted by a white or black ‘x’. CH3CCH, in the first column, has contour levels of 10, 30, 50, and 70 times the noise of the respective frequency range. c-C3H2, in the second column, H2CO, in the third, and CH3OH, in the fourth, have contour levels of 10, 40, 70, and 100 times the noise of the respective frequency range.

4.2.2 AFGL 2591 LTE model results

Initially, each molecule was modelled individually with one component, unless the CASSIS Gaussian fitting described suggested there were multiple components - either through different velocity components, line wings, or otherwise non-Gaussian shape. If multiple components were needed, certain values were fixed to balance the number of transition lines available and the number of free parameters in the model. The input parameters for all models are shown in Table A.1.

The initial one-component fit of c-C₃H₂ alone provided unrealistic excitation temperatures (minimising at the lower limit of 3 K) and corresponding large column densities that did not align with the observed emission. There is only weak detection of the lowest energy (E_up = 6 K) c-C₃H₂ transition line which dominates the fits. Given this, we were unable to model with c-C₃H₂ alone and it was modelled in the same component as CH₃CCH in the CH₃CCH one-component model. We assumed the two carbon chain species are originating from the same gas in the extended warm and cold extended envelope. The additional transition lines provided further constraints on the resulting excitation temperature and column density from the model. When combined with CH₃CCH, c-C₃H₂ reproduces the same structure of the integrated intensity maps of c-C₃H₂ in the column density maps, indicating that the column density reflects where we see emission in the data. Previous models, with c-C₃H₂ alone, produced unrealistic column density distributions. In addition, the model fits the spectral lines for each molecule well (example in Fig. A.4), supporting the assumption that these two molecules are coming from the same velocity and temperature gas. In further support, these are different conditions from the other two molecules, described next. As will be described in the following section for the results of IRAS 20126, c-C₃H₂ could also show some non-LTE behaviour.

The CH₃OH spectral line profiles contain more velocity structure than compared to the other molecules. CH₃OH was fit unsuccessfully by a one-component fit, as determined by non-converging models, high reduced χ² values, or unphysical values for the excitation temperature or column density. The low energy GBT lines were often over-represented compared to the higher energy IRAM lines, and produced very low excitation temperatures unexpected for LTE.

Given the possible non-LTE nature as indicated by the low excitation temperatures found, we tried to model CH₃OH in RADEX (van der Tak et al. 2007) with CASSIS. Unlike the LTE model, RADEX is too intensive to run over all pixels in the map. Thus, select points near the source VLA 3 and along the eastward extension of the molecules are modelled, along with one pixel off-source to the north, as seen in Fig. 7. Table 3 shows the results, where the model delineates the source into two components - one which traces the higher column density (10¹⁴ cm⁻²), warmer (20 K), gas near the source υ_lsr of −5.7 km s⁻¹ represented by the higher energy IRAM lines, and a second which traces the colder (10 K), lower column density (10¹³ cm⁻²) more negative velocity (−7 to −9 km s⁻¹) gas. Figure 8 shows the resulting model spectra for pixel (38,31).

The kinetic temperature from the RADEX model is then used to guide the input excitation temperature, a key parameter for this molecule, for the LTE model. Since the kinetic and excitation temperatures matched, the assumption of LTE was valid and we used the LTE model to produce parameter maps for CH₃OH. In subsequent LTE modelling, the latter cold component was able to be differentiated into two components in the LTE model - the CH₃OH ‘blob’ towards the north-east below in component two and the blue-shifted component three, both described below.

H₂CO, lastly, is well-represented by a one-component model where the parameters, aside from size, are free. Taking direction from the CCMs, we tried to fit H₂CO and CH₃OH in the same component of the model to represent these molecules originating from the same gas. However, none of the multiple colder components of CH₃OH aligned with the hotter H₂CO emission.

The parameter maps, for velocity, excitation temperature, column density, and line width (FWHM) of the four molecules are displayed in Fig. 9. For all parameter maps, the first component is shown for each molecule. This component may also be seen as the ‘main’ component, as it is most closely associated with the protostar and the envelope. The second and third components for CH₃OH are displayed in Fig. A.3.

All four molecules show a fairly constant velocity structure near the source υ_lsr of −5.7 (±0.03−0.5) km s⁻¹. Component three of CH₃OH, however, is notably blue-shifted and fixed at −9.4 km s⁻¹. Component two of CH₃OH, representing the offset methanol ‘blob,’ has no contribution towards the centre leaving a hole in the map. Towards the edges, where the detected lines become weaker, this value is slightly blue-shifted.

The carbon chain molecules show a fairly uniform temperature of 30 (±2–4) K across the main emitting areas of the source. The excitation temperature drops to near 10 K where the all lines are weaker and the c-C₃H₂ lines dominate. The temperature of H₂CO reaches towards 70 (±5) K slightly offset from the VLA 1–3 sources and shows a warm temperature across extended spatial scales. For CH₃OH, not displayed, the excitation temperatures were fixed at 18 K for component one and 8 K for components two and three.

For all components, the column density maps trace features of the integrated intensity maps in Fig. 5. CH₃CCH, at a peak of 4 (±0.4) × 10¹⁴ cm⁻², extends along and towards the east of AFGL 2591, with two distinct peaks. c-C₃H₂ has the lowest column density, 2 (±0.3) × 10¹³ cm⁻², and shows an extended distribution with a slight peak near the VLA sources. H₂CO peaks at 9 (±0.03) × 10¹³ cm⁻², and distinctly is distributed near the VLA sources. CH₃OH peaks at 5 (±0.2), 5 (±0.2), and 2 (±0.1) ×10¹⁴ cm⁻² for components one, two, and three, respectively. Component one has the highest column density of all molecules and peaks slightly offset of the protostars, similar to the other molecules. Component two shows the extended emission to the north-east found in the lower energy transition lines, while component three shows a less dense, blue-shifted component just north of the source.

The line widths of the carbon chain molecules are around 3 km s⁻¹ wide across the emitting region, and H₂CO is uniform at a slightly broader width near 4 km s⁻¹. Components two and three of CH₃OH were fixed at 3.9 km s⁻¹ and 2.5 km s⁻¹ respectively, while component one results in line widths of 2–3 km s⁻¹.

While all our detected lines represent extended envelope emission, the differing temperatures in LTE suggest the molecules are in different gas. As well, the multiple components reveal different structures in the region, which is examined in Sect. 5.1.

Fig. 7 Pixels fit with a RADEX model in AFGL 2591 marked on a Herschel PACS continuum image. The black and white numbers correspond to said pixels, with results listed in Table 3. The black ‘x’ represents the source VLA 3.

Fig. 8 RADEX fit for CH3OH in the central pixel of AFGL 2591, corresponding to pixel 7 in Fig. 7 and Table 3. The solid red line is the overall two-component model, the solid black line is the data, the dashed blue line is component one and the dashed yellow line is component two. Component two is only seen in the lower energy lines within the GBT range, thus component one and the overall model are the same for the lines in the IRAM range. The dashed grey line represents the observed transition frequency, with the energy level of this transition labelled.

Fig. 9 Parameters maps of AFGL 2951 showing the results of the LTE models. The black ‘x’ represents the source VLA 3. For CH3OH, Tex was fixed at 18 K.

Fig. 10 Parameters maps of IRAS 20126 showing the results of the LTE models. The black ‘x’ represents the protostellar source.

4.2.3 IRAS 20126 LTE model results

As in AFGL 2591, each molecule is initially modelled individually. For IRAS 20126, this produced the best results and the presented analysis keeps all molecular models separate. Contrary to AFGL 2591, the c-C₃H₂ and CH₃CCH combination did not represent the lines well. The c-C₃H₂ lines in IRAS 20126 are stronger, and drove the excitation temperatures too low for CH₃CCH. CH₃OH and H₂CO were also kept separate, as the slight disparities in temperature worsened the fits if combined. Given the limitations of the number of components we could discern with this data and this modelling, there may be shared components which we do not see.

Each molecule is modelled with one component unless the CASSIS Gaussian fit suggested there were multiple components. Figure 4 displays the model results for select transition lines covering a range of upper energy levels. For CH₃CCH and H₂CO, one component was adequate. For c-C₃H₂ and CH₃OH, a second narrow, cold component was added to better represent the detected lines. This component, in both molecules, is a foreground screen across the source in the lower energy transition lines - as seen in the c-C₃H₂ 2_1,2−1_0,1 (E_up = 6 K) and CH₃OH 2_0,2,0−1_0,1,0 (E_up = 7 K) lines in Fig. 4. In these cases of multiple components, certain values were fixed to balance the number of transition lines available and the number of free parameters in the model. In giving the errors of the varying components, we present the statistical errors from the model in the range seen across the map, excluding edges.

The parameter maps, for velocity, excitation temperature, column density, and line width of the four molecules are displayed in Fig. 10. For all parameter maps, the main or first component is shown for each molecule. The second components of c-C₃H₂ and CH₃OH are displayed in Fig. A.6.

There is a consistent velocity range for all species between −2.0 to −4.5 (±0.02−0.1) km s⁻¹. This is a range around the source V_lsr of −3.5 km s⁻¹. For the CCMs this is oriented with a red-shifted northern and blue-shifted southern feature. For the COMs, this is oriented with a red-shifted north-west and blue-shifted south-east feature. In c-C₃H₂, the second component was fixed at the same velocity as the first component, and in CH₃OH the second component was fixed at −0.2 km s⁻¹ from the first.

Similar to the velocity, the excitation temperature structure is differentiated between categories of molecules. The carbon chains have a roughly consistent temperature across the source, with CH₃CCH at 25–30 (±0.9–4) K and c-C₃H₂ at 8–12 (±0.8–3) K. The COMs show a warmer central region falling off in temperature towards the edges, with a range of 15–35 (±0.3–2) K for CH₃OH located just off the protostar, and 55–85 (±3–6) K for H₂CO located on the protostar. The second component of c-C₃H₂ varies from 4–8 (±0.4−6) K with few pixels at 10–12 K, while the second component of CH₃OH was fixed at 10 K.

Given the low excitation temperature of c-C₃H₂, it is possible that it is not in LTE and these results may not be a true representation of the physical gas temperature. To explore this, we used RADEX to model a pixel that was bright in c-C₃H₂, as shown in Fig. 11. We reproduced the spectrum using a temperature of 30 K, as expected from the CH₃CCH results, an ortho- to pararatio of 3, and a density of 10⁵ cm⁻³.

All column density maps peak near the site of the protostar, with the COMs showing a steeper falloff towards the edges. CH₃CCH, CH₃OH, and H₂CO peak near 6, 9, and 3 (±0.3, 0.2, and 0.1) × 10¹⁴ cm⁻² respectively, while c-C₃H₂ peaks near 5 (±0.3) × 10¹³ cm⁻². The second component of c-C₃H₂ was fixed at 1 × 10¹³ cm⁻² and for CH₃OH it was fixed at 2.1 × 10¹³ cm⁻².

The CCM lines are consistently narrower in line width, ranging from about 1–2.5 km s⁻¹, while the COM lines are 47 km s⁻¹ wide. The former show a slightly wider line profile near the protostar, while the latter show a widening towards the southeast of the source. The second component of c-C₃H₂ was fixed at a line width of 0.5 km s⁻¹ while for CH₃OH it was fixed at 1.2 km s⁻¹.

The visual and quantitative differences in the results, namely velocity and temperature, between the types of molecules suggests they are coming from different physical environments. We elaborate on the known and observed structure in Sect. 5.2.

Fig. 11 RADEX fit for c-C3H2 in the central pixel of IRAS 20126. The model used a temperature of 30 K, as expected from the CH3CCH results, an ortho- to para- ratio of 3, and a density of 105 cm−3. The model is in red while the data is in black. The upper state energy level of each transition is shown on each plot.

5 Discussion

Detections of carbon chains in high-mass star-forming regions have increased over recent years, with candidate, but not confirmed, WCCC sources (Saul et al. 2015; Taniguchi et al. 2018a). Mookerjea et al. (2012) present the first trace of WCCC in a high-mass star-forming region with detections of the simple chain C₃ in DR21(OH), another source with hot cores in Cygnus X. The C₃ abundance is consistent with the warm-up model of WCCC in Sakai et al. (2008). The emission and the resulting temperatures of 30–50 K – a ‘lukewarm corino’ – are consistent with the warm envelope around the hot core. The extended spatial distribution and the low, but warm, temperatures are also seen in our sources AFGL 2591 and IRAS 20126.

Comparing carbon chain molecules to CH₃OH is a recently available and practised method of analysing star-forming regions through molecular observations. Table 4 summarises several studies that used either CH₃CCH or c-C₃H₂, as in this study. At first glance, there is a spread of N(CH₃CCH)/N(CH₃OH), possibly due to the variety of studies and sources that examine this ratio. All CH₃CCH studies presented used single dish observations and resolved objects on scales of 1–7×10⁴ AU, comparable to this study, and consistent with this species being found in the cold and warm extended envelopes. In Fayolle et al. (2015), CH₃CCH is only detected with the IRAM 30m and not with higher resolution Submillimeter Array observations. There is also variation within the protostellar age – Fayolle et al. (2015) note their sources, showing N(CH₃CCH)/N(CH₃OH) > 1, are potentially younger than typical hot core sources, while Taniguchi et al. (2018b) find N(CH₃CCH)/N(CH₃OH) > 1 for their source with a UCH II region. For c-C₃H₂, Higuchi et al. (2018) look at scales of a few thousand AU, a smaller region, yet they find variation within the survey of protostars and comparable values to our sources.

In this work, at the continuum peak of IRAS 20126 N(CH₃CCH)/N(CH₃OH) = 0.99 while N(c-C₃H₂)/N(CH₃OH) = 0.06. For AFGL 2591, N(CH₃CCH)/N(CH₃OH) = 0.97 while N(c-C₃H₂)/N(CH₃OH) = 0.05. The column density ratios are shown across the source maps in Fig. 12. The ratio maps in IRAS 20126 are fairly uniform and structured compared to AFGL 2591, which could arise from the smaller spatial extension seen in IRAS 20126, at half the distance of AFGL 2591. However, there is a wide range of ratios within a single source and not just from source to source, making cloud to cloud comparisons somewhat unreliable.

Environmental differences cause the variation between CCMs and CH₃OH. In our data, the column densities of CH₃CCH and CH₃OH are of comparable magnitude, and CH₃CCH is typically seen on a uniform, extended scale in comparison to concentrated CH₃OH closer to protostellar sources. Santos et al. (2022) most recently compare N(CH₃CCH)/N(CH₃OH) in the massive hot molecular core G331.512-0.103 with single-point APEX observations and find a ratio of 0.42. The authors compare abundance ratios to that of other high-mass star-forming regions, with values across a range of 0.31–2.2 (with error, Taniguchi et al. 2018b; Giannetti et al. 2017; Fayolle et al. 2015). CH₃OH is likely to be more spatially compact than CH₃CCH, however. The latter is typically an envelope species, as in Fayolle et al. (2015) where the ratios are > 1. This study focuses on organic poor massive young stellar objects, rather than CH₃OH rich hot cores, though the presence and relative abundance of CH₃CCH and other complex species makes them comparable to, or possible precursors to, high mass hot cores.

c-C₃H₂ is often an order of magnitude lower in column density than CH₃OH, with no correlation between the two molecules. This is not surprising, due to the different formation routes but has been discussed in previous observations. Spezzano et al. (2016) find different environmental conditions for c-C₃H₂ and CH₃OH in the prestellar core L1544, with interstellar radiation keeping C in its atomic form and allowing carbon chains to form near the edge of the source - a new view of carbon chain chemistry. In the shielded region of the cloud, CO forms, leading to the production of complex organic molecules. Thus, there is a stark difference in spatial distribution between the two molecules (their Fig. A.1). L1554, however, is a much different environment than our two protostellar regions. In general, as a collapsing core evolves, the gas and dust environment will change. The original view of WCCC, which will be discussed for our sources in Sect. 5.3, is dependent upon the timescale, where carbon chains are reproduced in the warm protostellar environment. More specifically, AFGL 2591 and IRAS 20126 are relatively isolated sources and have no external irradiation sources to affect chemical processes. For a more direct comparison, Higuchi et al. (2018) find no correlation between c-C₃H₂ and CH₃OH for 36 protostars in the Perseus region. The column density ratio they find is on order of magnitude 10⁻²−10⁻¹ similar to what we find in both IRAS 20126 and AFGL 2591.

Fig. 12 Column density ratio maps for CH3CCH, c-C3H2, and H2CO with respect to CH3OH. The upper panel shows AFGL 2591 and the lower panel IRAS 20126. The white ‘x’ represents the protostellar sources within each region.

5.1 AFGL 2591

The environment of AFGL 2591 has well-known features: sources VLA 1–5, a warm inner envelope, and an extended methanol plume (Hasegawa & Mitchell 1995; van der Wiel et al. 2011; Jiménez-Serra et al. 2012; Gieser et al. 2019). An extended distribution around this region is also seen in smaller molecules such as CO, HCN, HCO+, H₂CO (van der Tak et al. 1999). Traces of these features are able to be seen from our LTE model (Fig. 9). The sources VLA 1–5 are not spatially resolved in our single dish observations, and the structures present in both our integrated intensity maps (Fig. 5) and column density maps (Fig. 9) - aside from the CH₃OH 2_0,2,0−1_0,1,0 (E_up = 7 K) line and component two CH₃OH maps, which trace the methanol ‘plume’ – are showing the smoothed out sources as well as the inner envelope. The dense inner envelope is shown to be a dominant component for transition lines with E_up ≤ 200 Κ (van der Wiel et al. 2013).

VLA 3, an early B-type star and the central heating source (Trinidad et al 2003), is also connected to a disk (van der Tak et al. 2006; Wang et al. 2012; Gieser et al. 2019) as well as large scale (>1′) outflow blue-shifted towards the south-west and red-shifted towards the north-east, traced in simple molecules (Hasegawa & Mitchell 1995; Jiménez-Serra et al. 2012; van der Wiel et al. 2013; Gieser et al. 2019). van der Wiel et al. (2011) detect the methanol plume extending past this inner envelope towards the north-east. The CH₃OH 2_0,2,0−1_0,1,0 (E_up = 7 K) (Fig. 5) line, which is blended with the 2_1,2,2−1_1,1,2 and 2_0,2,1−1_0,1,1 (E_up = 13, 20 K) lines, clearly traces this plume in multiple energy levels and at greater intensity than around VLA 1–5. The 6_0,6,0−5_0,5,0 (E_up = 48 K) line shows greatest intensity around VLA 1–5 but still extends into the NE, further supporting the findings from van der Wiel et al. (2011).

The colder temperatures, <20 K, for gas-phase CH₃OH indicate there are likely non-thermal desorption methods liberating this molecule after formation on dust grain surfaces. Cold methanol, at different scales and differing spatial features, is seen in numerous sources with differing mechanisms proposed. Bouvier et al. (2020), with single dish observations of Orion Molecular Cloud 2/3, conclude the detected CH₃OH is originating from the photo-dissociation regions rather than the hot corino.

Vastel et al. (2014), with single dish observations of the prestellar core L1544, conclude that non-thermal desorption methods, specifically photo-desorption, is responsible for the cold CH₃OH detected. Soma et al. (2015), similarly in Taurus Molecular Cloud 1, offer several possibilities for non-thermal CH₃OH desorption: weak shocks, photoevaporation, and chemical desorption. Favre et al. (2020) find a methanol ‘blob’ at smaller scales (1000 AU) using interferometric observations of L1521F – a very low luminosity object transitioning from prestellar to protostellar stages – with the likeliest mechanism for desorption being a slow, gentle shock. In AFGL 2591, van der Wiel et al. (2011) propose this large-scale ‘plume’ could be CH₃OH liberated by a shock front. The outflow (Fig. 14 of van der Wiel et al. 2011) coincides spatially with the production of this ‘plume’. In addition, it is seen in a range of E_up transition lines and thus is not discriminatory by temperature. Our component two aligns spatially with this feature. Our component three, also at T_ex = 8 K, represents a blue-shifted structure to the northeast of VLA 3. Gieser et al. (2019), at small scales, suggest there is a young secondary outflow traced in SiO, and reason it is possible if there are numerous young stellar objects within VLA 3. Suri et al. (2021), subsequently, present NOEMA observations showing fragmentation of VLA 3 into three low-mass cores. In our observations, however, we do not have enough information to determine if CH₃OH component three is due to an associated outflow.

CH₃CCH has a range of E_up transition lines giving confidence in the LTE excitation temperature results. Similar to other star-forming regions, and as seen in WCCC sources, the carbon chains are on an extended scale with gas temperatures of around 30 K. In comparison, the COMs are not tracing the hot core at extended scales as such, and represent other mechanisms leading to their presence in the gas phase. Chemical modelling is needed to determine the formation routes, and is presented in Sect. 5.3.

The comprehensive large scale structure of AFGL 2591, as suggested by our results and the literature, is shown graphically in Fig. 13. The brown circles represent the sources VLA 1–5, with the dominant and most massive hot core VLA 3 highlighted. The warm, dense inner envelope is in light grey, while the methanol plume is in dark grey. The general direction of the outflows are represented in red and blue, marking the shift in velocity as well.

Fig. 13 Cartoon depiction of AFGL 2591. The brown circles represent the sources VLA 1–5, with the dominant protostellar source VLA 3 highlighted. The warm, dense inner envelope is in light grey, while the methanol plume is in dark grey. The general direction of the outflows are represented in red and blue, marking the shift in velocity as well.

5.2 IRAS 20126

IRAS 20126 is a relatively simple high-mass star-forming region. The peak of molecular emission, in both our integrated intensity (Fig. 6) and column density maps (Fig. 10), is located near the protostellar source. The central molecular gas peaks at the site of the protostar and is oriented perpendicu1ar to a small-scale outflow, involving a known disk (Cesaroni et al. 1997, 1999, 2005; Zhang et al. 1998). The increased line widths (Fig. 10) towards the central protostellar region can be hypothesised to be due to increased turbulence, as in Cesaroni et al. (1997). The found temperatures in this source are not high enough to produce thermal broadening – it would be, for a methanol molecule at 30 K, only 0.4 km s⁻¹.

The protostar is young and embedded in a dense hot core (Cesaroni et al. 1997), reaching up to 200 K in the core where COMs are known to exist (Cesaroni et al. 2005; Xu et al. 2011). Our single dish observations are insufficient to spatially resolve the protostar, and we are seeing the protostar smoothed into the surrounding dense gas. Cesaroni et al. (1997), also with IRAM 30m observations, found a rotational temperature of 50 K for CH₃OH, similar to the peak in our maps of both CH₃OH and CH₃CCH.

In our observations, CH₃CCH and c-C₃H₂ have visibly different velocity orientations than CH₃OH and H₂CO; however, neither group are oriented with the dense molecular gas surrounding the disk but with each of the two known outflows. The disk orientation, at scales of about 1000 AU (Sridharan et al. 2005), is also not seen in our beam of 16000 AU. A large-scale (up to 2′) outflow traced in ¹²CO has a red-shifted northern flow and a blue-shifted southern flow (Wilking et al. 1990; Shepherd et al. 2000), corresponding to the same gradient and orientation seen in the υ_1sr of CH₃CCH with velocities increasing southward. The small-scale (<30″) bipolar outflow oriented north-west to south-east is traced in simple molecules (Cesaroni et al. 1997, 1999), suggested to be a jet-driven bow shock (Su et al. 2007).

In Cesaroni et al. (1997), there are curious dynamics seen in the CS(3−2) and HCO+(1−0) transition lines with blue-shifted outer wings and red-shifted inner wings in the north-west, and red-shifted outer wings and blue-shifted inner wings in the south-east. They interpret the structure as due to the orientation of the outflow axis relative to the plane of the sky. This outflow is traced by CH₃OH and H₂CO in our observations with a blue-shifted south-east flow, indicating it arises similarly to the inner wings. Palau et al. (2017) also find that COMs are seen extended towards the south-east in IRAS 20126, suggesting they arise from the warm cavity walls or post-shock region. Palau et al. (2017), in their modelling, find that H₂CO, while possible to form in the gas phase, is co-spatial with larger COMs and thus likely arises from the same mechanism. Chemical modelling for our results is presented in Sect. 5.3.

As in AFGL 2591, at kpc distances our single dish observations are showing large-scale emission and not the hot core. While we cannot discern molecular presence at the hot core scale, and thus the differentiation between COMs and WCCC species, there is a notable difference in the physical conditions that the tracers CH₃OH and H₂CO arise in versus CH₃CCH and c-C₃H₂. These components, as seen in our observations and past studies, are summarised graphically in Fig. 14.

In the above discussion for both IRAS 20126 and AFGL 2591, our results are limited by the number of lines detected compared to the number of physical components present in the sources. A one- or two-component fit is representative in certain cases but may not reveal composite structure. Bouvier et al. (2020) are careful to note that a lot of complex molecule studies, which would include our own, are single dish observations and cannot spatially differentiate sub-structures of the regions.

These results would also have increased confidence with better S/N data, especially at high frequencies. There is a thorough E_up range, however, covered in the transition lines within our frequency ranges which lends confidence to the excitation temperatures found, especially when molecules are in LTE.

Fig. 14 Cartoon depiction of IRAS 20126. The orange circle represent the protostellar source, and the brown ellipse represents the central molecular gas and disk. The small-scale south-east-north-west outflow, traced by simple molecules, is in light red where red-shifted, and light blue where blue-shifted. The large-scale north-south outflow is similarly represented in dark red and dark blue.

5.3 Chemical model

In order to provide concrete support for the chemical formation of these molecules in the observed environments, we modelled the species abundances of both high-mass star-forming regions using the NAUTILUS gas-grain code (Ruaud et al. 2016). NAUTILUS simulates the chemical evolution of a region in three phases – gas, grain surface, and grain mantle – using coupled differential equations to describe a network of over 10 000 reactions connecting 489 species. With NAUTILUS, we model how the abundance of any given species changes over time as a function of different astrophysical conditions. These conditions can be static or dynamic throughout the simulation, and include properties such as gas density, dust and gas temperature, visual extinction, ultraviolet flux and cosmic ray ionisation.

In this discussion, we lead with and display the models for CH₃CCH and CH₃OH. These two species have more numerous strong detections relative to c-C₃H₂ and H₂CO and are representative of CCM and COM evolution, respectively.

5.3.1 AFGL 2591 chemical model results

CH₃CCH and CH₃OH

We first determined the fractional abundances of CH₃CCH and CH₃OH in the cloud – the ratio between the column density of the analysed chemical to the column density of molecular hydrogen, H₂. We calculated the H₂ column density map for AFGL 2591 using archival JCMT SCUBA2 850 μm emission (JCMTCAL observations on 2022/09/05). These maps were regridded in CASA⁷ to match the spatial resolution of the GBT and IRAM data.

The 850 μm continuum is caused by dust emission, which correlate to the H₂ column densities as molecular hydrogen generally forms on dust grains (Wakelam et al. 2017). We used a gas to dust mass ratio of 100, a constant dust mass opacity coefficient of κ =1 cm² g⁻¹ appropriate for MRN dust with a thin ice mantle (MRN dust describes the dust size distribution in the Galaxy; Mathis et al. 1977; Ossenkopf & Henning 1994), a mean mass per particle of μ = 2.32, and a temperature derived from the LTE modelling of CH₃CCH and CH₃OH (typically ≈30 K). We divided the CH₃CCH LTE model column density map by the H₂ column density map and found CH₃CCH abundances varying from 5 × 10⁻¹⁰ to 5 × 10⁻⁹ with the higher abundances in the northernmost ‘blob’ in Fig. 9.

CH₃OH is more complicated with three distinct kinematic components. As we could not, a priori, determine how much of the total H₂ column density is associated with each component, we explored the CH₃OH abundances in three small regions where we could reasonably isolate each component – a position in one component where there is no emission from the other two components. Since the LTE modelling of these three components utilised fixed temperatures of 18 K for component one and 8 K for components two and three, we used these temperatures to calculate the H₂ column densities in each of these components from the 850 μm map. We divided the CH₃OH column density maps for each component by the appropriate H₂ column density map and found CH₃OH abundances of 6.3 × 10⁻¹⁰ for component one (18 K), 3.7 × 10⁻¹⁰ for component two (8 K), and 2.3 × 10⁻¹⁰ for component three (8 K).

In NAUTILUS, we started with a two-phase model. First, a cold quiescent cloud was allowed to evolve for 10⁵ yr. This stage was run with initial abundances appropriate for the diffuse ISM (see Table 5), a temperature of 10 K, a density of 10⁴ cm⁻³, and a visual extinction of 50 magnitudes. Then, we assume that embedded protostars warm up the cloud, a stage evolved for another 10⁶ yr. In the second part of this simulation, the observed temperature of each pixel found from the LTE models is used as the dust and gas temperatures (typically 20–30 K), while the visual extinction is calculated from the average H₂ column density. Studies of the H₂ to A_V conversion range from (Rachford et al. 2009; Zhu et al. 2017) to (Bohlin et al. 1978; Whittet 1981). In this paper, we assumed an ‘average’ value of . 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. In following models, we modified both the CR ionisation rate and the UV scaling factor.

The volume densities of massive star-forming regions, on the scales that we present here, have been shown to have a centre-to-edge gradient of n ∝ r^−a where α ranges from 1.0 to 1.6 (van der Tak et al. 2000; Beuther et al. 2002). To account for such density gradients, we calculated models for three different densities: 10⁴ which is appropriate for the cloud edges, 5 × 10⁴ an intermediate value, and 10⁵ cm⁻³, which simulates conditions in the cloud cores.

Figure 15 shows the results of these connected simulations. The horizontal dashed blue and red lines show the minimum and maximum CH₃OH and CH₃CCH abundances, respectively. In general, this simple simulation can produce the observed CH₃CCH abundances for the majority of the cloud at about 10⁶ yr. This is true regardless of position – at the cloud edges where n = 10⁴ or the cores where n = 10⁵. The one caveat to this statement, however, occurs in the northern part of the map, where the CH₃CCH abundances are the highest. If the density in this region is as low as 10⁴ cm⁻³, then our models cannot account for these high observed abundances. Figure 15 also shows that we cannot reproduce the observed CH₃OH abundances at anytime for any density.

Given the poor match to our CH₃OH observations, we modified the chemical model to attempt to restrict the conversion of CH₃OH from the gas phase to the dust phase, or to help thermally desorb CH₃OH off the dust grains. A variety of methods were attempted, including increasing the ambient ultraviolet flux and cosmic ray ionisation rates. Increasing the scaling factor of the UV flux to S = 1000 had no effect upon the chemistry of either species. This is expected since the UV reaction rates take the form k ∝ e^{−S A}_v and A_V is generally so high that the UV field is completely attenuated. To simulate the effects of an external radiation field on the lower volume and column density cloud edges, we also ran a model with n = 10⁴ and A_V = 5. This simulation results in CH₃OH and CH₃CCH abundances that are respectively two and four orders of magnitude lower than those shown in Fig. 15 due to the increased rate of photodestruction of both molecules.

Increasing the cosmic ray ionisation rate by a factor of three increases the abundances of both species by a very small amount, but not enough to increase the CH₃OH abundances to observable values. Continued increases in the cosmic ray ionisation rate, however, have the opposite effect and result in lower abundances than those seen in Fig. 15, again due to an increased destruction rate.

As thermal desorption mechanisms cannot produce the observed CH₃OH abundances, we turn our attention to non-thermal ones. A shockwave passing through a cloud filled with dust may induce sputtering, which causes the dust grains to shatter (Jones et al. 1994), liberating the grain species back into the gas phase and causes a sharp increase in the gas phase abundances. An outflow, such as the ones seen in our sources, can produce such a shock (Zhang et al. 2005; Herbst & van Dishoeck 2009; Palau et al. 2017). Since it only includes thermal desorption processes, NAUTILUS cannot simulate sputtering directly. We simulated the effects of shock induced, non-thermal desorption (i.e. sputtering) by implementing a three stage model.

The first stage, cold cloud evolution, is the same as our prior model. The second stage, a shock, simulated the passage of a shock wave. We took the abundances at the end of stage one, removed all the grain surface species and added them back into the gas phase. The shock stage begun with enhanced gas phase abundances of all species and no grain surface species. During the shock we utilised the density and temperature evolution of the gas calculated by Palau et al. (2017) for IRAS 20126 (see their Fig. 5, top). This shock model follows the parametric approximation of Jiménez-Serra et al. (2008) for C-type shocks, assuming a shock velocity of 40 km s⁻¹ and a pre-shock density of 10⁴ cm⁻³. In our chemical models of this stage, we fixed the dust temperature to 80 K. The duration of the shock is 10⁴ yr, during which the temperature sharply increased to 2500 K before decreasing to 50 K, while the density slowly increased from 10⁴ to 8.2 × 10⁴ cm⁻³, where it stabilised for the rest of the simulation. In the post shock stage, we returned the gas to the temperature used in the warm-up stage above, and allowed the cloud to evolve for a final 10⁶ yr at a constant gas density of 8.2 × 10⁴ cm⁻³.

Figure 16 shows the results of the three stage chemical evolution model. The panel on the left shows the shock phase, and the panels on the right show the post shock evolution for temperatures of 8 K, appropriate for component two and component three of CH₃ OH, 18 K, appropriate for component one of CH₃OH, and 30 K, appropriate for CH₃CCH. As can be seen in the right panels, all of the observed abundances can be reached after about 10⁵ yr of postshock evolution. For CH₃CCH, the higher abundance limit is reached much earlier in the simulation during the shock phase. In AFGL 2591, while the majority of the cloud can be modelled as a simple cold cloud followed by a warm up phase, as is appropriate for WCCC, the northern region may require a non-thermal desorption mechanism for CH₃CCH if the density is around 10⁴ cm⁻³.

Fig. 15 Results from the NAUTILUS two stage chemical evolution model. The solid blue and red lines represent the time evolution of the CH3 OH and CH3CCH abundances, respectively. The blue and red dashed horizontal lines indicate the minimum and maximum observed abundances for CH3 OH and CH3 CCH, respectively, in AFGL 2591. The dotted black line indicates the minimum and maximum observed abundances in IRAS 20126. The left panel shows the cold cloud evolution stage and the right panels show the warm-up stage for three different densities.

H₂CO and c-C₃H₂

To model H₂CO and c-C₃H₂ we first assumed that the physical conditions were identical to those for CH₃CCH and CH₃OH; namely that the cloud experienced a 10 K cold stage followed by a 30 K warm-up stage. We also used the H₂ column density map calculated at 30 K to determine the H₂CO and c-C₃H₂ abundances. Under these conditions, the H₂CO abundance ranges from 3 × 10⁻⁹ in the cloud core to 10⁻⁸ in the north western cloud edges. The c-C₃H₂ abundances range from 4 × 10⁻¹¹ in the cloud core to 2 × 10⁻¹⁰ along the western and northern edges. Under these conditions, the observed c-C₃H₂ abundances are reached in the cold cloud stage and remain high well into the warm-up stage. Note that the same is true if we model the c-C₃H₂ abundances using a temperature of 20 K to calculate the H₂ column density. H₂CO, however, requires the 30 K warm-up stage in order to reach the observed abundances.

While an ‘average’ temperature of 30 K may be appropriate for c-C₃H₂ (see Fig. 9), it is not for H₂CO where the average temperature is about 40–50 K. For H₂CO we used temperatures of 50 K to calculate the H₂ column density and then the abundance map. Given the higher temperature, the resultant H₂ column densities are lower, resulting in slightly increased H₂CO abundances which now vary from 5 × 10⁻⁹ to 2 × 10⁻⁸. Under these conditions, the model cannot produce the observed abundances in the cold cloud stage, but can easily do so early in the 50 K warm-up stage. Therefore, no shock liberation of either molecule off the grain surfaces is required to produce the observed abundances.

5.3.2 IRAS 20126 chemical model results

CH₃CCH and CH₃OH

In IRAS 20126, we only considered the main component of the cloud since, other than the cold ‘foreground screen’ described in Sect. 4.2.3, the CH₃CCH and CH₃OH LTE model is dominated by one main component. To produce the corresponding H₂ column density map, we used archival JCMT SCUBA-1 850 μm continuum observations (Proposal ID M98BU24). These maps were re-gridded to match the spatial resolution of the GBT and IRAM data. We divided the CH₃CCH and CH₃OH LTE model column density maps by the H₂ column density map and found abundances varying from 10⁻¹⁰ and 10⁻⁹ for both species.

We followed the same chemical modelling approach in NAUTILUS as was used for AFGL 2591. We allowed the cloud to evolve for 10⁵ yr under cold conditions (n = 10⁴ cm⁻³, T = 10 K) starting with general ISM abundances, and then warmed the gas to 30 K to evolve for another 10⁶ yr.

Figure 15 shows the results of these simulations, along with the results for AFGL 2591. The left panel displays the cold cloud stage with the parameters listed above and the right panel shows the warm-up stage for a pixel with T = 30 K, A_V = 50 (i.e. N(H₂) = 10²³ cm⁻²), and, from the top, n = 10⁴ cm⁻³, n = 5 × 10⁴ cm⁻³, and n = 10⁵ cm⁻³. The dotted black horizontal lines indicate the observed abundance limits of 10⁻¹⁰ and 10⁻⁹. Figure 15 shows that although the observed CH₃CCH abundances can be reproduced after a few ×10⁵ yr in the warm up phase, the model CH₃OH abundances never reach the observed values. While we show the results for one specific pixel, these results are generally true in each pixel across the entire source. The main reason the modelled CH₃OH abundances are so low is due to freeze out of gas phase methanol onto the dust grains.

As discussed previously, a number of outflows and associated shocks have been identified and observed within IRAS 20126 (e.g. Shepherd et al. 2000; Zhang et al. 2001). Similarly to AFGL 2591, we simulated the effects of shock induced desorption with a three stage model. We included a shock stage in which all of the molecules frozen onto the dust (after the cold cloud evolution stage) are sputtered back into the gas phase, using the same gas density and temperature evolution. Then, after the shock stage, the cloud is allowed to continue to evolve at a density of 8.2 × 10⁴ cm⁻³ and at the currently observed temperatures derived from the LTE model.

Figure 16 presents the results from the shock on the left and post shock stages on the right, for our test pixel with a temperature of 30 K. The corresponding cold cloud stage is identical to Fig. 15. The dashed horizontal lines indicate the observed abundance limits of 10⁻¹⁰ and 10⁻⁹. Throughout the duration of the shock, and for the first 10⁵ yr of the post shock evolution, the CH₃OH abundance remains above the observed values. After this, CH₃OH begins to refreeze onto the dust grain surfaces and, in doing so, passes through the range of observed abundances. While CH₃CCH does not require a shock to produce the observed abundances (see Fig. 15), the passage of a shock can also produce the observed abundances after 10⁶ yr.

Again, while the results presented are for the physical parameters in a single pixel, these results are generally true across the entire region. Shocks, which we suggest are produced by the observed outflows, are required to explain the observed CH₃OH abundances over the entirety of IRAS 20126. They are not, however, needed to produce the observed abundances of CH₃CCH. This is as seen in AFGL 2591, with the exception of the northern region. Our results for CH₃CCH in IRAS 20126 are uniformly consistent with Warm Carbon Chain Chemistry.

Fig. 16 Results from the NAUTILUS three stage chemical evolution model. The solid blue and red lines represent the time evolution of the CH3OH and CH3CCH abundances, respectively. The blue and red dashed horizontal lines indicate the minimum and maximum observed abundances for CH3OH and CH3CCH, respectively, in AFGL 2591. The dotted black line indicates the minimum and maximum observed abundances in IRAS 20126. The left panel shows the shock stage and the right panels show the post shock stage for three different temperatures.

H₂CO and c-C₃H₂

The modelling of H₂CO and c-C₃H₂ proceeded in a fashion similar to that done for AFGL 2591. We first assumed the same conditions as for CH₃OH and CH₃CCH, and then used the temperatures determined from our LTE model (see Fig. 10). With the H₂ column density map calculated at 30 K we derived H₂CO abundance ranges from 6 × 10⁻¹¹ around the cloud edges to 3 × 10⁻¹⁰ in the cloud core. The c-C₃H₂ abundance ranges from 3 × 10⁻¹¹ in the cloud core and eastern edges to 10⁻¹⁰ along the western edges. Under these conditions, the abundances of both molecules are reached in the cold cloud stage and remain high well into the warm-up stage.

From Fig. 10 we see that the temperature derived for c-C₃H₂ is about 8 K and fairly uniform across the map. For H₂CO, the temperature is about 80 K and also fairly uniform. We used these two temperatures to calculate the respective H₂ column density maps and found that, for T = 8 K, the c-C₃H₂ abundance varies from 4 × 10⁻¹² in the cloud core and eastern edges to 10⁻¹¹ along the western edges. For T = 80 K, the H₂CO abundance ranges from 2 × 10⁻¹⁰ around the cloud edges to 10⁻⁹ in the cloud core. The chemical models are able to produce these observed abundances as well. For c-C₃H₂ at 8 K, the models match the observations early in the cold cloud stage. For H₂CO at 80 K, the model matches the observations early in the 80 K warm-up stage – although, the lowest abundances can be reached in the cold cloud stage as well. Therefore, as with AFGL 2591, no shock liberation of either molecule off the grain surfaces is required to produce the observed abundances.

6 Summary

This paper covers a targeted survey of carbon chain molecules in two high-mass star-forming regions - AFGL 2591 and IRAS 20126 - towards the Cygnus X star-forming complex. With observations from the IRAM 30 m telescope and the GBT 100 m dish, we present:

1. Detections of numerous rotational transition lines of two carbon chain molecules, CH₃CCH and c-C₃H₂, and two molecules related to the complex organic family, CH₃OH and H₂CO, in both sources;

2. An LTE model, developed in python to loop over several spectra in a map, that determines the physical environment the detected molecular species arise from;

3. Findings of carbon chain and complex organic molecules originating from different gases. The excitation temperature of carbon chains is typically 20–30 K while for complex organics can range from 8–85 K at these large scales (smoothing out the hot core). The velocity structures of each molecular type vary, and trace known structures in the regions;

4. Chemical evolution modelling with NAUTILUS that demonstrates that the observed abundances of CH₃CCH can be produced in the warm-up stage, where the embedded protostar has warmed the cloud to temperatures of 20–30 K. This is consistent with Warm Carbon Chain Chemistry. CH₃OH, on the other hand, requires a non-thermal desorption method to produce the observed abundances. A shock, sputtering the molecules off the dust grain surface, adequately does this. These mechanisms are found for both high-mass star-forming regions.

These results expand on current knowledge of carbon chain chemistry in high-mass star-forming regions. Warm carbon chain chemistry, and other envelope scale chemical mechanisms, can readily be studied with single dish telescopes. For an investigation of hot core chemistry, higher spatial resolution observations are required. These results also demonstrate the use of LTE models for mapping the physical environment using wide-band data cubes - the type of data many modern telescopes are capable of providing.

Acknowledgements

This work is based on observations carried out under project numbers 021-20 and 122-20 with the 30m 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. P.F. would like to thank Larry Morgan and Dave Frayer at the GBT/NRAO for their help developing observing scripts as well as GBTIDL scripts for ARGUS calibration and reduction. P.F. would also like to thank Pablo Torne and Monica Rodriguez at IRAM for their help observing with the IRAM telescope. P.F. and R.P. 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, P.F. and R.P. 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 Region 3.

Appendix A Additional LTE model description

Fig. A.1 Model spectrum, in colour, over the observed data, in black, of AFGL 2591 for transition lines of (a) CH3CCH on the brightest pixel of the CH3CCH 50−40 (Eup = 12 K), line, and (b) c-C3H2 on the brightest pixel of the c-C3H2 21,2−10,1 (Eup = 6 K), line. One component is used to fit the two carbon chain molecules together. The dashed grey line represents the observed transition frequency, with the energy level of this transition labelled.

Fig. A.2 Model spectrum, in colour, over the observed data, in black, of AFGL 2591 for transition lines of (a) CH3OH on the brightest pixel of the CH3OH on the brightest pixel of the CH3OH 60,6,0−50,5,0 (Eup = 48 K) line, and (b) H2CO on the brightest pixel of the H2CO 40,4−30,3 (Eup = 35 K) line. Three components are used to fit CH3OH while one component is used to fit H2CO. The dashed grey line represents the observed transition frequency, with the energy level of this transition labelled.

Fig. A.3 Parameters maps of AFGL 2951 for components two and three of CH3OH in AFGL 2591 from the LTE model. The black or white ‘x’ in each plot represents the source VLA 3. For both components, Tex is fixed at 8 K. For component two, the FWHM is fixed at 3.9 km s−1. For component three, the vlsr is fixed at −9.4 km s−1 and the FWHM at 2.5 km s−1

Fig. A.4 Same as Fig. A.1 for (a) CH3CCH and (b) c-C3H2 in IRAS 20126. One component is used to fit CH3CCH while two components are used for c-C3H2.

Fig. A.5 Same as Fig. A.2 for (a) CH3OH and (b) H2CO in IRAS 20126. Two components are used to fit CH3OH while one component is used to fit H2CO.

Fig. A.6 Parameters maps for component two of c-C3H2 and CH3OH in IRAS 20126 from the LTE model. The black or white ‘x’ in each plot represents the protostar. For c-C3H2, Ntot is fixed at 1 × 1013 cm−2 and the FWHM at 0.5 km s−1. For CH3OH, Ntot is fixed at 2.1 × 1013 cm−2, Tex at 10 K the FWHM at 1.2 km s−1.

References

All Tables

All Figures

	Fig. 1 Herschel PACS 160 μπι image of AFGL 2591 (proposal ID KPGT_fmotte_1, PI F. Motte). The colour scale is logarithmic from 0.8 to 80 Jy pixel−1, while the contours are 1, 5, 10, 30, 40, and 60 Jy pixel−1. The black ‘x’ represents the dominant protostellar source VLA 3.
In the text

	Fig. 2 Berschel PACS 160 μm image of IRAS 20126+4104 (proposal ID OT1_rcesaron_1, PI R. Cesaroni). The colour scale is logarithmic from 0.4 to 40 Jy pixel−1, while the contours are 1, 5, 10, 20, and 30 Jy pixel−1. The black ‘x’ represents the protostellar source.
In the text

Fig. 3 Select transition lines of the carbon chain molecules, CH3CCH and c-C3H2, the complex organic molecules CH3OH, and the precursor complex organic molecule H2CO, at the central pixel of the AFGL 2591 GBT and IRAM maps. The chosen lines represent a range of different upper energy levels, which are indicated in black on the plot. The dashed grey line represents the observed transition frequency at this location. The observed data is shown in black. The model spectra, to be discussed in Sect. 4.2 are shown in colour. The red solid line represents the total model, and the sole component if only one component is fit. For the multi-component fit in CH3OH the separate components are shown in blue, yellow, and brown dashed lines.

In the text

	Fig. 4 Same as Fig. 3 for the central pixel of the IRAS 20126 GBT and IRAM maps. For the multi-component fits in c-C3H2 and CH3OH the separate components are shown in blue and yellow dotted lines.
In the text

Fig. 5 AFGL 2591 integrated intensity, , maps of rotational transition lines spanning a range of energy levels. The upper state energy level is indicated in the top right-hand corner and the beam size in the bottom left-hand corner. The source VLA 3 is noted by a white or black ‘x’. CH3CCH, in the first column, c-C3H2, in the second, H2CO, in the third, have contour levels of 10, 20, 30, and 40 times the noise level of the respective frequency range. CH3OH, in the fourth column, has contour levels of 10, 20, 40, 80, and 120 times the noise level of the respective frequency range.

In the text

	Fig. 6 Same as Fig. 5 for IRAS 20126, the source noted by a white or black ‘x’. CH3CCH, in the first column, has contour levels of 10, 30, 50, and 70 times the noise of the respective frequency range. c-C3H2, in the second column, H2CO, in the third, and CH3OH, in the fourth, have contour levels of 10, 40, 70, and 100 times the noise of the respective frequency range.
In the text

	Fig. 7 Pixels fit with a RADEX model in AFGL 2591 marked on a Herschel PACS continuum image. The black and white numbers correspond to said pixels, with results listed in Table 3. The black ‘x’ represents the source VLA 3.
In the text

Fig. 8 RADEX fit for CH3OH in the central pixel of AFGL 2591, corresponding to pixel 7 in Fig. 7 and Table 3. The solid red line is the overall two-component model, the solid black line is the data, the dashed blue line is component one and the dashed yellow line is component two. Component two is only seen in the lower energy lines within the GBT range, thus component one and the overall model are the same for the lines in the IRAM range. The dashed grey line represents the observed transition frequency, with the energy level of this transition labelled.

In the text

	Fig. 9 Parameters maps of AFGL 2951 showing the results of the LTE models. The black ‘x’ represents the source VLA 3. For CH3OH, Tex was fixed at 18 K.
In the text

	Fig. 10 Parameters maps of IRAS 20126 showing the results of the LTE models. The black ‘x’ represents the protostellar source.
In the text

	Fig. 11 RADEX fit for c-C3H2 in the central pixel of IRAS 20126. The model used a temperature of 30 K, as expected from the CH3CCH results, an ortho- to para- ratio of 3, and a density of 105 cm−3. The model is in red while the data is in black. The upper state energy level of each transition is shown on each plot.
In the text

	Fig. 12 Column density ratio maps for CH3CCH, c-C3H2, and H2CO with respect to CH3OH. The upper panel shows AFGL 2591 and the lower panel IRAS 20126. The white ‘x’ represents the protostellar sources within each region.
In the text

	Fig. 13 Cartoon depiction of AFGL 2591. The brown circles represent the sources VLA 1–5, with the dominant protostellar source VLA 3 highlighted. The warm, dense inner envelope is in light grey, while the methanol plume is in dark grey. The general direction of the outflows are represented in red and blue, marking the shift in velocity as well.
In the text

	Fig. 14 Cartoon depiction of IRAS 20126. The orange circle represent the protostellar source, and the brown ellipse represents the central molecular gas and disk. The small-scale south-east-north-west outflow, traced by simple molecules, is in light red where red-shifted, and light blue where blue-shifted. The large-scale north-south outflow is similarly represented in dark red and dark blue.
In the text

Fig. 15 Results from the NAUTILUS two stage chemical evolution model. The solid blue and red lines represent the time evolution of the CH3 OH and CH3CCH abundances, respectively. The blue and red dashed horizontal lines indicate the minimum and maximum observed abundances for CH3 OH and CH3 CCH, respectively, in AFGL 2591. The dotted black line indicates the minimum and maximum observed abundances in IRAS 20126. The left panel shows the cold cloud evolution stage and the right panels show the warm-up stage for three different densities.

In the text

Fig. 16 Results from the NAUTILUS three stage chemical evolution model. The solid blue and red lines represent the time evolution of the CH3OH and CH3CCH abundances, respectively. The blue and red dashed horizontal lines indicate the minimum and maximum observed abundances for CH3OH and CH3CCH, respectively, in AFGL 2591. The dotted black line indicates the minimum and maximum observed abundances in IRAS 20126. The left panel shows the shock stage and the right panels show the post shock stage for three different temperatures.

In the text

Fig. A.1 Model spectrum, in colour, over the observed data, in black, of AFGL 2591 for transition lines of (a) CH3CCH on the brightest pixel of the CH3CCH 50−40 (Eup = 12 K), line, and (b) c-C3H2 on the brightest pixel of the c-C3H2 21,2−10,1 (Eup = 6 K), line. One component is used to fit the two carbon chain molecules together. The dashed grey line represents the observed transition frequency, with the energy level of this transition labelled.

In the text

Fig. A.2 Model spectrum, in colour, over the observed data, in black, of AFGL 2591 for transition lines of (a) CH3OH on the brightest pixel of the CH3OH on the brightest pixel of the CH3OH 60,6,0−50,5,0 (Eup = 48 K) line, and (b) H2CO on the brightest pixel of the H2CO 40,4−30,3 (Eup = 35 K) line. Three components are used to fit CH3OH while one component is used to fit H2CO. The dashed grey line represents the observed transition frequency, with the energy level of this transition labelled.

In the text

	Fig. A.3 Parameters maps of AFGL 2951 for components two and three of CH3OH in AFGL 2591 from the LTE model. The black or white ‘x’ in each plot represents the source VLA 3. For both components, Tex is fixed at 8 K. For component two, the FWHM is fixed at 3.9 km s−1. For component three, the vlsr is fixed at −9.4 km s−1 and the FWHM at 2.5 km s−1
In the text

	Fig. A.4 Same as Fig. A.1 for (a) CH3CCH and (b) c-C3H2 in IRAS 20126. One component is used to fit CH3CCH while two components are used for c-C3H2.
In the text

	Fig. A.5 Same as Fig. A.2 for (a) CH3OH and (b) H2CO in IRAS 20126. Two components are used to fit CH3OH while one component is used to fit H2CO.
In the text

	Fig. A.6 Parameters maps for component two of c-C3H2 and CH3OH in IRAS 20126 from the LTE model. The black or white ‘x’ in each plot represents the protostar. For c-C3H2, Ntot is fixed at 1 × 1013 cm−2 and the FWHM at 0.5 km s−1. For CH3OH, Ntot is fixed at 2.1 × 1013 cm−2, Tex at 10 K the FWHM at 1.2 km s−1.
In the text