11institutetext: University of Calgary, Calgary, Canada
11email<EMAIL_ADDRESS>
22institutetext: IRAP, Université de Toulouse, CNRS, UPS, CNES, Toulouse,
France
33institutetext: Tata Institute of Fundamental Research, Mumbai, India

# Parsec scales of carbon chain and complex organic molecules in AFGL 2591 and
IRAS 20126

P. Freeman 11 S. Bottinelli 22 R. Plume 11 E. Caux 22 C. Monaghan 11 B.
Mookerjea 33

(Submitted 22 June 2023; Revised 04 August 2023; Accepted 08 August 2023)

###### Abstract

Context. There is a diverse chemical inventory in protostellar regions leading
to the classification of extreme types of systems. Warm carbon chain chemistry
sources, for one, are the warm and dense regions near a protostar containing
unsaturated carbon chain molecules. Since the presentation of this definition
in 2008, there is a growing field to detect and characterise these sources.
The details are lesser known in relation to hot cores and in high-mass star-
forming regions – regions of great importance in galactic evolution.

Aims. To investigate the prevalence of carbon chain species and their
environment in high-mass star-forming regions, we have conducted targeted
spectral surveys of two sources in the direction of Cygnus X – AFGL 2591 and
IRAS 20126+4104.

Methods. We observed these sources in frequency ranges around 85, 96, and 290
GHz with the Green Bank Telescope and the IRAM 30m Telescope. We have
constructed a Local Thermodynamic Equilibrium (LTE) model using the observed
molecular spectra to determine the physical environment in which these
molecules originate. We map both the observed spatial distribution and the
physical parameters found from the LTE model. We also determine the formation
routes of these molecules in each source using the three-phase NAUTILUS
chemical evolution code.

Results. We detect several lines of propyne, CH3CCH, and cyclopropenylidene,
$c$-C3H2 as tracers of carbon chain chemistry, as well as several lines of
formaldehyde, H2CO, and methanol, CH3OH, as a precursor and a tracer of
complex organic molecule chemistry, respectively. We find excitation
temperatures of 20-30 K for the carbon chains and 8-85 K for the complex
organics. The observed abundances, used as input for the chemical evolution
code, are 10-9 to 10-10 for both CH3CCH and CH3OH. The CH3CCH abundances are
reproduced by a warm-up model, consistent with warm carbon chain chemistry,
while the observed CH3OH abundances require a shock mechanism sputtering the
molecules into the gas phase.

Conclusions. Single-dish observations are useful for studying the envelope-
scale chemistry of star-forming regions, including mechanisms such as warm
carbon chain chemistry. As well, LTE models lend well to the wide-band maps
obtained from these telescopes. The physical and chemical environment
determined for complex hydrocarbons and complex organics lends understanding
to high-mass star formation.

###### Key Words.:

stars: formation – ISM: molecules – Astrochemistry – Submillimeter: ISM

##  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 ($n_{\text{H}_{2}}$ $>$ 106
cm-3) 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, G. H. et al. 1996; Helmich & van Dishoeck 1997) as are their
compact ($\leq$ 100 AU in radius) equivalent in low-mass regions, dubbed hot
corinos (Cazaux et al. 2003; Ceccarelli 2004; Bottinelli et al. 2004, 2007).
In star-forming 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 CH3OH. 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. CnH, HCnN) that
are well-known in cold (T $\leq$ 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 CH4. As the
gas warms in the star formation process and CH4 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$\times$ 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/CH3OH 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 C4H (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 CH3OH, 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 CH4 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, CH3CCH, 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 high-mass 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).

CH3CCH 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,
CH3CCH develops from the successive hydrogenation of C3 (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 CH4, forming on then desorbing from grain
surfaces in the warm-up 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
CH3CCH.

Thaddeus et al. (1985) identified cyclopropenylidene, $c$-C3H2, 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$-C3H2 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 CH4, 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-C3H2 correlate
with CH4.

###  2.2 Complex organic molecules – formaldehyde and methanol

Formaldehyde, H2CO, 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. H2CO was the first polyatomic organic molecule
detected in space, Snyder et al. (1969) observed absorption spectra towards
both galactic and extragalactic sources. Methanol, CH3OH, 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.

H2CO and CH3OH 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 H2CO
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 CH3OH transitions are masers, pumped into an excited state either
collisionally (Class I) or by far-infrared radiation (Class II). CH3OH 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).

###  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 $\mu m$ continuum image of the source. At a distance of 3.33$\pm$0.11 kpc
(Rygl et al. 2012), AFGL 2591 has a total mass of about 2$\times 10^{4}$ M⊙
and a measured total IR luminosity of 2$\times 10^{5}$ 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$\times 10^{4}$ 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 (CH3OH, H2CO, 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.

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

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 $\mu 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$\times 10^{4}$ 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 C34S 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 $\sim$ 1000 AU). This was confirmed by Chen et al. (2016)
whose modelling of SMA observations of CH3CN 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, CH3CCH). 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.

Figure 2: Herschel PACS 160 $\mu 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, while the contours are 1, 5, 10, 20, and 30 Jy/pixel. 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^{\prime}\times 1.5^{\prime}$ region in both AFGL 2591 and IRAS 20126 in the
frequency range 287.22 - 295 GHz. The offset position was -120$\arcsec$
horizontally and -240$\arcsec$ vertically. The phase centre was
$\alpha$J(2000)$=20^{\mathrm{h}}29^{\mathrm{m}}24\aas@@fstack{s}9$ and
$\delta$J(2000)$=+40\degr 11\arcmin 20\aas@@fstack{\prime\prime}0$ for AFGL
2591 and $\alpha$J(2000)$=20^{\mathrm{h}}14^{\mathrm{m}}26\aas@@fstack{s}04$
and $\delta$J(2000)$=+41\degr 13\arcmin 32\aas@@fstack{\prime\prime}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 minutes or less, was done using either K3-50A or NGC 7027 and
resulted in typical corrections of 1-2′′. The focus, completed every three
hours or after a sunset or sunrise, was done with the same source as for
pointing and resulted in typical corrections of $<$ 0.3$\arcsec$. 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$\arcsec$ at 290 GHz, with a pixel size of
4.4$\arcsec$. Data reduction, and the production of maps, was completed with
GILDAS/CLASS111https://www.iram.fr/IRAMFR/GILDAS.

The 100m 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$\arcmin$ 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$\arcmin$ in azimuth for the reference measurements. Each petal,
with a map radius of 0.5$\arcmin$, took 54.7 seconds with a scan rate of
2.3$\arcsec$ per second and an integration time of 1 second. The pointing and
focus was completed every 30-40 minutes 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-200 K. With spectrometer mode 2 we obtained a spectral sampling of 92 kHz
in both frequency ranges, or 0.32 km s-1 at 85 GHz and 0.28 km s-1 at 96 GHz.
The beam size is 9.2$\arcsec$ or 10.0$\arcsec$ at 96 or 85 GHz respectively,
with a pixel size of 2$\arcsec$. The data were reduced and calibrated using
GBTIDL222https://gbtidl.nrao.edu/index.shtml.

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-1. 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 CASA333https://casa.nrao.edu/ to match the 85 GHz
beam size and spatial grid. The final restored beam is 10.0$\arcsec$ with a
pixel size of 4$\arcsec$. Average rms noise levels are shown in Table 1. In
the Results, Sect. 4, the data are in units of TA*, adjusted during
calibration. In the LTE Model, the data are converted to Tmb using the
telescope Beff/Feff 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.

Table 1: Details of the data set used in this paper.

Telescope | Frequency | Spectral Sampling | Spatial Sampling | IRAS 20126 RMS noise | AFGL 2591 RMS noise
---|---|---|---|---|---
| [GHz] | [km s-1] | [$\arcsec$] | [TA*, mK] | [TA*, mK]
GBT | 84.50-85.75 | 0.32 | 4 | 27.20 | 25.91
GBT | 95.55-96.80 | 0.28 | 4 | 25.99 | 23.61
IRAM | 287.00-295.00 | 0.80 | 4 | 37.49 | 75.28

444The spectral sampling is the channel width and the spatial sampling is the
pixel size.

##  4 Results

###  4.1 Detected lines

Species were identified through the CASSIS555http://cassis.irap.omp.eu/
LTE+RADEX666http://cassis.irap.omp.eu/help/?page=html/m_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 $v_{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 CH3CCH (CDMS tag 40502 for the
vibrational ground state), $c$-C3H2 (CDMS tag 38508), CH3OH (CDMS tag 32504),
and H2CO (CDMS tag 30501) lines with Aij $>$ 1$\times 10^{-6}$ s-1 and Eup $<$
150 K. For CH3OH the Aij limit was set to 2$\times 10^{-6}$ s-1 to remove the
40 K maser line. For CH3CCH the Eup limit was increased to 200 K to account
for the a-CH3CCH line that has Eup $<$ 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-CH3CCH, o- and
p-$c$-C3H2, A- and E-CH3OH, o- and p-H2CO – 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 Eup a-CH3CCH 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
CH3CCH 50-40 (Eup = 12 K) line and the CH3OH 20-10, 2-1-1-1, 20-10 triplet
(Eup = 7, 13, 20 K). Figures 17, 18, 20, and 21 display all transition lines
of each molecule. The propyne detections with Eup 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.

Figure 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. Figure 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. Table 2: Rotational transition lines of
selected molecules in the observed frequency ranges.

Molecule | Quantum Numbers | Frequency | $E_{\rm up}$ | $A_{\rm ij}$ | Quantum Numbers | Frequency | $E_{\rm up}$ | $A_{\rm ij}$
---|---|---|---|---|---|---|---|---
|  | [GHz] | [K] | [10-5 s-1] |  | [GHz] | [K] | [10-5 s-1]
CH3CCH | 50-$4_{0}$ | e | 85.4573 | 12.31 | 0.2 | 170-170 | e | 290.5020 | 125.50 | 8.5
| 51-41 | e | 85.4557 | 19.53 | 0.2 | 171-171 | e | 290.4965 | 132.73 | 8.5
| 52-42 | e | 85.4508 | 41.21 | 0.2 | 172-172 | e | 290.4799 | 154.40 | 8.4
| 53-43 | a | 85.4426 | 77.34 | 0.1 | 173-173 | a | 290.4522 | 190.52 | 8.3
$c$-C3H2 | 21,2-10,1 | o | 85.3389 | 6.45 | 2.3 | 43,2-42,3 | o | 85.6564 | 29.07 | 1.5
| 32,2-31,3 | p | 84.7277 | 16.14 | 1.0 | 75,3-66,0 | p | 294.0356 | 82.91 | 2.4
CH3OH | 20,2,0-10,1,0 | A+ | 96.7413 | 6.96 | 0.3 | 43,2,1-52,4,1 | E | 288.7056 | 70.93 | 1.8
| 21,2,2-11,1,2 | E | 96.7394 | 12.54 | 0.3 | 62,5,1-52,4,1 | E | 290.3077 | 71.01 | 9.3
| 20,2,1-10,1,1 | E | 96.7445 | 20.09 | 0.3 | 62,4,2-52,3,2 | E | 290.3073 | 74.66 | 9.5
| 21,2,0-11,1,0 | A+ | 95.9143 | 21.45 | 0.2 | 62,5,0-52,4,0 | A- | 290.1847 | 86.46 | 9.5
| 21,1,1-11,0,1 | E | 96.7555 | 28.01 | 0.3 | 62,4,0-52,3,0 | A+ | 290.2641 | 86.47 | 9.5
| 51,5,2-40,4,1 | E | 84.5212 | 40.39 | 0.2 | 63,4,1-53,3,1 | E | 290.2132 | 96.47 | 8.0
| 60,6,0-50,5,0 | A+ | 290.1106 | 48.74 | 10.6 | 63,4,0-53,3,0 | A+ | 290.1895 | 98.55 | 7.9
| 32,1,0-41,4,0 | A+ | 293.4641 | 51.64 | 2.9 | 63,3,0-53,2,0 | A- | 290.1905 | 98.55 | 7.9
| 61,6,2-51,5,2 | E | 290.0697 | 54.32 | 10.3 | 63,3,2-53,2,2 | E | 290.2097 | 111.47 | 8.0
| 60,6,1-50,5,1 | E | 289.9394 | 61.86 | 10.6 | 64,3,0-54,2,0 | A- | 290.1613 | 129.10 | 5.9
| 61,6,0-51,5,0 | A+ | 287.6708 | 62.87 | 10.1 | 64,3,0-54,2,0 | A+ | 290.1614 | 129.10 | 5.9
| 61,5,0-51,4,0 | A- | 292.6729 | 63.71 | 10.6 | 64,3,2-54,2,2 | E | 290.1624 | 136.66 | 5.9
| 61,5,1-51,4,1 | E | 290.2487 | 69.81 | 10.6 | 64,2,1-54,1,1 | E | 290.1833 | 144.76 | 5.9
H2CO | 40,4-30,3 | p | 290.6234 | 34.90 | 69 | 43,2-33,1 | o | 291.3804 | 140.94 | 30.4
| 42,3-32,2 | p | 291.2378 | 82.07 | 52.1 | 43,1-33,0 | o | 291.3844 | 140.94 | 30.4
| 42,2-32,1 | p | 291.9481 | 82.12 | 52.5 |  |  |  |  |

777The line properties are taken from the CDMS and VASTEL catalogues. Limits
were set for Eup $<$ 150 K (200 K for CH3CCH) and Aij $>$ 1.0$\times 10^{-6}$
(2.0 $\times 10^{-6}$ for CH3OH).

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 CH3OH, the warm lines of Eup =
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 CH3OH. 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 CH3CCH and $c$-C3H2, respectively.

The beam size, 10.0′′, as displayed in the integrated intensity maps,
corresponds to 3.3$\times$104 AU (0.16 pc) in AFGL 2591 and 1.6$\times$104 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).

Figure 5: AFGL 2591 integrated intensity $\int T_{A}\mathrm{d}v$ (K km s-1)
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. Figure 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 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
software888http://cassis.irap.omp.eu/docs/RadiativeTransfer.pdf. 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 Levenberg-Marquardt 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 $\chi^{2}$, 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 CH3OH 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 CH3OH 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 et al. 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 CH3CCH 12 K 50-40 (Eup = 12 K) and 51-41 (Eup =
19 K) lines at 85 GHz are included in the CH3CCH and $c$-C3H2 LTE models, and
any pixel that had the same threshold for the CH3OH 20,2,0-10,1,0 (Eup = 7 K)
and 21,2,2-11,1,2 (Eup = 12 K) CH3OH lines at 96 GHz or the 60,6,0-50,5,0 (Eup
= 48 K) and 61,6,2-51,5,2 (Eup = 54 K) lines at 290 GHz (to capture lines in
both GBT and IRAM data) are used for the CH3OH and H2CO 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$-C3H2
and H2CO masks produced more extended maps. To allow comparison between
molecules we continued with the more restrictive CH3CCH and CH3OH maps with a
loosened S/N threshold.

To begin the analysis, we extracted the spectrum from the brightest pixel of
certain molecular transition lines – CH3CCH 50-40 (Eup = 12 K); $c$-C3H2
21,2-10,1 (Eup = 6 K); CH3OH 20,2,0-10,1,0, 21,2,2-11,1,2, 20,2,1-10,1,1
triplet (Eup = 7, 12, 20 K), CH3OH 60,6,0-50,5,0 (Eup = 48 K); H2CO 40,4-30,3
(Eup = 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_{\rm up,max}$ and $A_{\rm 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$\arcsec$, 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$\arcsec$) beam sizes,
the resulting column density changed but remained within a factor of 3. With
larger (30$\arcsec$) 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.

####  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 6.

The initial one-component fit of $c$-C3H2 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 (Eup = 6 K)
$c$-C3H2 transition line which dominates the fits. Given this, we were unable
to model with $c$-C3H2 alone and it was modelled in the same component as
CH3CCH in the CH3CCH 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 CH3CCH, $c$-C3H2 reproduces the same structure
of the integrated intensity maps of $c$-C3H2 in the column density maps,
indicating that the column density reflects where we see emission in the data.
Previous models, with $c$-C3H2 alone, produced unrealistic column density
distributions. In addition, the model fits the spectral lines for each
molecule well (example in Fig. 20), 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$-C3H2 could also show some non-LTE behaviour.

The CH3OH spectral line profiles contain more velocity structure than compared
to the other molecules. CH3OH was fit unsuccessfully by a one-component fit,
as determined by non-converging models, high reduced $\chi^{2}$ 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 CH3OH 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 (1014 cm-2), warmer (20 K), gas near the source vlsr of $-$5.7
km s-1 represented by the higher energy IRAM lines, and a second which traces
the colder (10 K), lower column density (1013 cm-2) more negative velocity
($-$7 to $-$9 km s-1) 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 CH3OH. In
subsequent LTE modelling, the latter cold component was able to be
differentiated into two components in the LTE model – the CH3OH ‘blob’ towards
the north-east below in component two and the blue-shifted component three,
both described below.

H2CO, 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 H2CO and CH3OH in the same component of the model to represent
these molecules originating from the same gas. However, none of the multiple
colder components of CH3OH aligned with the hotter H2CO 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 CH3OH
are displayed in Fig. 19.

All four molecules show a fairly constant velocity structure near the source
vlsr of $-$5.7 ($\pm$0.03-0.5) km s-1. Component three of CH3OH, however, is
notably blue-shifted and fixed at $-$9.4 km s-1. Component two of CH3OH,
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 ($\pm$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$-C3H2 lines
dominate. The temperature of H2CO reaches towards 70 ($\pm$5) K slightly
offset from the VLA 1-3 sources and shows a warm temperature across extended
spatial scales. For CH3OH, 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. CH3CCH, at a peak of 4 ($\pm$0.4) $\times$1014 cm-2,
extends along and towards the east of AFGL 2591, with two distinct peaks.
$c$-C3H2 has the lowest column density, 2 ($\pm$0.3) $\times$1013 cm-2, and
shows an extended distribution with a slight peak near the VLA sources. H2CO
peaks at 9($\pm$0.03) $\times$1013 cm-2, and distinctly is distributed near
the VLA sources. CH3OH peaks at 5 ($\pm$0.2), 5 ($\pm$0.2), and 2 ($\pm$0.1)
$\times$1014 cm-2 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-1 wide across
the emitting region, and H2CO is uniform at a slightly broader width near 4 km
s-1. Components two and three of CH3OH were fixed at 3.9 km s-1 and 2.5 km s-1
respectively, while component one results in line widths of 2-3 km s-1.

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.

Figure 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. Table 3:
RADEX results for CH3OH in AFGL 2591.

Pixel | Tkin | Vlsr | Ntot | NH2 | FWHM
---|---|---|---|---|---
| (K) | (km s-1) | ($\mathrm{cm^{-2}}$) | ($\mathrm{cm^{-3}}$) | (km s-1)
1 (32,34) | 12.1 | -6.15 | 4.1$\times 10^{14}$ | 6.4$\times 10^{7}$ | 2.8
| 11.0 | -7.30 | 12.9$\times 10^{13}$ | 1.0$\times 10^{5}$ | 3.2
2 (35,34) | 17.5 | -6.35 | 4.1$\times 10^{14}$ | 1.5$\times 10^{7}$ | 3.1
| 10.6 | -9.50 | 6.7$\times 10^{13}$ | 1.0$\times 10^{5}$ | 3.2
3 (38,34) | 19.5 | -5.95 | 3.5$\times 10^{14}$ | 3.9$\times 10^{8}$ | 2.8
| 11.2 | -9.50 | 10.0$\times 10^{13}$ | 1.0$\times 10^{5}$ | 3.2
4 (41,34) | 16.5 | -6.00 | 2.3$\times 10^{14}$ | 1.0$\times 10^{7}$ | 1.9
| 10.9 | -9.30 | 3.2$\times 10^{13}$ | 1.0$\times 10^{5}$ | 3.2
5 (32,31) | 14.5 | -5.30 | 2.4$\times 10^{14}$ | 3.8$\times 10^{6}$ | 3.1
| 10.2 | -7.10 | 3.1$\times 10^{13}$ | 1.0$\times 10^{5}$ | 3.2
6 (35,31) | 16.5 | -5.70 | 3.6$\times 10^{14}$ | 7.8$\times 10^{6}$ | 2.6
| 10.8 | -8.90 | 14.9$\times 10^{13}$ | 1.0$\times 10^{5}$ | 3.2
7 (38,31) | 21.4 | -5.75 | 2.6$\times 10^{14}$ | 3.8$\times 10^{7}$ | 2.6
| 10.3 | -9.60 | 2.6$\times 10^{13}$ | 1.0$\times 10^{5}$ | 3.2
8 (41,31) | 19.1 | -6.15 | 2.0$6\times 10^{14}$ | 1.1$\times 10^{7}$ | 2.4
| 10.3 | -9.00 | 3.4$\times 10^{13}$ | 1.5$\times 10^{5}$ | 3.2
9 (39,42) | 15.3 | -4.85 | 6.6$\times 10^{13}$ | 1.9$\times 10^{6}$ | 3.2
| 11.6 | -8.8 | 19.2$\times 10^{13}$ | 1.0$\times 10^{5}$ | 3.2
999Eight pixels in E-W strips across the central source of AFGL 2591, and one
off-source, were taken and fit with RADEX, a non-LTE model. These are labelled
in black and white in Fig. 7. Two components were fit to each pixel’s spectra.
A- and E-CH3OH are used in this model due to the availability of collisional
co-efficients. There are not enough lines to determine the A-/E- ratio, thus
it is assumed to be 1.

Figure 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. Figure 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.

####  4.2.3 IRAS 20126 LTE model results

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

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$-C3H2 and CH3CCH
combination did not represent the lines well. The $c$-C3H2 lines in IRAS 20126
are stronger, and drove the excitation temperatures too low for CH3CCH. CH3OH
and H2CO 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
CH3CCH and H2CO, one component was adequate. For $c$-C3H2 and CH3OH, 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$-C3H2 21,2-10,1 (Eup = 6 K)
and CH3OH 20,2,0-10,1,0 (Eup = 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$-C3H2 and CH3OH are displayed in Fig. 22.

There is a consistent velocity range for all species between $-$2.0 to $-$4.5
($\pm$0.02-0.1) km s-1. This is a range around the source Vlsr of $-$3.5 km
s-1. 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$-C3H2, the second
component was fixed at the same velocity as the first component, and in CH3OH
the second component was fixed at $-$0.2 km s-1 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 CH3CCH at 25-30
($\pm$0.9-4) K and $c$-C3H2 at 8-12 ($\pm$0.8-3) K. The COMs show a warmer
central region falling off in temperature towards the edges, with a range of
15-35 ($\pm$0.3-2) K for CH3OH located just off the protostar, and 55-85
($\pm$3-6) K for H2CO located on the protostar. The second component of
$c$-C3H2 varies from 4-8 ($\pm$0.4-6) K with few pixels at 10-12 K, while the
second component of CH3OH was fixed at 10 K.

Given the low excitation temperature of $c$-C3H2, 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$-C3H2, as shown in Fig. 11. We reproduced the spectrum using 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.

All column density maps peak near the site of the protostar, with the COMs
showing a steeper falloff towards the edges. CH3CCH, CH3OH, and H2CO peak near
6, 9, and 3 ($\pm$ 0.3, 0.2, and 0.1) $\times$1014 cm-2 respectively, while
$c$-C3H2 peaks near 5 ($\pm$0.3) $\times$1013 cm-2. The second component of
$c$-C3H2 was fixed at 1$\times$1013 cm-2 and for CH3OH it was fixed at
2.1$\times$1013 cm-2.

The CCM lines are consistently narrower in line width, ranging from about
1-2.5 km s-1, while the COM lines are 4-7km s-1 wide. The former show a
slightly wider line profile near the protostar, while the latter show a
widening towards the south-east of the source. The second component of
$c$-C3H2 was fixed at a line width of 0.5 km s-1 while for CH3OH it was fixed
at 1.2 km s-1.

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.

Figure 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 C3 in DR21(OH), another source with hot cores in Cygnus X. The C3
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 CH3OH is a recently available and
practised method of analysing star-forming regions through molecular
observations. Table 4 summarises several studies that used either CH3CCH or
$c$-C3H2, as in this study. At first glance, there is a spread of
N(CH3CCH)/N(CH3OH), possibly due to the variety of studies and sources that
examine this ratio. All CH3CCH studies presented used single dish observations
and resolved objects on scales of 1-7$\times$104 AU, comparable to this study,
and consistent with this species being found in the cold and warm extended
envelopes. In Fayolle et al. (2015), CH3CCH 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(CH3CCH)/N(CH3OH) $>$ 1, are potentially younger than
typical hot core sources, while Taniguchi et al. (2018b) find
N(CH3CCH)/N(CH3OH) $>$ 1 for their source with a UCH II region. For $c$-C3H2,
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(CH3CCH)/N(CH3OH) = 0.99
while N($c$-C3H2)/N(CH3OH) = 0.06. For AFGL 2591, N(CH3CCH)/N(CH3OH) = 0.97
while N($c$-C3H2)/N(CH3OH) = 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.

Table 4: Column density ratios from the literature of carbon chain molecules
with respect to methanol.

Source | Source Type | N(CH3CCH)/N(CH3OH) | N($c$-C3H2)/N(CH3OH)
---|---|---|---
G331.51-0.103a | hot molecular core | 0.42$\pm$0.05 | -
G12.89+0.49b | massive young stellar object | 0.34${}^{+0.28}_{-0.15}$ | -
G16.86-2.16b | massive young stellar object | 0.36${}^{+0.28}_{-0.18}$ | -
G28.28-0.36b | massive young stellar object | 1.61${}^{+1.4}_{-0.86}$ | -
Inner galaxy c | ATLASGAL high mass clumps, averaged | 0.31 | -
NGC 7538 IRS 9d | organic-poor massive young stellar object | 1.3$\pm$0.4 | -
W3 IRS5d | organic-poor massive young stellar object | 2.2$\pm$0.7 | -
AFGL 490d | organic-poor massive young stellar object | 1.8$\pm$0.8 | -
Perseus molecular cloud e | Class 0/1 protostars | - | 0.009$\pm$0.003-0.80$\pm$0.34
L1527e | low-mass star-forming region | - | 0.60$\pm$0.07
IRAS 20126f | high-mass star-forming region | 0.99 | 0.06
AFGL 2591f | high-mass star-forming region | 0.97 | 0.05

aSantos et al. (2022), bTaniguchi et al. (2018b), cGiannetti et al. (2017), as
averaged in a, dFayolle et al. (2015), eHiguchi et al. (2018) and references
therein, fthis work.

Figure 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.

Environmental differences cause the variation between CCMs and CH3OH. In our
data, the column densities of CH3CCH and CH3OH are of comparable magnitude,
and CH3CCH is typically seen on a uniform, extended scale in comparison to
concentrated CH3OH closer to protostellar sources. Santos et al. (2022) most
recently compare N(CH3CCH)/N(CH3OH) 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). CH3OH is likely to be more
spatially compact than CH3CCH, 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 CH3OH rich
hot cores, though the presence and relative abundance of CH3CCH and other
complex species makes them comparable to, or possible precursors to, high mass
hot cores.

$c$-C3H2 is often an order of magnitude lower in column density than CH3OH,
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$-C3H2 and CH3OH 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$-C3H2 and CH3OH for 36 protostars in the Perseus
region. The column density ratio they find is on order of magnitude 10-2-10-1,
similar to what we find in both IRAS 20126 and AFGL 2591.

###  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+, H2CO (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 CH3OH 20,2,0-10,1,0 (Eup = 7 K) line and component
two CH3OH 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${}_{\mathrm{up}}\leq$
200 K (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\arcmin$) 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
CH3OH 20,2,0-10,1,0 (Eup = 7 K) (Fig. 5) line, which is blended with the
21,2,2-11,1,2 and 20,2,1-10,1,1 (Eup = 13, 20 K) lines, clearly traces this
plume in multiple energy levels and at greater intensity than around VLA 1-5.
The 60,6,0-50,5,0 (Eup = 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 CH3OH 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 CH3OH 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 CH3OH detected. Soma et al. (2015),
similarly in Taurus Molecular Cloud 1, offer several possibilities for non-
thermal CH3OH 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 CH3OH
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 Eup transition lines and thus is not discriminatory by
temperature. Our component two aligns spatially with this feature. Our
component three, also at Tex = 8 K, represents a blue-shifted structure to the
north-east 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 CH3OH component three is due to an associated outflow.

CH3CCH has a range of Eup 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.

Figure 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 perpendicular
to a small-scale outflow, involving a known disk (Cesaroni et al. 1997; Zhang
et al. 1998; Cesaroni et al. 1999, 2005). 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 -1.

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
CH3OH, similar to the peak in our maps of both CH3OH and CH3CCH.

In our observations, CH3CCH and $c$-C3H2 have visibly different velocity
orientations than CH3OH and H2CO; 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 16 000 AU. A large-scale (up to
2$\arcmin$) outflow traced in 12CO 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 $v_{lsr}$ of
CH3CCH with velocities increasing southward. The small-scale ($<$ 30$\arcsec$)
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 CH3OH and H2CO 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 H2CO, 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 CH3OH and H2CO arise in versus CH3CCH and $c$-C3H2. 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 Eup 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.

Figure 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 CH3CCH and CH3OH.
These two species have more numerous strong detections relative to $c$-C3H2
and H2CO and are representative of CCM and COM evolution, respectively.

Table 5: Initial ISM abundances in the NAUTILUS model. Species | Abundance
---|---
He | 9.0 $\times$ 10-2
N | 6.2 $\times$ 10-5
O | 2.4 $\times$ 10-4
H | 0.0
H2 | 0.5
C+ | 1.7 $\times$ 10-4
S+ | 1.5 $\times$ 10-5
Si+ | 8.0 $\times$ 10-9
Fe+ | 3.0 $\times$ 10-9
Na+ | 2.0 $\times$ 10-9
Mg+ | 7.0 $\times$ 10-9
P+ | 2.0 $\times$ 10-10
Cl+ | 1.0 $\times$ 10-9
F | 6.7 $\times$ 10-9

####  5.3.1 AFGL 2591 chemical model results

#### CH3CCH and CH3OH

We first determined the fractional abundances of CH3CCH and CH3OH in the cloud
– the ratio between the column density of the analysed chemical to the column
density of molecular hydrogen, H2. We calculated the H2 column density map for
AFGL 2591 using archival JCMT SCUBA2 850 $\mu m$ emission (JCMTCAL
observations on 2022/09/05). These maps were re-gridded in
CASA101010https://casa.nrao.edu/ to match the spatial resolution of the GBT
and IRAM data.

The 850 $\mu m$ continuum is caused by dust emission, which correlate to the
H2 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 $\kappa$ = 1 cm2/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 $\mu=2.32$, and a temperature derived from the LTE modelling of
CH3CCH and CH3OH (typically $\approx$ 30 K). We divided the CH3CCH LTE model
column density map by the H2 column density map and found CH3CCH abundances
varying from $5\times 10^{-10}$ to $5\times 10^{-9}$ with the higher
abundances in the northernmost ‘blob’ in Fig. 9.

CH3OH is more complicated with three distinct kinematic components. As we
could not, a priori, determine how much of the total H2 column density is
associated with each component, we explored the CH3OH 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 H2 column densities in each of these components
from the 850 $\mu m$ map. We divided the CH3OH column density maps for each
component by the appropriate H2 column density map and found CH3OH abundances
of $6.3\times 10^{-10}$ for component one (18 K), $3.7\times 10^{-10}$ for
component two (8 K), and $2.3\times 10^{-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 105 years. This stage was run with initial
abundances appropriate for the diffuse ISM (see Table 5), a temperature of 10
K, a density of $10^{4}$ cm-3, and a visual extinction of 50 magnitudes. Then,
we assume that embedded protostars warm up the cloud, a stage evolved for
another 106 years. 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 H2 column density. Studies of the H2 to AV
conversion range from AV = $\frac{N_{H_{2}}}{2.1\times 10^{21}}$ (Rachford et
al. 2009; Zhu et al. 2017) to AV = $\frac{N_{H_{2}}}{1.9\times 10^{21}}$
(Bohlin et al. 1978; Whittet 1981). In this paper, we assumed an ‘average’
value of AV = $\frac{N_{H_{2}}}{2\times 10^{21}}$. We used a standard cosmic
ray ionisation rate of $1.3\times 10^{-17}$ s-1. For the UV field, the
NAUTILUS code provided a function of the form $S\times 10^{8}$ photons cm-2
s-1 (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\propto
r^{-\alpha}$ where $\alpha$ 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: 104 which is appropriate for the cloud
edges, 5 $\times$ 104 an intermediate value, and 105 cm-3 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 CH3OH and CH3CCH
abundances, respectively. In general, this simple simulation can produce the
observed CH3CCH abundances for the majority of the cloud at about $10^{6}$
years. This is true regardless of position – at the cloud edges where n = 104
or the cores where n = 105. The one caveat to this statement, however, occurs
in the northern part of the map, where the CH3CCH abundances are the highest.
If the density in this region is as low as $10^{4}$ cm-3, then our models
cannot account for these high observed abundances. Figure 15 also shows that
we cannot reproduce the observed CH3OH abundances at anytime for any density.

Given the poor match to our CH3OH observations, we modified the chemical model
to attempt to restrict the conversion of CH3OH from the gas phase to the dust
phase, or to help thermally desorb CH3OH 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\propto e^{-S\leavevmode\nobreak\
A_{V}}$ and AV 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 = 104
and AV = 5. This simulation results in CH3OH and CH3CCH 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 CH3OH 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 CH3OH 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-1 and a pre-shock density of
$10^{4}$ cm-3. In our chemical models of this stage, we fixed the dust
temperature to 80 K. The duration of the shock is $10^{4}$ years, during which
the temperature sharply increased to 2500 K before decreasing to 50 K, while
the density slowly increased from $10^{4}$ to $8.2\times 10^{4}$ cm-3, 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^{6}$ years at a constant gas
density of $8.2\times 10^{4}$ cm-3.

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 CH3OH, 18 K, appropriate for component one of CH3OH,
and 30 K, appropriate for CH3CCH. As can be seen in the right panels, all of
the observed abundances can be reached after about $10^{5}$ years of post
shock evolution. For CH3CCH, 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 CH3CCH if the density is around $10^{4}$
cm-3.

Figure 15: Results from the NAUTILUS two 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 cold cloud
evolution stage and the right panels show the warm-up stage for three
different densities. Figure 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.

#### H2CO and $c$-C3H2

To model H2CO and $c$-C3H2 we first assumed that the physical conditions were
identical to those for CH3CCH and CH3OH; namely that the cloud experienced a
10 K cold stage followed by a 30 K warm-up stage. We also used the H2 column
density map calculated at 30 K to determine the H2CO and $c$-C3H2 abundances.
Under these conditions, the H2CO abundance ranges from $3\times 10^{-9}$ in
the cloud core to $10^{-8}$ in the north western cloud edges. The $c$-C3H2
abundances range from $4\times 10^{-11}$ in the cloud core to $2\times
10^{-10}$ along the western and northern edges. Under these conditions, the
observed $c$-C3H2 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$-C3H2 abundances using a temperature of 20 K to calculate the H2 column
density. H2CO, 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$-C3H2 (see
Fig. 9), it is not for H2CO where the average temperature is about 40-50 K.
For H2CO we used temperatures of 50 K to calculate the H2 column density and
then the abundance map. Given the higher temperature, the resultant H2 column
densities are lower, resulting in slightly increased H2CO abundances which now
vary from $5\times 10^{-9}$ to $2\times 10^{-8}$. 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 produced the observed
abundances.

####  5.3.2 IRAS 20126 chemical model results

#### CH3CCH and CH3OH

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 CH3CCH and
CH3OH LTE model is dominated by one main component. To produce the
corresponding H2 column density map, we used archival JCMT SCUBA-1 850 $\mu 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 CH3CCH
and CH3OH LTE model column density maps by the H2 column density map and found
abundances varying from 10-10 and 10-9 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^{5}$ years under cold
conditions (n$=10^{4}$ cm-3, T = 10 K) starting with general ISM abundances,
and then warmed the gas to 30 K to evolve for another $10^{6}$ years.

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, AV = 50 (i.e. N(H2) $=10^{23}$ cm-2), and, from the top, $n=10^{4}$
cm-3, $n=5\times 10^{4}$ cm-3, and $n=10^{5}$ cm-3. The dotted black
horizontal lines indicate the observed abundance limits of $10^{-10}$ and
$10^{-9}$. Figure 15 shows that although the observed CH3CCH abundances can be
reproduced after a few $\times 10^{5}$ years in the the warm up phase, the
model CH3OH 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 CH3OH 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\times 10^{4}$ cm-3 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^{-10}$ and $10^{-9}$.
Throughout the duration of the shock, and for the first $10^{5}$ years of the
post shock evolution, the CH3OH abundance remains above the observed values.
After this, CH3OH begins to refreeze onto the dust grain surfaces and, in
doing so, passes through the range of observed abundances. While CH3CCH 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^{6}$
years.

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 CH3OH abundances over the entirety of IRAS 20126. They
are not, however, needed to produce the observed abundances of CH3CCH. This is
as seen in AFGL 2591, with the exception of the northern region. Our results
for CH3CCH in IRAS 20126 are uniformly consistent with Warm Carbon Chain
Chemistry.

#### H2CO and $c$-C3H2

The modelling of H2CO and $c$-C3H2 proceeded in a fashion similar to that done
for AFGL 2591. We first assumed the same conditions as for CH3OH and CH3CCH,
and then used the temperatures determined from our LTE model (see Fig. 10).
With the H2 column density map calculated at 30 K we derived H2CO abundance
ranges from $6\times 10^{-11}$ around the cloud edges to $3\times 10^{-10}$ in
the cloud core. The $c$-C3H2 abundance ranges from $3\times 10^{-11}$ in the
cloud core and eastern edges to $10^{-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$-C3H2 is about 8 K and
fairly uniform across the map. For H2CO, the temperature is about 80 K and
also fairly uniform. We used these two temperatures to calculate the
respective H2 column density maps and found that, for T = 8 K, the $c$-C3H2
abundance varies from $4\times 10^{-12}$ in the cloud core and eastern edges
to $10^{-11}$ along the western edges. For T = 80 K, the H2CO abundance ranges
from $2\times 10^{-10}$ around the cloud edges to $10^{-9}$ in the cloud core.
The chemical models are able to produce these observed abundances as well. For
$c$-C3H2 at 8 K, the models match the observations early in the cold cloud
stage. For H2CO 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. 1.

detections of numerous rotational transition lines of two carbon chain
molecules, CH3CCH and $c$-C3H2, and two molecules related to the complex
organic family, CH3OH and H2CO, in both sources.

  2. 2.

a 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. 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. 4.

chemical evolution modelling with NAUTILUS that demonstrates that the observed
abundances of CH3CCH 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. CH3OH, 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. PF 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. PF would also like to thank Pablo Torne
and Monica Rodriguez at IRAM for their help observing with the IRAM telescope.
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 Region 3.

## References

  * Aikawa et al. (2020) Aikawa, Y., Furuya, K., Yamamoto, S., & Sakai, N. 2020, ApJ, 897, 110
  * Aikawa et al. (2008) Aikawa, Y., Wakelam, V., Garrod, R. T., & Herbst, E. 2008, ApJ, 674, 984
  * Andron et al. (2018) Andron, I., Gratier, P., Majumdar, L., et al. 2018, MNRAS, 481, 5651
  * Askne et al. (1984) Askne, J., Hoglund, B., Hjalmarson, A., & Irvine, W. M. 1984, A&A, 130, 311
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167
  * Avery et al. (1976) Avery, L. W., Broten, N. W., MacLeod, J. M., Oka, T., & Kroto, H. W. 1976, ApJL, 205, L173
  * Ball et al. (1970) Ball, J. A., Gottlieb, C. A., Lilley, A. E., & Radford, H. E. 1970, ApJ, 162, L203
  * Batrla & Menten (1988) Batrla, W. & Menten, K. M. 1988, ApJ, 329, L117
  * Bergin et al. (1994) Bergin, E. A., Goldsmith, P. F., Snell, R. L., & Ungerechts, H. 1994, ApJ, 431, 674
  * Beuther et al. (2002) Beuther, H., Schilke, P., Menten, K. M., et al. 2002, ApJ, 566, 945
  * Blake et al. (1987) Blake, G. A., Sutton, E. C., Masson, C. R., & Phillips, T. G. 1987, ApJ, 315, 621
  * Bohlin et al. (1978) Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132
  * Bottinelli et al. (2004) Bottinelli, S., Ceccarelli, C., Neri, R., et al. 2004, ApJ, 617, L69
  * Bottinelli et al. (2007) Bottinelli, S., Ceccarelli, C., Williams, J. P., & Lefloch, B. 2007, A&A, 463, 601
  * Bouvier et al. (2020) Bouvier, M., López-Sepulcre, A., Ceccarelli, C., et al. 2020, A&A, 636, A19
  * Broten et al. (1978) Broten, N. W., Oka, T., Avery, L. W., MacLeod, J. M., & Kroto, H. W. 1978, ApJL, 223, L105
  * Calcutt et al. (2019) Calcutt, H., Willis, E. R., Jørgensen, J. K., et al. 2019, A&A, 631, A137
  * Caselli et al. (1993) Caselli, P., Hasegawa, T. I., & Herbst, E. 1993, ApJ, 408, 548
  * Cazaux et al. (2003) Cazaux, S., Tielens, A. G. G. M., Ceccarelli, C., et al. 2003, ApJ, 593, L51
  * Ceccarelli (2004) Ceccarelli, C. 2004, in Astronomical Society of the Pacific Conference Series, Vol. 323, Star Formation in the Interstellar Medium: In Honor of David Hollenbach, ed. D. Johnstone, F. C. Adams, D. N. C. Lin, D. A. Neufeeld, & E. C. Ostriker, 195
  * Cesaroni et al. (1999) Cesaroni, R., Felli, M., Jenness, T., et al. 1999, A&A, 345, 949
  * Cesaroni et al. (1997) Cesaroni, R., Felli, M., Testi, L., Walmsley, C. M., & Olmi, L. 1997, A&A, 325, 725
  * Cesaroni et al. (2005) Cesaroni, R., Neri, R., Olmi, L., et al. 2005, A&A, 434, 1039
  * Charnley et al. (1997) Charnley, S. B., Tielens, A. G. G. M., & Rodgers, S. D. 1997, The Astrophysical Journal, 482, L203
  * Chen et al. (2016) Chen, H.-R. V., Keto, E., Zhang, Q., et al. 2016, ApJ, 823, 125
  * Favre et al. (2020) Favre, C., Vastel, C., Jimenez-Serra, I., et al. 2020, A&A, 635, A189
  * Fayolle et al. (2015) Fayolle, E. C., Öberg, K. I., Garrod, R. T., van Dishoeck, E. F., & Bisschop, S. E. 2015, A&A, 576, A45
  * Garrod & Widicus Weaver (2013) Garrod, R. T. & Widicus Weaver, S. L. 2013, Chemical Reviews, 113, 8939
  * Giannetti et al. (2017) Giannetti, A., Leurini, S., Wyrowski, F., et al. 2017, A&A, 603, A33
  * Gieser et al. (2019) Gieser, C., Semenov, D., Beuther, H., et al. 2019, A&A, 631, A142
  * Graninger et al. (2016) Graninger, D. M., Wilkins, O. H., & Öberg, K. I. 2016, ApJ, 819, 140
  * Harvey-Smith et al. (2008) Harvey-Smith, L., Soria-Ruiz, R., Duarte-Cabral, A., & Cohen, R. J. 2008, MNRAS, 384, 719
  * Hasegawa & Mitchell (1995) Hasegawa, T. I. & Mitchell, G. F. 1995, ApJ, 451, 225
  * Helmich & van Dishoeck (1997) Helmich, F. P. & van Dishoeck, E. F. 1997, A&AS, 124, 205
  * Herbst & van Dishoeck (2009) Herbst, E. & van Dishoeck, E. F. 2009, ARA&A, 47, 427
  * Hickson et al. (2016) Hickson, K. M., Wakelam, V., & Loison, J.-C. 2016, Molecular Astrophysics, 3-4, 1
  * Higuchi et al. (2018) Higuchi, A. E., Sakai, N., Watanabe, Y., et al. 2018, ApJS, 236, 52
  * Imai et al. (2022) Imai, M., Oya, Y., Svoboda, B., et al. 2022, The Astrophysical Journal, 934, 70
  * Imai et al. (2016) Imai, M., Sakai, N., Oya, Y., et al. 2016, ApJL, 830, L37
  * Irvine et al. (1981) Irvine, W. M., Hoglund, B., Friberg, P., Askne, J., & Ellder, J. 1981, ApJ, 248, L113
  * Isokoski et al. (2013) Isokoski, K., Bottinelli, S., & van Dishoeck, E. F. 2013, A&A, 554, A100
  * Jiménez-Serra et al. (2008) Jiménez-Serra, I., Caselli, P., Martín-Pintado, J., & Hartquist, T. W. 2008, A&A, 482, 549
  * Jiménez-Serra et al. (2012) Jiménez-Serra, I., Zhang, Q., Viti, S., Martín-Pintado, J., & de Wit, W. J. 2012, ApJ, 753, 34
  * Johnston et al. (2013) Johnston, K. G., Shepherd, D. S., Robitaille, T. P., & Wood, K. 2013, A&A, 551, A43
  * Jones et al. (1994) Jones, A. P., Tielens, A. G. G. M., Hollenbach, D. J., & McKee, C. F. 1994, ApJ, 433, 797
  * Jørgensen et al. (2020) Jørgensen, J. K., Belloche, A., & Garrod, R. T. 2020, ARA&A, 58, 727
  * Kalenskii et al. (2002) Kalenskii, S. V., Slysh, V. I., Val’tts, I. E., Winnberg, A., & Johansson, L. E. B. 2002, in Cosmic Masers: From Proto-Stars to Black Holes, ed. V. Migenes & M. J. Reid, Vol. 206, 195
  * Kaźmierczak-Barthel et al. (2014) Kaźmierczak-Barthel, M., van der Tak, F. F. S., Helmich, F. P., et al. 2014, A&A, 567, A53
  * Kroto et al. (1978) Kroto, H. W., Kirby, C., Walton, D. R. M., et al. 1978, ApJL, 219, L133
  * Kuiper et al. (1984) Kuiper, T. B. H., Kuiper, E. N. R., Dickinson, D. F., Turner, B. E., & Zuckerman, B. 1984, ApJ, 276, 211
  * Lada et al. (1984) Lada, C. J., Thronson, H. A., J., Smith, H. A., Schwartz, P. R., & Glaccum, W. 1984, ApJ, 286, 302
  * Lindberg et al. (2016) Lindberg, J. E., Charnley, S. B., & Cordiner, M. A. 2016, The Astrophysical Journal Letters, 833, L14
  * Little et al. (1978) Little, L. T., MacDonald, G. H., Riley, P. W., & Matheson, D. N. 1978, MNRAS, 183, 45P
  * Lovas et al. (1976) Lovas, F. J., Johnson, D. R., Buhl, D., & Snyder, L. E. 1976, ApJ, 209, 770
  * Macdonald, G. H. et al. (1996) Macdonald, G. H., Gibb, A. G., Habing, R. J., & Millar, T. J. 1996, A&AS, 119, 333
  * Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425
  * Mookerjea et al. (2012) Mookerjea, B., Hassel, G. E., Gerin, M., et al. 2012, A&A, 546, A75
  * Moscadelli et al. (2011) Moscadelli, L., Cesaroni, R., Rioja, M. J., Dodson, R., & Reid, M. J. 2011, A&A, 526, A66
  * Newville et al. (2014) Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2014, LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python
  * Okoda et al. (2023) Okoda, Y., Oya, Y., Francis, L., et al. 2023, FAUST VII. Detection of A Hot Corino in the Prototypical Warm Carbon-Chain Chemistry Source IRAS 15398-3359
  * Ossenkopf & Henning (1994) Ossenkopf, V. & Henning, T. 1994, A&A, 291, 943
  * Oya et al. (2017) Oya, Y., Sakai, N., Watanabe, Y., et al. 2017, ApJ, 837, 174
  * Palau et al. (2017) Palau, A., Walsh, C., Sánchez-Monge, Á., et al. 2017, MNRAS, 467, 2723
  * Park et al. (2006) Park, I. H., Wakelam, V., & Herbst, E. 2006, A&A, 449, 631
  * Plambeck & Menten (1990) Plambeck, R. L. & Menten, K. M. 1990, ApJ, 364, 555
  * Rachford et al. (2009) Rachford, B. L., Snow, T. P., Destree, J. D., et al. 2009, ApJS, 180, 125
  * Rodgers & Charnley (2001) Rodgers, S. D. & Charnley, S. B. 2001, ApJ, 546, 324
  * Rodríguez-Garza et al. (2017) Rodríguez-Garza, C. B., Kurtz, S. E., Gómez-Ruiz, A. I., et al. 2017, The Astrophysical Journal Supplement Series, 233, 4
  * Ruaud et al. (2016) Ruaud, M., Wakelam, V., & Hersant, F. 2016, MNRAS, 459, 3756
  * Rygl et al. (2012) Rygl, K. L. J., Brunthaler, A., Sanna, A., et al. 2012, A&A, 539, A79
  * Sakai et al. (2009) Sakai, N., Sakai, T., Hirota, T., Burton, M., & Yamamoto, S. 2009, ApJ, 697, 769
  * Sakai et al. (2008) Sakai, N., Sakai, T., Hirota, T., & Yamamoto, S. 2008, ApJ, 672, 371
  * Sakai et al. (2010) Sakai, N., Shiino, T., Hirota, T., Sakai, T., & Yamamoto, S. 2010, ApJL, 718, L49
  * Sakai & Yamamoto (2013) Sakai, N. & Yamamoto, S. 2013, Chemical Reviews, 113, 8981
  * Sanna et al. (2012) Sanna, A., Reid, M. J., Carrasco-González, C., et al. 2012, ApJ, 745, 191
  * Santos et al. (2022) Santos, J. C., Bronfman, L., Mendoza, E., et al. 2022, ApJ, 925, 3
  * Saul et al. (2015) Saul, M., Tothill, N. F. H., & Purcell, C. R. 2015, ApJ, 798, 36
  * Schiff & Bohme (1979) Schiff, H. I. & Bohme, D. K. 1979, ApJ, 232, 740
  * Shepherd et al. (2000) Shepherd, D. S., Yu, K. C., Bally, J., & Testi, L. 2000, ApJ, 535, 833
  * Snyder & Buhl (1973) Snyder, L. E. & Buhl, D. 1973, Nature Physical Science, 243, 45
  * Snyder et al. (1969) Snyder, L. E., Buhl, D., Zuckerman, B., & Palmer, P. 1969, Phys. Rev. Lett., 22, 679
  * Soma et al. (2015) Soma, T., Sakai, N., Watanabe, Y., & Yamamoto, S. 2015, The Astrophysical Journal, 802, 74
  * Spezzano et al. (2016) Spezzano, S., Bizzocchi, L., Caselli, P., Harju, J., & Brünken, S. 2016, A&A, 592, L11
  * Sridharan et al. (2005) Sridharan, T. K., Williams, S. J., & Fuller, G. A. 2005, ApJ, 631, L73
  * Su et al. (2007) Su, Y.-N., Liu, S.-Y., Chen, H.-R., Zhang, Q., & Cesaroni, R. 2007, The Astrophysical Journal, 671, 571
  * Suri et al. (2021) Suri, S., Beuther, H., Gieser, C., et al. 2021, A&A, 655, A84
  * Taniguchi et al. (2019) Taniguchi, K., Herbst, E., Caselli, P., et al. 2019, ApJ, 881, 57
  * Taniguchi et al. (2021) Taniguchi, K., Majumdar, L., Takakuwa, S., et al. 2021, ApJ, 910, 141
  * Taniguchi et al. (2018a) Taniguchi, K., Miyamoto, Y., Saito, M., et al. 2018a, ApJ, 866, 32
  * Taniguchi et al. (2018b) Taniguchi, K., Saito, M., Majumdar, L., et al. 2018b, ApJ, 866, 150
  * Thaddeus et al. (1985) Thaddeus, P., Vrtilek, J. M., & Gottlieb, C. A. 1985, ApJ, 299, L63
  * Trinidad et al. (2003) Trinidad, M. A., Curiel, S., Cantó, J., et al. 2003, ApJ, 589, 386
  * Turner et al. (1999) Turner, B. E., Terzieva, R., & Herbst, E. 1999, ApJ, 518, 699
  * van der Tak et al. (2007) van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627
  * van der Tak et al. (1999) van der Tak, F. F. S., van Dishoeck, E. F., Evans, Neal J., I., Bakker, E. J., & Blake, G. A. 1999, ApJ, 522, 991
  * van der Tak et al. (2000) van der Tak, F. F. S., van Dishoeck, E. F., Evans, Neal J., I., & Blake, G. A. 2000, ApJ, 537, 283
  * van der Tak et al. (2006) van der Tak, F. F. S., Walmsley, C. M., Herpin, F., & Ceccarelli, C. 2006, A&A, 447, 1011
  * van der Wiel et al. (2013) van der Wiel, M. H. D., Pagani, L., van der Tak, F. F. S., Kaźmierczak, M., & Ceccarelli, C. 2013, A&A, 553, A11
  * van der Wiel et al. (2011) van der Wiel, M. H. D., van der Tak, F. F. S., Spaans, M., et al. 2011, A&A, 532, A88
  * Vastel et al. (2014) Vastel, C., Ceccarelli, C., Lefloch, B., & Bachiller, R. 2014, ApJ, 795, L2
  * Vrtilek et al. (1987) Vrtilek, J. M., Gottlieb, C. A., & Thaddeus, P. 1987, ApJ, 314, 716
  * Wakelam et al. (2017) Wakelam, V., Bron, E., Cazaux, S., et al. 2017, Molecular Astrophysics, 9, 1
  * Wang et al. (2012) Wang, K. S., van der Tak, F. F. S., & Hogerheijde, M. R. 2012, A&A, 543, A22
  * Whittet (1981) Whittet, D. C. B. 1981, MNRAS, 196, 469
  * Wilking et al. (1990) Wilking, B. A., Blackwell, J. H., & Mundy, L. G. 1990, AJ, 100, 758
  * Wu et al. (2019) Wu, Y., Lin, L., Liu, X., et al. 2019, A&A, 627, A162
  * Xu et al. (2011) Xu, J.-L., Wang, J.-J., & Ning, C.-C. 2011, The Astrophysical Journal, 744, 175
  * Zhang et al. (2001) Zhang, Q., Hunter, T. M., Sridharan, T. K., & Kawamura, J. H. 2001, in Astronomical Society of the Pacific Conference Series, Vol. 235, Science with the Atacama Large Millimeter Array, ed. A. Wootten, 191–194
  * Zhang et al. (2005) Zhang, Q., Hunter, T. R., Brand, J., et al. 2005, ApJ, 625, 864
  * Zhang et al. (1998) Zhang, Q., Hunter, T. R., & Sridharan, T. K. 1998, The Astrophysical Journal, 505, L151
  * Zhu et al. (2017) Zhu, H., Tian, W., Li, A., & Zhang, M. 2017, MNRAS, 471, 3494

##  Appendix A Additional LTE model description

Table 6: Parameter inputs for the LTE model.

Species | Ntot | Tex | vlsr | FWHM
---|---|---|---|---
| (cm-2) | (K) | (km s-1) | (km s-1)
AFGL 2591 |  |  |  |
CH3CCH | 1 $\times$ 1011, 1 $\times$ 1017, 1 $\times$ 1014 | 3, 100, 40 | -10, -2, -5.7 | 0.5, 4, 2
$c$-C3H2 | 1 $\times$ 1011, 1 $\times$ 1017, 1 $\times$ 1014 | 3, 100, 40 | -10, -2, -5.7 | CH3CCH FWHM
H2CO | 1 $\times$ 1011, 1 $\times$ 1017, 1 $\times$ 1014 | 3, 150, 90 | -10, -2, -6 | 1, 7, 4
CH3OH C1 | 1 $\times$ 1011, 1 $\times$ 1017, 1 $\times$ 1014 | 18 | -10, -1, -6 | 1, 7, 3
CH3OH C2 | 1 $\times$ 1011, 1 $\times$ 1017, 1 $\times$ 1014 | 8 | -10, -2, -6 | 3.9
CH3OH C3 | 1 $\times$ 1011, 1 $\times$ 1017, 1 $\times$ 1014 | 8 | -9.4 | 2.5
IRAS 20126 |  |  |  |
CH3CCH | 1 $\times$ 1011, 1 $\times$ 1017, 1 $\times$ 1014 | 3, 100, 30 | -5, -1, -3.1 | 0.5, 4, 2
$c$-C3H2 C1 | 1 $\times$ 1011, 1 $\times$ 1017, 1 $\times$ 1014 | 3, 100, 10 | -5, -1, -3.1 | 0.5, 4, 2.4
$c$-C3H2 C2 | 1.0 $\times$ 1013 | 3, 100, 10 | C1 vlsr | 0.5
H2CO | 1 $\times$ 1011, 1 $\times$ 1017, 1 $\times$ 1014 | 3, 150, 50 | -6, -1, -3.6 | 1, 8, 6.9
CH3OH C1 | 1 $\times$ 1011, 1 $\times$ 1017, 1 $\times$ 1014 | 3, 150, 50 | -6, -1, -3.6 | 1, 8, 6.9
CH3OH C2 | 2.1 $\times$ 1013 | 10 | C1 vlsr-0.2 | 1.2
111111Where the parameter is variable, the inputs are minimum, maximum, and
starting value, and where the parameter is fixed, the input is the starting
value. If a parameter is linked to another, the reference parameter is noted
in the table. The final parameter, source size, is always fixed at
20$\arcsec$.

(a)

(b)

Figure 17: Spectral line fits from the LTE model in AFGL 2591. All transition
lines are shown 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. The model spectrum is in colour while the
observed data is in black. 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.

(a)

(b)

Figure 18: Spectral line fits from the LTE model in AFGL 2591. All transition
lines are shown 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. The model
spectrum is in colour while the observed data is in black. 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. Figure 19: Parameters maps 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

(a)

(b)

Figure 20: Same as Fig. 17 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.

(a)

(b)

Figure 21: Same as Fig. 18 for (a) CH3OH and (b) H2CO in IRAS 20126. Two
components are used to fit CH3OH while one component is used to fit H2CO.
Figure 22: 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 $\times$ 1013 cm-2 and the FWHM at
0.5 km s-1. For CH3OH, Ntot is fixed at 2.1 $\times$ 1013 cm-2, Tex at 10 K
the FWHM at 1.2 km s-1.