© L. Bouscasse et al. 2022

1 Introduction

The molecular richness of the star-forming gas was first recognised towards sites of high-mass star formation. Hot molecular cores have been identified as birthplaces of high-mass stars and clusters where high densities and temperatures prevail. Temperatures above ~100K lead to the sublimation of icy dust mantles, which together with the high temperatures give rise to an emergence of molecular complexity (e.g. Kurtz et al. 2000; Belloche et al. 2016; Jørgensen et al. 2020). The molecular complexity manifests in the detection of O-bearing, N-bearing, and S-bearing deuterated molecules, carbon chains, and complex organic molecules¹ (COMs; for a review, see e.g. Herbst & van Dishoeck 2009). Even though hot cores show similarities in their molecular content, they have been found to exhibit a large diversity in terms of relative molecular abundances (Walmsley & Schilke 1993; Bisschop et al. 2007; Belloche et al. 2013; Widicus Weaver et al. 2017; Minh et al. 2018). To reveal the origin of this diversity, the next step is to investigate the molecular composition of recently recognised (high-mass) protostellar envelopes in a young evolutionary stage in great detail, where the early warm-up phase chemistry could reveal the first steps towards a more complex chemical evolution.

Unbiased spectral surveys have been proven to be extremely useful to pinpoint the physical and chemical conditions around (high-mass) protostars. Simultaneously observable large band-widths grant access to several transitions of a given molecule, allowing us to cover a broad range of upper-level energies to estimate the physical conditions giving rise to the molecular emission. This could potentially reveal the origin of the emitting gas even if the region is unresolved in the beam of single-dish telescopes. The large bandwidth of modern receivers meant that spectral surveys became relatively easy to perform, with a growing number of such programs tracing a range of star formation environments from the quiescent to the extreme regions: for example, the TIMASS survey (The IRAS 16293-2422 Millimeter And Submillimeter Spectral Survey, Caux et al. 2011) and ASAI (Astrochemical Surveys At IRAM, Lefloch et al. 2018), which targeted mainly low-mass protostars, Belloche et al. (2013), who targeted the extreme star formation region in Sgr B2, and several others listed in Herbst & van Dishoeck (2009). There is a growing number of such studies at high angular resolution as well, such as the PILS survey (the ALMA Proto-stellar Interferometric Line Survey, Jørgensen et al. 2016, SOLIS (Seeds Of Life In Space, Ceccarelli et al. 2017), and the EMoCA (Exploring Molecular Complexity with ALMA) and ReMoCA (Re-exploring molecular complexity with ALMA) surveys of Sgr B2(N) (Belloche et al. 2016, 2019).

The high-mass protostar G328.2551-0.5321 was first identified in the SPARKS (Search for High-mass Protostars with ALMA revealed up to kilo-parsec scales) project (Csengeri et al. 2017b; Csengeri et al. 2018) based on the ATLASGAL (APEX Telescope Large Area Survey of the Galaxy) survey (Schuller et al. 2009; Csengeri et al. 2014). It is located at a distance of ${2.5}_{-0.5}^{+1.7}$ kpc with a systemic velocity of −43.1 kms⁻¹ (Csengeri et al. 2018). Csengeri et al. (2018) found observational signatures of accretion shocks associated with a putative accretion disk witnessed for the first time towards a high-mass protostar. The physical properties such as the temperature and the molecular characteristics of these accretion shocks were different compared to the rest of the envelope (Csengeri et al. 2019). Towards the inner envelope at Τ ~ 100 K, about ten COMs were identified using the 7.5 GHz limited bandwidth that is available in a single setup with ALMA. The protostar is expected to be in a stage before the onset of strong ionising radiation (see more details in Csengeri et al. 2018), and the radiatively heated inner region is still compact (≤1000au). Because the molecular emission is dominated by the accretion shocks and the cold gas, Csengeri et al. (2019) suggested that this object might be a precursor to a classical, radiatively heated hot core. This object has been found to be isolated from 0.3 pc down to 400au scales with no other compact and bright sub-millimetre source in its vicinity. It is therefore an appropriate target for single-dish observations and is an excellent laboratory for studying the onset of radiative heating in the chemistry of the early warm-up phase.

Here we study the full molecular composition of G328.2551-0.5321 using the Atacama Pathfinder EXperiment telescope (APEX, Gusten et al. 2006). We performed an unbiased spectral survey in the 2 mm, 1.2 mm, and 0.8 mm atmospheric bands in order to study the physical conditions and the molecular composition of the dense gas in this hot core precursor. The paper is organised as follows: in Sect. 2 we present the observations and the data reduction. In Sect. 3, we show the results and the molecular composition of the protostellar envelope. In Sect. 4, we describe the envelope structure. In Sect. 5, we discuss the composition of light molecules, and in Sect. 6 we discuss the detection and properties of COMs in the envelope. In Sect. 7, we compare this emerging hot core with hot cores from the literature. Finally, in Sect. 8 we present our conclusions.

Fig. 1 Continuous frequency coverage of the spectral survey. The black curve shows the atmospheric transmission as computed by the ATM model written by Juan Pardo (Pardo et al. 2002) with a water vapour of 1.2 mm.

2 Observations and data reduction

We performed an unbiased spectral survey from 159 GHz to 374 GHz towards the protostellar object associated with the ATLASGAL clump G328.2551-0.5321 (RA 15h57m59.791s and Dec. –53°58′00.56″ (J2000); Csengeri et al. 2018)² with the APEX telescope on the Chajnantor plateau, in the Chilean Andes.

The observations cover the 0.8 mm band using the FLASH⁺ (First Light APEX Submillimeter Heterodyne) receiver (Klein et al. 2014), the 1.2 mm band with the PI230 receiver, and finally the 2 mm band with the SEPIA180 (Swedish-ESO PI Instrument for APEX; Belitsky et al. 2018) and PI230 receivers³ (Fig. 1). In addition, we also observed additional setups as a bonus around 408 GHz and 465 GHz with the FLASH460 channel in parallel of the FLASH345 observations. The observations were made in position-switching mode with a reference position offset of 600″, −600″ in right ascension and declination with respect to the position of G328.2551-0.5321. The reference position was chosen to be a nearby position to ensure spectral baseline stability. The telescope focus was checked with Mars and Jupiter. The pointing accuracy is ~2″, and was checked with regular pointings on IRAS 15194-5115, NGC 6072 and RAFGL2135. All three receivers are sideband separated with a sideband rejection of 15 dB for PI230 and FLASH345 (Güsten priv. comm., Klein et al. 2014) and 20 dB for SEPIA180 (Belitsky et al. 2018). The PI230 receiver was used covering 4–12GHz of the intermediate frequency (IF) band between 200 GHz and 270 GHz, and 4–8GHz between 186 GHz and 210 GHz. SEPIA180 and FLASH345 were both used in their original configuration, providing an IF band of 4–8 GHz. In order to identify line contamination from the other sideband (ghost lines) and occasionally appearing spikes, we covered the frequency range twice with tunings shifted by 4 GHz for the PI230 and 2 GHz for FLASH345, SEPIA180, and PI230 in the 4–8GHz configuration, and hence used 12, 24, and 12 setups for the 2 mm, 1.2 mm, and 0.8 mm atmospheric windows, respectively (Fig. 2). The half-power beam width varies from 39″ at 159 GHz to 16” at 374 GHz, corresponding to 0.5 pc and 0.2 pc physical scales at the distance of the source, respectively.

The median system temperature is 179 K over the complete spectral survey, and the noise varies between 11 and 25 mK at 0.7 km s⁻¹ on the main-beam temperature scale (T_mb) (Table 1), the variation depending on the receiver temperature and the observing conditions. In total we cover a frequency range of 206 GHz with this spectral survey, which was observed for a total of 7 h (including the reference position).

The data were calibrated with the APEX online calibrator. We performed the data reduction with the Aug. 18 and March 19b version of GILDAS⁴ softwares. For the data reduction, we first averaged each subscan of every backend. Then, we compared the intensities of the lines in each polarisation (for SEPIA180 and PI230), which were then averaged together. For each backend of every setup, we removed a zeroth-order baseline selecting emission free channels. Finally, we compared the two frequency coverages with different tunings and identified about ten ghost lines, which were blanked in the data, while channels with spikes were masked. The FLASH345 data from April 2019 were shifted by 0.25 km s⁻¹ to fix a problem related to the receiver. We corrected this shift in the data and then averaged them with the data from September 2019.

We converted the data from antenna temperature $(T_{\text{A}}*)$ to main-beam temperature (T_mb) using the telescope efficiencies determined from planet observations that were regularly observed during the project. These efficiencies were then applied to the whole survey using the Ruze formula: ${\eta }_{\text{mb}}={\eta }_0\times {\text{e}}^{-{\left(\frac{4\pi \delta \nu }{c}\right)}^2}$, where η_mb is the main-beam efficiency, η₀ depends on the receiver, δ is the surface accuracy, ν is the frequency, and c is the speed of light. The main-beam efficiencies for some frequencies are given in Table 1. A forward efficiency of 95% was assumed for all receivers. Finally, we smoothed the native velocity resolution to a common resolution of 0.7 km s⁻¹. The average rms noise level (σ) over the entire band is 18 mK. The parameters of the observations are detailed in Table 1.

Fig. 2 Observed spectral windows for the 1 mm band with the PI230 receiver at APEX. The solid lines represent the inner backend units, and the dashed lines show the outer backend units.

3 Results

We show the spectrum towards G328.2551-0.5321 from 252 GHz to 258 GHz in Fig. 3. The complete spectrum is shown in Appendix C. In total, 3055 lines were detected over the 206 GHz covered by the survey above the 3σ level. Assuming Gaussian statistics, this gives us about ninelines as false-positive detections when we consider a 3σ threshold for the peak line temperature. The overall line density is found to be ~ 15 lines/GHz above 3σ. This line density is significantly lower than the 3 mm band spectrum of the bright hot core, Sgr B2(N), which shows a line density of 102 lines/GHz⁵ above 4σ (Belloche et al. 2013) with a noise of 22 mK $(T_{\text{A}}*)$ at a spectral resolution of 1.2 km s⁻¹, corresponding to a noise of 28 mK at a spectral resolution of 0.7 km s⁻¹.

3.1 Line identification

We identified the molecules at the origin of the observed spectral lines using the CDMS⁶ (Müller et al. 2005) and JPL⁷ (Pickett et al. 1998) catalogues with the help of the Weeds package for CLASS (Maret et al. 2011). These two databases provide the rest frequencies, the upper-level energies, E_up, and the Einstein coefficients, A_ij, for many molecules. To identify the emitting molecules and their transitions, we first visually identified the brightest lines and looked up the most probable molecular line with the highest A_ij value at the corresponding frequency in the spectroscopic databases. Then we searched for all the transitions from this molecule that we expect to detect considering their upper-level energies and Einstein coefficients, or synthetic local thermal equilibrium (LTE) models (Sect. 3.4) in the band above a minimum signal-to-noise ratio of 3 (Table 3). If all the expected transitions were detected, we considered this molecule as identified. We repeated this process with the next brightest line. A final iteration for the robustness of our detections was made against the LTE synthetic spectrum (Sect. 3.4), where we verified that all expected transitions were detected and no line was present in the synthetic spectrum that was not consistent with the data.

This procedure was efficient for molecules with a few atoms having up to 60 transitions in the band. We found a few molecules with only one transition above our detection threshold of a minimum signal-to-noise ratio of 3. While the transitions of HCNH⁺ (3–2), DNC (3–2), NS⁺ (4–3), and PN (5–4) are clearly detected above the 3σ threshold (Fig. 4), we claim these molecules as tentatively detected because they only have one transition.

The heavier, i.e. higher-mass, COMs have a considerably larger number of transitions in the band due to their large partition function. Therefore, we used the LTE modelling (see Sect. 3.4) of their spectra and visually checked whether all the expected lines were detected with intensities consistent with the synthetic spectrum. Altogether, we were able to firmly identify 2935 lines from all species. About 120 weak lines above 3σ remained unidentified.

In total, we securely detect 39 species plus 22 isotopo-logues (Tables 2 and 3), and we tentatively detect four additional molecules. In Table 2 we group the molecules according to their atomic content, such as carbon chains, S-bearing, O-bearing, and N-bearing deuterated species, and COMs. The groups with the largest number of molecules are the COMs and S-bearing molecules. About 1600 of the 2935 identified lines are from O-bearing species (including O-bearing COMs), ~300 lines from S-bearing molecules, ~70 lines from carbon chains, and ~900 lines from N-bearing molecules (including N-bearing COMs). Rotational transitions of HC₃N from within its first vibrationally excited state (υ₇ = 1) and CH₃OH from within its torsionally excited state (υ_t = 1) are detected.

Fig. 3 Spectrum towards G328.2551-0.5321 between 252 GHz and 258 GHz. The filled grey histogram shows the spectra, and the coloured labels show the detected transitions of all the light molecules and some transitions of the COMs. The green line represents the LTE model obtained with Weeds. The upper panel for each frequency range shows the labels for the small molecules, and the lower panel corresponds to a zoom in temperature for better visibility of the COMs which are label.

Fig. 4 Spectra, in black, of DNC(3–2), NS+ (4–3), PN(5–4), HCNH+ (3–2) shown as filled grey histograms. The dashed red line shows the systemic velocity of the source at –43.13 kms−1 (Csengeri et al. 2018). The blue line represents the best Gaussian fit for each line. Eup represents the upper-level energy of the transition in temperature units.

3.2 Line profiles and fitting

Our spectral resolution of 0.7 km s⁻¹ is sufficient to resolve the lines well, and we note a variety of spectral shapes associated with different molecules; some examples are shown in Fig. 5. Most of the molecules with optically thin emission, for example HC¹⁸O⁺, SO₂, and ²⁹SiO, show a single Gaussian component. However, a fraction of the lines, such as CH₃OH, HNCO, HC₃N, and OCS, requires a two-component Gaussian fit using a narrow and a broad component. Methanol lines with a high upper-level energy (E_up ≥ 150 K) exhibit a profile consisting of two velocity components that are offset from the source v_lsr velocity by about ±4kms⁻¹, similarly as reported in (Csengeri et al. 2018) towards the accretions shocks on the basis of interferometric data. Finally, a few lines show non-Gaussian profiles that are heavily affected by self-absorption, such as the lines of HCO⁺, HCN, HNC, and H₂S. Broad line-wings associated with high-velocity gas are observed in the CO transitions, for example. The outflow emission was associated by eye because the outflow appears as non-Gaussian residual in the fitting process. Table 3 indicates the molecules in which some transitions exhibit outflow wings. We distinguish the cold and warm components in a quantitative way (see Sect. 4).

For the lines showing Gaussian profiles, we used the fitting procedure of CLASS and derived the line peak temperature, T_mb, velocity at line peak, υ, full width at half maximum, Δυ, and integrated intensity, W, of each unblended line. Molecules with a resolved hyperfine structure (HFS), such as NS, NO, and CCH (see Fig. 6 for an example and Table 4 for the list of the fitted molecules with an HFS procedure) were fitted based on the relative intensities of the HFS components from the JPL or CDMS catalogues using the HFS procedure of CLASS. An exception is N₂H⁺, where we marginally resolve the HFS lines and see potentially different velocity components that we do not aim to decompose here and hence used a simplified fitting approach using a single Gaussian. The physical properties of the N₂H⁺ emitting gas and the average line width have considerable uncertainties, therefore we do not use the results obtained for this species in Sect. 4.

Fig. 5 Examples of the different line profiles observed in the spectral survey. From the top to the bottom, we show in the left column a narrow line profile for HC18O+ (3–2) and HC3N (25–24) with a two-component Gaussian fit and HC3N (33–32) with a single broad component. In the right column,HCO+ (3–2) has a line profile with redshifted self-absorption, SiO (4–3) has a broad component coming from the outflow, and the profile of CH3OH (11−5–12−4) shows two velocity components tracing the accretion shocks. All the lines are plotted in a window of 70 km s−1 centred on the rest velocity. The dotted red line shows the systemic velocity of the source at –43.13 km s−1.

3.3 Rotational diagrams

The large frequency coverage of the survey gives access to a large number of rotational transitions for the different species. Assuming LTE conditions, this allows us to estimate the total column density, N, of each species and the rotational temperature, T_rot. Therefore, in order to obtain a first estimate of the physical and chemical conditions of the emitting gas, we used rotational diagrams to estimate N and T_rot that then served as an input for a more detailed LTE modelling using Weeds following a similar approach as described in (Belloche et al. 2016) (Sect. 3.4).

To do this, we use the integrated intensity, W, of each unblended line at a frequency, ν, from the Gaussian fitting. The population of each rotational level of a molecule is given by $$N_u^{\text{thin}}=\frac{8\pi k{v}^{2}W}{h{c}^3{A}_{\text{ul}}{B}_{\text{dil}}}$$(1)

for optically thin lines, where k is the Boltzmann constant, B_dil is the beam dilution factor depending on the size of the emitting region, θ_S, and the beam size, θ_beam, A_ul is the Einstein coefficient, h is the Planck constant, and c is the speed of light (see e.g. Goldsmith & Langer 1999). Because we assumed that most of the emission is unresolved within the beam, the beam dilution factor, defined as $B_{\text{dil}}=\frac{{\theta }_{\text{S}}^2}{{\theta }_{\text{S}}^2+{\theta }_{\text{beam}}^2}$, needs to be taken into account and shows a large variation due to the large frequency coverage of our dataset. The size of the emitting region is determined based on the rotational diagrams and the LTE modelling using Weeds in an iterative process (see Sect. 3.4).

In the case of optically thick lines, the line opacity needs to be taken into account, and the population of each level becomes $N_u=N_u^{\text{thin}}{C}_{\tau }$ where $C_{\tau }=\frac{\tau }{1-e^{-\tau }}$ The parameter τ corresponds to the total line opacity, from which C_τ, the optical depth correction factor is computed. Following (Goldsmith & Langer 1999), to estimate N_tot, we used the expression $$ln\left(\frac{N_{\text{u}}}{g_{\text{u}}}\right)=ln\left(\frac{N_{\text{tot}}}{Z}\right)-\frac{E_{\text{up}}}{k{T}_{\text{rot}}},$$(2)

where Ζ is the rotational partition function, g_u is the statistical weight, and E_up is the upper-level energy of each transition.

The rotational diagrams are then displayed as ln$\left(\frac{N_{\text{u}}}{g_{\text{u}}}\right)$ versus E_up, and by fitting a linear function, we derive the total column density N_tot and the rotational temperature T_rot. The partition function was interpolated from the tables given in the CDMS and JPL databases. We perform a linear least-square fit to the rotational diagrams, and from these parameters, we computed the rotational temperature and the total column density of each temperature component of the molecule. We estimated uncertainties on the rotational temperature and the column density using a Monte Carlo approach. We varied each measurement in the rotational diagram within its error bars assuming a uniform distribution and repeated the least-square fit. The final results correspond to the average of the estimated column densities and rotational temperatures, while their uncertainties correspond to the standard deviation of the values determined with our Monte Carlo method.

Several of the molecules we cover exist in different symmetries. We treat them together as single species, such as the forms E and A of propyne (CH₃CCH), methyl cyanide (CH₃CN), and methanol (CH₃OH).

We detect several transitions for C₂H₃CN, but all except four lines are found to be blended. Therefore, we only used these lines for the rotational diagram covering the range between 160 Κ and 230 K, which represents half of the complete range of detected transitions for C₂H₃CN (160–350 K).

As an example, we show the rotational diagram of the HC₃N molecule in Fig. 7, for which the rotational temperature increases from 74 Κ to 78 K, corrected for optical depth effects. We show all rotational diagrams used for this study in Appendix B.

Fig. 6 Hyperfine structure fitting of the transition CCH (2–1) in the 2 mm band. The result of the fitting procedure is indicated in blue, and the systemic velocity corresponds to the frequency of the most intense component hyperfine transition. The quantum numbers follow the notation from (Gottlieb et al. 1983), with Ν corresponding to the rotational levels, J to the spin doubling, and F to the hyperfine levels.

Fig. 7 Rotational diagram of HC3N for the broad velocity component towards G328.2551-0.5321. The rotational diagram when we assume that HC3N lines are optically thin is shown in blue. We represent the rotational diagram corrected for the opacity of each line in red. Rotational temperatures and column densities are indicated with the same colour code. The beam dilution was calculated for each transition for a source size of 4″. The error bars represent a lσ- statistical noise plus a 20% calibration error. To calculate the uncertainties on the rotational temperature and the column density, we used the Monte Carlo method and assumed a uniform distribution of the error for each data point.

3.4 LTE modelling

3.4.1 Method

In order to fit the observed spectrum, we performed an LTE modelling with Weeds (Maret et al. 2011). The input parameters for each molecule were the column density (N_tot [cm⁻²]), the size of the emission (θ_S [″]), the excitation temperature (T_ex [K])⁸, the averaged full-width at half maximum (fWhM, Δυ [km s⁻¹]), and the velocity offset (υ_off [km s⁻¹]) with respect to the systemic velocity of the source.

For the light molecules, we first used the excitation temperature and the column density from their rotational diagram as an input parameter for the modelling. For a first guess of the sizes, we used the ALMA images shown in (Csengeri et al. 2019), which were available for a limited number of molecules. For the remaining molecules, we estimated the sizes where possible based on the iterative approach described below. We used the averaged FWHM and velocity offset from the Gaussian line-fitting of the unblended transitions. Using these first-guess parameters, we then visually inspected the obtained model and extracted the residual. When the residual was lower than the 3σ threshold, we considered the fit to be satisfactory, otherwise, we modified the input parameters and repeated the fit. To improve the model, we varied the column density and the temperature within the error bars given by the fit to the rotational diagrams. When we failed to obtain a satisfactory fit within this parameter range, we varied the size of the emitting region in the rotational diagrams and used the new column density and temperature as input. With this new set of parameters (column density, temperature, and size), we repeated the process until we found a satisfactory fit within the 3σ threshold for all the transitions. Our uncertainties depend on the spectroscopic data of each molecule and on the uncertainties due to the source size. For molecules that transition at similar upper-level energies over the spectral band at different frequencies, for example, CH₃OCH₃ or SO₂, the source size is relatively well constrained. For linear molecules with increasing upper-level energies as a function of frequency, our size estimates are limited to those with optically thin transitions. Molecules with highly optically thick emission in all their transitions, such as HCO⁺, HCN, HNC, H₂S (cold phase), and CS, could not be modelled by this method, hence we used their optically thin isotopologues assuming that they originate from the same gas.

In all cases, the rare isotopologues of each main species were modelled separately in order to determine their column density independently from the main isotopologue, allowing us to derive isotopic ratios (see Sect. 5.3). Because all SiO transitions show a complex velocity profile, we used the canonical ²⁸Si/²⁹Si ratio (19.6; Penzias 1981a) to estimate the SiO column density.

In order to identify all the transitions for each molecule, we first initiated a model with a first guess of the parameters for the complex molecules, identified the unblended transitions, and derived a rotational diagram. We then obtained the LTE model and corrected the list of the detected transitions if necessary in an iterative fashion. We applied this method to all the COMs.

The different envelope components identified by their kinematics (i.e. velocity profile) were treated independently at first and then modelled together to inspect the final model of the molecule. This approach assumes a uniform spatial distribution of the molecules and ignores their spatial segregation, as suggested by ALMA observations (Csengeri et al. 2019), on scales smaller than our angular resolution.

After we obtained a satisfactory fit within our quantitative criterion (unless stated otherwise, see Sect. 3.4.2) for each molecule separately, we modelled all the species together to determine whether some remaining lines were left and if the blended lines were properly modelled. Where necessary, we adjusted the model of each molecule. For each species, we were able to determine the column density, the excitation temperature, and the size of the emission for each part of the envelope from which the emission arises. The results of the modelling are summarised in Table 4.

3.4.2 Exceptions

For molecules with K-ladder transitions (H₂CO, H₂CS, and CH₃CN), we used a different approach. For (ortho and para) H₂CO and H₂CS, we estimated the temperature following the method of (Mangum & Wootten 1993) using the ratio between the $J=3_{2,1}-2_{2,0}\left(E_{\text{up}}/k=68\text{K}\right)$ and $J=3_{0,3}-2_{0,2}\left(E_{\text{up}}/k=21\text{K}\right)$ transitions, and find a temperature of 89 K, corresponding to the warm component of the envelope. The estimated column density has a larger uncertainty, because the emission sizes are poorly constrained. Therefore, we estimated a column density using sizes of 4–7″ and 2–5″, which is the ortho/para-H₂CO and ortho/para H₂CS, respectively, which correspond to the typical sizes for the warm component.

We find, however, that this component does not fit all the K transitions well, suggesting the presence of another, cold component for both H₂CO and H₂CS. We estimated the temperature of this component based on the rotational diagram of the K = 2 transitions, and similarly then used the LTE models to fit the emission. We took 2O–25″ and 30–35″ for the size of the emitting region for these components, respectively. Similarly as for the warm component, the uncertainty on the size of the emitting region and the column density for these species is considerably. Consequently, we do not use these values for the discussion in Sect. 4.

For CH₃CN, we estimated the temperature, column density, and size of the emission using the rotational diagrams for the K=0 component. We find a moderately extended size of 15″ and a temperature of 40 K. Our LTE model is not satisfactory, however, suggesting the need of a warm component. Because the APEX data do not give further constraints on this warm component, we simply added a 200 K temperature component with a size of 1”, which is constrained by the ALMA data in (Csengeri et al. 2019). This composite model fits the K = 1–6 transitions rather well and underestimates the higher excitation K = 7,8, 9 transitions in the J = 17–16 and J = 19–18 lines. An additional hot component could explain this excess emission, but we cannot constrain its properties. The ¹³CH₃CN isotopologue is not detected in our survey.

Although the rotational diagram of C₂H₅CN shows some scatter (Fig. B.9), we fit only one temperature component, which is found to have a temperature around 90 K. Because the ALMA observations (Csengeri et al. 2019) suggest the presence of a compact warm component, we tested whether our LTE fitting using two temperature components (25 K and 200 K) instead of a single one (90 K) would provide a better fit. We found that we cannot distinguish between these models; both of them could be consistent with our data. For simplicity, we therefore used the values of a single-component fit.

H¹³CN, HC¹⁵N, and HN¹³C have only three transitions in our band, and the rotational diagram is consistent with a single gas temperature component. However, the LTE model reveals that one temperature component is not satisfactory. Therefore, we need to use a model with two temperature components: one at 10 K with a size of ~35″, and a second one at ~80K with a size of 3″.

The rotational diagrams of HDO and HC(0)NH₂ show only one temperature component, but the line profiles reveal the dynamics of accretion shocks similarly to CH3OH (Fig. 5). We added a shock component to our model using the excitation temperature and size of the emitting region of CH₃OH,v_t = 1.

3.5 Outflow contribution

The line profiles of several molecules show high velocity emission. The most pronounced, highest velocity wings are seen in CO, but SiO, H₂S, HCO⁺ HNC, HCN, and CS also indicate non-Gaussian line-wings that are likely associated with the high-velocity outflowing gas. A prominent bipolar molecular outflow has been imaged in the CO J = 3–2 line by ALMA (Csengeri et al. 2018). Owing to their high optical depth, many of these molecules exhibit self-absorbed line profiles. For these molecules, the column densities and the temperature were determined from the rarer isotopologues. It is beyond the scope of this work to estimate the properties of the gas associated with the outflowing material. In the following, we therefore omit the outflow components from the LTE modelling and fit a single Gaussian to the line to take only the emission from the envelope into account.

4 Structure of the G328.2551-0.5321 protostar

The protostar embedded in the G328.2551-0.5321 clump has been shown to be dominated by a single massive fragment down to a size scale of ~400 au (Csengeri et al. 2018). Therefore, we assumed that all the emission detected towards this source originates from a single protostellar envelope that is spatially unresolved within the beam sizes of our observations (0.2–0.5 pc at the distance of the source).

4.1 Observational constraints

In Sect. 3.2, we identified broad and narrow line profiles, outflow wings, and for some species, multiple (typically two) velocity components. This suggests that although it remains spatially unresolved in our observations, the protostellar envelope has an internal structure in which the physical parameters may change.

From the LTE fitting of simple species (Sect. 3.4), we derived the distribution of measured line widths versus T_ex (Fig. 8). We show the results of a K-mems clustering analysis (MacQueen 1967), which shows a preference for two groups of line widths corresponding to a low-temperature and more quiescent gas (T_ex < 50 K, Δv < 6 km s⁻¹), and a warmer component with typically broader line widths (50 K< T_ex < 150 K, 6 km s⁻¹ < Δv < l2kms~¹). For simplicity, we do not include in Fig. 8 the contribution from the accretion shocks, which correspond to high T_ex > 120 Κ and exhibit Δv < 3–1 km s″¹. This is discussed in more detail in Sect. 6.1. The broader velocity dispersion close to the protostar can be explained by a combination of larger thermal broadening, a higher level of micro-turbulence due to the mechanical feedback of the protostellar embryo, and to other kinematic effects such as infall and rotation. The contribution of thermal broadening to the line width is estimated by ${\sigma }_{\text{th}}=\sqrt{k_{B}T/\left(\mu m_{\text{H}}\right)}$, where k_b is the Boltzmann constant, Τ is the gas temperature, μ is the molecular mass, and oth is the mass of the hydrogen atom. This gives a thermal velocity dispersion of 0.07 km s⁻¹ for the second lightest molecule, CCH, in the cold gas and 0.17 km s⁻¹ for HCN in the warm gas. When σ_th is converted into the FWHM of a Gaussian line using $FWHM=2\sqrt{2ln2}{\sigma }_{\text{th}}$, the thermal broadening is still below the 0.7 km s⁻¹ velocity resolution of the final data product. The thermal broadening has only a minor effect in this case. We do not discuss the relative contribution of these effects here.

The LTE fitting also provides the distribution of the obtained T_ex as a function of the size of the emitting region (Fig. 9). Although here (T_ex) corresponds to a spatially averaged excitation temperature over the size of the emission (θ_s), it is expected to represent the temperature structure of the source because under LTE conditions, T_ex is equivalent to the kinetic temperature. We find that even this spatially averaged excitation temperature exhibits a rather gradual decrease in temperature with the radius. We fit a power-law distribution of $\langle T_{\text{ex}}\rangle \left({\theta }_{\text{s}}\right)\propto {\theta }_{\text{s}}^{\beta }$ , and find aβ of –0.6. This result is intriguing because we did not consider any a priori temperature structure in our line fitting and decomposition analysis.

Although we observe a rather gradual decrease in the temperature, we can distinguish two physical components for the bulk emission of the envelope based on the line-width analysis. Therefore, we discuss in the following the warm component of the envelope (together with the accretion shocks) within a typical size of 0.06 pc and T_ex > 50 Κ and the cold component of the envelope, which extends up to a size of 0.5 pc and exhibits T_ex <50K.

Fig. 8 Excitation temperature vs. mean line width. A statistical study of K-Means (MacQueen 1967) was performed to determine the limit between the different components of the envelope. The blue data points represent the cold component, and the warm component of the envelope is represented by the red data points. Single transitions fitted with multiple velocity components are considered as individual measurements in this figure.

Fig. 9 Excitation temperature vs. the size of the emission region for each molecule whose emission arises from the envelope. The black curve represents the power-law fit of the data. The grey area indicates power-law temperature profiles with an exponent between –0.3 and –0.9. This profile does not include the emission originating from the accretion shocks that correspond to a local increase in temperature. It also excludes species for which no strong constraints on the size of the emitting region could be derived (see Sect. 3.4.2).

Fig. 10 Column density profile (dashed black line) and averaged column density (blue line) vs. radius (r). The blue triangle shows the averaged column density estimate from ACA array observations from Csengeri et al. (2017a). The green circle corresponds to the beam-averaged column density estimated from ATLASGAL (Csengeri et al. 2014). The yellow square shows the column density estimates roughly at the position of the accretion shocks with ALMA (Csengeri et al. 2019). The error bars indicate a 20% error on these measurements.

4.2 Toy model for the source structure

To constrain the physical origin and the abundance variations of the molecules detected in this survey, we used a toy model (see Fig. 10) to describe the physical components within the envelope. (Csengeri et al. 2018, 2019) showed that in the close vicinity of the protostellar embryo, a picture emerges that is qualitatively similar to that observed for some low-mass protostars, with a compact accretion disk and shocks at the centrifugal barrier. These components close to the vicinity of the protostar are confined to a compact region of <1″that is considerably smaller than our beam, hence we did not model it in detail. For simplicity, we also omit here the discussion of the high-velocity emission originating from the outflowing gas and fast shocks, and instead focus our discussion on the bulk of the envelope and its physical properties.

Csengeri et al. (2018) reported that the source structure can be fitted with a power-law profile and a compact source, the latter likely corresponding to an accretion disk. Such a profile is consistent with the results of Csengeri et al. (2019), who reported similar values of H₂ column densities at three positions with a projected distance of ≤ 1000 au from the protostar, suggesting a flattened inner region. Therefore, we describe the H₂ volume density as a function of radius using a Plummer-like function, which corresponds to a flattened density profile towards the innermost regions, and a power-law profile in the outer regions. We used a Plummer-like column density profile of $$N_{{\text{H}}_2}\left(r\right)=N_0\times {\left(1+r^2/a^2\right)}^{-1/2}$$(3)

(Fig. 10), where a and N₀ are the Plummer radius and the normalisation factor (respectively). We chose a = 750 au, and computed N₀ such that the observed $N_{{\text{H}}_2}$ at radius a is consistent with the observed value of $N_{{\text{H}}_2}\left(a\right)=1.9\times {10}^{24}{\text{cm}}^{-2}$ from Csengeri et al. (2019). This corresponds to a density profile following the form of n∝ r⁻² at r ≥ a. The molecular column density estimates provided by the LTE modelling with Weeds correspond to a certain size of the emitting region (θ_s). To be able to infer molecular abundances relative to H₂, we therefore computed the average H₂ column density $\left(\langle N_{{\text{H}}_2}\rangle \right)$ for this corresponding size. We computed an average H₂ column density,$\langle N_{{\text{H}}_2}\rangle$, as a function of physical (and angular) radius, r, and θ_s, respectively, as shown in Fig. 10. Weeds assumes a Gaussian brightness distribution for the source, where the size of the emitting region, θ_s, corresponds to the FWHM of the Gaussian (Maret et al. 2011). We related our source radius, r to θ_s using r = θ_s/235. In Table 4, we indicate the source size, θ_s for which the averaged H₂ column density is calculated.

In Fig. 10, we compare the average H₂ column density of our model to the estimates of the H₂ column density from APEX/LABOCA, and ALMA data from the literature. We used a dust emissivity of 0.0185 cm² g⁻¹ at 345 GHz (which already includes a gas-to-dust ratio of 100) to compute the H₂ column density from the flux density using the standard formulae. Taking the beam-averaged flux density values from the ATLASGAL survey (Csengeri et al. 2014) with the temperature of ~22K estimated in Csengeri et al. (2018) for the cold gas component, we estimate a beam averaged column density of $N_{{\text{H}}_2}=1.8\times {10}^{23}{\text{cm}}^{-2}$ within the APEX beam at 345 GHz of ${19}_{.}^{\text{'}\text{'}}2$. We performed the same calculations taking the ALMA 7m array measurements from Csengeri et al. (2017a) and estimate an averaged column density of $N_{{\text{H}}_2}=1.3\times {10}^{24}{\text{cm}}^{-2}$ within the ALMA 7 m synthesised beam. These estimates are broadly consistent with the $N_{{\text{H}}_2}$ profile we used.

5 Molecular composition: Simple molecules

This unbiased spectral line survey allows us to detect 29 species, 54 including isotopologues, which are simple molecules, that is, have fewer than six atoms. We associate their origin with the physical components of the envelope derived above and investigate the potential chemical origin of these species. In particular, we discuss here the properties of carbon-chain molecules, O-,N-,and S-bearing molecules, and deuterated species in the cold and warm components of the envelope. The sketch of the molecular composition of the various physical components of the envelope is shown in Fig. 11.

5.1 Sulphur chemistry

5.1.1 Richness in sulphur-bearing species

One of the most remarkable results concerning the molecular composition of G328.2551-0.5321 is the detection of a large number of sulphur-bearing molecules (10 species plus eight isotopologues) with high column densities and abundances. In addition to the most abundant S-bearing molecules in star-forming regions, such as SO, SO₂, OCS, H₂S, CS, and H₂CS (e.g Hatchell et al. 1998; Li et al. 2015; Widicus Weaver et al. 2017), we detect here several transitions from SO⁺, HCS⁺, NS, and CH₃SH and show a tentative detection of NS⁺. We list the estimated parameters (T_kin,N, and θ_s) from the LTE modelling (Sect. 3.4) in Table4 and show the column densities and abundances of these species in the cold and warm components of the envelope in Fig. 12.

The undepleted sulphur abundance in the diffuse interstellar medium has been found to be around $S/H\sim \left(3.5\pm 1.5\right)\times {10}^{-6}$ (Goicoechea et al. 2006). In dark clouds and hot cores, however, the total abundance of S-bearing molecules represents only 0.1% of the cosmic value (Tieftrunk et al. 1994; Charnley 1997). It might be argued that sulphur is present in its atomic form and thus cannot be observed. However, this assumption does not hold in hot core models (Vidal et al. 2017). Towards G328.2551-0.5321, the highest molecular column densities are found to be up to 10¹⁷ cm⁻² in the warm component of the envelope, and abundances relative to H₂ as high as 8 × 10⁻⁸ for OCS and SO₂ are reached. In the cold gas, the highest column density is found for H₂S with a value of 2.6 × 10¹⁶ cm⁻² and an abundance relative to H₂ of 4.6 × 10⁻⁸. These high molecular abundances of S-bearing species correspond to ≥2% of the cosmic value in the warm component of the envelope and to ≥1% of the cosmic value in the cold component of the envelope. Therefore, this object has a less pronounced sulphur depletion than is commonly observed in the dense cold gas. However, high abundances of H₂S have been observed towards some high-mass star-forming regions, such as the Orion hot core (Blake et al. 1996), and high values of SO and SO₂ were observed in W3 IRS5 (Wang et al. 2013).

Similar to other species, the sulphur-bearing species show different line profiles, suggesting that these molecules are present in various physical components of the envelope. For example, a line profile with a redshifted self-absorption suggesting infall motions (e.g. Myers et al. 1995) is observed in the ground-state transition of ortho-H₂S. Broad, non-Gaussian line-wings suggest that CS, SO, and SO₂ are also present in the outflowing gas. We associate the Gaussian, relatively narrow velocity components to the bulk of the envelope.

Based on the kinetic temperatures estimated from the LTE modelling, we distinguish between an origin in the cold or warm components of the envelope. For temperatures below 70 K, sulphur-bearing molecules might desorb through non-thermal processes, and gas-phase processes dominate the chemical reactions (van der Tak et al. 2003; Wakelam et al. 2011). We find that this temperature corresponds rather well to the separation between the cold and warm components of the envelope. The largest variety of sulphuretted molecules (nine, plus eight isotopologues) is detected towards the cold component of the envelope, while in the warm component, we identify only five molecules (Fig. 13). In general, we find higher column densities and molecular abundances in the warm component than in the cold gas, which is consistent with the expected more efficient thermal desorption at higher temperatures. It is very intriguing that S-bearing molecules such as OCS, H₂S, SO₂, CS, and SO exhibit abundances as high as methanol, suggesting a particular enhancement of sulphur chemistry in this object. In particular, methanol is only one order of magnitude more abundant than methanethiol with a ratio CH3OH/CH3SH of 25, while in the well-studied hot core Sgr B2(N2) (Müller et al. 2016), the methanol abundance relative to H₂ is two orders of magnitude higher and the CH₃OH/CH₃SH ratio is one order of magnitude higher. Towards another classical hot core associated with G327.3–0.6, the CH₃OH/CH₃SH ratio also appears higher, with a value of 60 (Gibb et al. 2000). In the following, we therefore investigate the sulphur-bearing species in the cold and warm components of the envelope in more detail.

Fig. 11 Sketch of G328.2551-0.5321 ignoring projection effects. The full source structure is represented on the left, and we zoom into the envelope of the protostar on the right.

Fig. 12 Column densities (upper panel), abundances relative to H2 (middle panel), and relative column densities to methanol (lower panel) of the light molecules and methanol. Only the main isotopologues are represented.

Fig. 13 Abundances of the S-bearing molecules relative to H2 in the cold and warm phases. Each molecule is represented by one colour.

5.1.2 Sulphuretted species in the cold component of the envelope

The velocity profiles suggest that both the S-bearing ions and the neutral molecules originate from the same gas that is associated with the cold quiescent envelope. Other studies based on mapping of sulphuretted molecules suggest, however, that despite a similar line width compared to quiescent gas tracers (e.g. H¹³CO⁺ c-C₃H₂, and HDCO), for example, SO and SO₂ do not necessarily originate from the undisturbed gas (e.g. Chernin et al. 1994; Bachiller et al. 2001; Buckle & Fuller 2002). Although these species may chemically form purely via gas-phase reactions in the undisturbed gas (e.g. Pratap et al. 1997; Dickens et al. 2000), we cannot exclude a contribution from the outflow cavity walls and thus from the disturbed gas (Codella et al. 2013) even for the narrow line-width component. The excitation temperatures in the cold component of the envelope are all around 10–18 Κ for sulphur-bearing molecules, and we find a cold-gas component at a temperature of 10 Κ traced by SO, ³⁴SO, and NS. Most of the molecules, such as H₂CS, CS, HCS⁺, and SO⁺, show a temperature of 25–30 K, while we also trace a warmer component at about ~40K with SO, H₂S, SO₂, and OCS.

The highest column density in the cold gas is observed for ${\text{H}}_2\text{S}\left(2.6\times {10}^{16}{\text{cm}}^{-2}\right)$ and corresponds to a molecular abundance relative to ${\text{H}}_2\text{\hspace{0.17em}of\hspace{0.17em}}4\text{.6}\times {\text{10}}^{-8}$. This is a significantly higher abundance than the median values calculated for hot cores and hot corinos (e.g. Li et al. 2015), suggesting that a large fraction of sulphur is released into the gas phase here. It is, however, close to the value reported towards Orion KL (Blake et al. 1996; Esplugues et al. 2014; Crockett et al. 2014). The second most abundant species are CS and OCS, followed by SO and SO₂, while the lowest abundances are at the order of a few times 10⁻¹⁰, corresponding to molecular ions, such as HCS⁺ and SO⁺. The abundances of S-bearing molecules found in the cold component of G328.2551-0.5321 are higher than the median values found in hot cores and hot corinos (e.g. Charnley 1997; Hatchell et al. 1999; Li et al. 2015).

5.1.3 Evidence for shock chemistry in the cold gas?

We also detect the reactive molecular ion SO⁺ in the cold gas phase with an abundance relative to H₂ of 1.9 × 10⁻¹⁰. This molecule has been detected in the interstellar medium in various environments (Turner 1992, 1994); for example towards photodissociation regions (PDRs; Fuente et al. 2003), and star-forming regions as well (Podio et al. 2014; Nagy et al. 2015). Towards G328.2551-0.5321, we do not detect other species that would be characteristic of PDR chemistry, however, such as HOC⁺, or CO⁺ suggesting that the origin of SO⁺ could be rather related to ion-neutral chemistry or dissociative shocks (Turner 1994). Further evidence for this scenario is the abundance ratio of X(SO⁺)/X(SO₂)~0.05, which is below the value of 0.4–1 measured for PDRs (Fuente et al. 2003; Ginard et al. 2012), and is similar to values determined for other high-mass star-forming regions (Nagy et al. 2015). Gas-phase chemical models predict a high X(SO⁺)/X(SO) ratio of ~l for dense and ionised regions (Ginard et al. 2012). Our value of 0.04 is similar to the 0.06–1 found towards MonR2 (Ginard et al. 2012) and cold dense cores (Agúndez et al. 2019) and translucent clouds (Turner 1995, 1996). Neufeld & Dalgarno (1989) have shown that the release and ionisation of sulphur into the gas phase through shocks leads to an enhancement of SO⁺ with abundances relative to H₂ reaching 10⁻¹⁰–10⁻⁹. Shock chemistry could therefore explain the higher abundance of SO⁺ towards G328.2551-0.5321.

Similarly, HCS⁺, H₂CS, and CS have abundances of about 10⁻¹⁰, 10⁻⁹, and 10⁻⁸, respectively, which are higher than what is observed around high-mass protostars (Hatchell et al. 1998; van der Tak et al. 2003; Li et al. 2015) and are similar to what is observed in shocked regions (Minh et al. 1991; Turner 1992, 1994; Podio et al. 2014). Our observations of sulphur-bearing species in the cold component of the envelope therefore suggest a potentially significant impact or contribution from shock chemistry that leads to enhanced abundances of sulphur-bearing molecules. While the exact origin of shocks is difficult to elucidate because a powerful ejection of material in the form of an outflow has been observed towards this object (Csengeri et al. 2018), the outflow cavity walls, where the velocity profile resembles quiescent gas, might be a possible origin.

5.1.4 Sulphur in the warm envelope

In contrast to the cold component of the envelope, the warm gas shows less diversity in sulphur-bearing molecules. We detect SO₂, OCS, H₂S, H₂CS, and CH₃SH in decreasing abundance in the warm component of the envelope. The abundances of OCS, SO₂, H₂CS, and CH₃SH relative to H₂ are higher in the warm (Fig. 13, upper panel) than in the cold component of the envelope, reaching values as high as 7–8 x 10⁻⁸ for SO₂ and OCS, which are among the most abundant sulphur-bearing molecules in star-forming regions (e.g. Sakai et al. 2010). Towards typical hot cores, abundances of S-bearing molecules of about 10⁻⁸ have been reported for SO₂, OCS, and CS and abundances of 10⁻⁹ for H₂CS (e.g. Hatchell et al. 1998; van der Tak et al. 2003; Herpin et al. 2009). The excitation temperatures are similarly high for these molecules, around 100 K, except for H₂S, which is at 166 K.

Towards Orion KL, the carbon-free sulphur-bearing species show an order-of-magnitude increase in molecular abundances towards the hot regions (Luo et al. 2019). In our case, however, the carbon-sulphur bearing molecule OCS also shows a significant enhancement in molecular abundance towards the warm region. This prevents us from making a similar distinction between the behaviour of carbon-free and carbon-bearing sulphuretted species.

While the main sulphur reservoir on the grains is still poorly known, several candidates have been proposed, such as H₂S (Charnley 1997; Holdship et al. 2016; Vidal et al. 2017), SH (Vidal et al. 2017), and/or OCS (Hatchell et al. 1998; van der Tak et al. 2003). On the other hand, SO, SO2, and other simple S-bearing molecules are thought to be formed in the gas phase (Podio et al. 2014; Esplugues et al. 2014; Holdship et al. 2016). OCS and SO₂ have been firmly or tentatively detected in interstellar ices (Palumbo et al. 1995; Boogert et al. 1997), but not H2S. The high H₂S abundances observed in hot cores (Hatchell et al. 1998) and in the present object are inconsistent with the gas-phase formation. Gas-phase models can reproduce observed abundances of ~10⁻¹⁰ in the quiescent cold dark cloud TMC-1 (Millar et al. 1989). Charnley (1997) and Millar et al. (1997) suggested that in hot cores, H₂S must be produced in significant amount in the grain mantles. When it is released into the cold gas-phase component of the envelope, it leads to the formation of SO and SO₂, and ultimately, to CS and H₂CS. The molecules OCS and SO₂ represent the most important reservoir of sulphur among the detected molecules in the warm gas phase, corresponding to 40% each (Fig. 13, upper panel). H₂S is the main reservoir in the cold gas phase (~50%), while it represents only 15% in the warm gas phase. Chen et al. (2015) reported that uV radiation could open a chemical reaction channel that allows the formation of OCS from H₂S, which could explain the destruction of H₂S in the warm component of the envelope. The weak UV radiation from the protostar could influence the warm component of the envelope without having a significant impact on the cold component of the envelope. Another possible explanation is a reaction on the grain surface of H₂S to form ethanethiol C₂H₅SH), which was proposed by Müller et al. (2016). The authors claimed, however, that this reaction has only a minor effect on the H₂S abundance. This product is not detected here with an upper limit of 2.7 × 10⁻⁸ on its abundance relative to H₂. van der Tak et al. (2003) argued that OCS is the main sulphur reservoir on the grains, which could naturally provide an explanation for the high abundances of OCS observed in the warm component of the envelope.

5.1.5 Role of S-bearing species as chemical clocks

Chemical models predict that the relative abundances of sulphur-bearing species may be used to estimate the ages based on the chemistry of dense clouds (Buckle & Fuller 2003). For this purpose, the abundance ratios of several sulphur-bearing species, for example, OCS/SO₂, SO/SO₂, and H₂S/SO₂, have been invoked as tracing the chemical age of the gas (e.g. Hatchell et al. 1998; Herpin et al. 2009; Wakelam et al. 2011). Here we investigate these ratios for molecules that we find in the warm component of the envelope, that is, OCS/SO₂ and H₂S/SO₂. The models from Wakelam et al. (2011) predict a decrease in these ratios for sources with more evolved chemistry. In Fig. 14, we compare our measurements with the values of Herpin et al. (2009), which were estimated for three sources that are located within a heliocentric distance of 5 kpc. The authors proposed a rough evolutionary sequence from the youngest to the oldest for these sources mainly based on the flux ratio of the hot to cold gas in the spectral energy distribution (SED) of this sample, as of IRAS 18264–1152 (infrared quiet dense core with a dynamical age of ~0.5 × 10⁴ yr), IRAS 05358+3543 (massive dense core associated with an outflow and a dynamical age of ~3.6 x 10⁴ yr), and IRAS 18162–2048 (infrared bright dense core with a powerful outflow and dynamical age of ~10⁶ yr). This comparison shows moderate variation in the OCS/SO₂ ratio, but the ratio of H₂S over SO₂ is clearly higher towards G328.2551-0.5321, which is consistent with a scenario in which it is younger than the other sources.

5.2 Deuterated molecules

Deuterium enhancement (deuteration) occurs in cold, dense interstellar clouds, that is, preferentially in prestellar cores (Aikawa et al. 2012; Taquet et al. 2012, 2014). It has been invoked as a potential age tracer for pre- and protostellar evolution.

Towards G328.2551-0.5321, we do not see many deuter-ated molecules in general, only four deuterated molecules are firmly detected in the cold component of the envelope (DCN, DCO⁺, HDCO, and HDCS), while in the warm component of the envelope and in the accretion shocks, only HDO is detected. We list their estimated parameters based on LTE modelling (Sect. 3.2) in Table4 and discuss upper limits for frequently detected deuterated molecules, such as HDS and N₂D⁺, and deuterated methanol in the following.

We derived the deuterium fractionation (D/H ratio) based on the abundance ratios (Χ_χD/Χ_χη) of the species for which the non-deuterated form (or one of its rarer isotopologues) is also detected. When the chemical group containing the deuterium atom is composed of several hydrogen atoms, we need to take the different arrangements between Η and D within the group into account in order to infer the D/H ratio from the abundance ratios (Χ_χD/Χ_χη) (for more details, see Manigand et al. 2019). In our case, this concerns HDO, HDCO, HDCS, and the upper limits estimated for HDS, where D/H ratio is defined by (1/2) × Χ_χD/Χ_χη- For the upper-limit estimates of CH₂DOH, we need to use $\text{D}/\text{H}=\left(1/3\right)\times X_{{\text{CH}}_2\text{DOH}}/X_{{\text{CH}}_3\text{OH}}$ In Table 5 we list the two [XD]/[XH] abundance ratios and the deuterium fractionation (D/H ratio). In the following, we discuss the deuterium fractionation as the D/H ratio.

We find a deuterium fractionation between 0.002 and 0.03 for the detected deuterated molecules in the cold component of the envelope. The strongest deuterium enhancement of 0.03 is seen for HDCS, while DCN is at 0.005–0.01. DCO⁺ and HDCO exhibit a deuterium fractionation of 0.002–0.003 and 0.02, respectively. The cosmic deuterium abundance is 1.5 × 10⁻⁵ (Linsky 2003), hence the deuterium enhancement in general is significant towards this object, although it is typically lower than the values of 0.01–0.3 found for low-mass star-forming regions (e.g. Turner 2001; Marcelino et al. 2005).

Deuterium fractionation is efficient in the interstellar medium. The larger proton affinity of D compared to Η leads to an increase in the fraction of D in neutrals and ions, but the zero-point energy for D containing molecules is also lower than that of their main isotopologue. One of the main reaction pathways is through the reaction with H₂D⁺ (Millar et al. 1989), which produces for example DCO⁺ and N₂D⁺ (Aikawa et al. 2018). This reaction is efficient at low temperatures. An alternative route for deuteration to proceed is efficient even at higher temperatures (>60 K) in the gas phase by reaction with CH₂D⁺ (Turner 2001; Roueff et al. 2007; Parise et al. 2009). This chemical pathway indicates that even after the sublimation of CO, deuterium chemistry may remain active. For example both DCN and DCO⁺ could be formed in an efficient way in the warm component of the envelope (Parise et al. 2009), and they indeed show a high deuteration ratio in the warm gas phase towards hot cores and disks (Ren et al. 2012; Zinchenko et al. 2012; Favre et al. 2015). Towards G328.2551-0.5321, both DCN and DCO⁺ are only detected in the cold component of the envelope at low temperatures below 20 K, favouring a pathway through reaction starting from H₂D⁺.

Similarly, the singly deuterated thioformaldehyde, HDCS, could also form in the gas phase and might even reach a high fractionation of >0.02 according to the models of Minowa et al. (1997) and Marcelino et al. (2005). This is in fact consistent with our observational results. These models also give a similarly high fractionation for HDCO, where our observational results suggest a somewhat higher fractionation of 0.04.

The only deuterated molecule that we identified in the warm component of the envelope is HDO. While we did not detect ${\text{H}}_2{}^{18}\text{O}$, which would allow an estimate of the H₂O abundance, our modelling allows us to estimate an upper limit of 1.0 × 10¹⁸ cm⁻² for ${\text{H}}_2{}^{18}\text{O}$ in the warm component of the envelope assuming the same size, excitation temperature, and line width as HDO. The canonical ^l6O/¹⁸O ratio of 500 implies a water deuteration of >0.0002 towards G328.2551-0.5321. Towards hot cores, this abundance ratio is typically found to be between 5 × 10⁻⁵ and 0.0003 (Jacq et al. 1990; Comito et al. 2003; Liu et al. 2013; Coutens et al. 2014; Kulczak-Jastrzębska 2017) in the hot envelope and up to 0.001 for Orion-KL (Jacq et al. 1990). Our observations suggest a somewhat stronger enhancement of HDO compared to other high-mass star-forming sites, while it is similar to that of Orion-KL. HDO can form in the cold gas on the grains (see van Dishoeck et al. 2013 for a review on its formation) and also by gas-phase reactions at high temperatures or in shocks (e.g. Codella et al. 2010). Water deuteration is therefore strongly linked to physical conditions like the dust temperature during the freeze-out in the quiescent cloud (Cazaux et al. 2011).

Interestingly, we do not detect the deuterated form of several of the most abundant molecules in the cold envelope, for example H₂S and CH₃OH. Using the same source sizes, excitation temperatures, and line widths as for their non-deuterated species, we estimate upper limits for the column density of HDS and the deuterated forms of CH₃OH (CH₃OD and CH₂DOH), which gives upper limits of 0.001–0.002 for H2S and 0.002 for CH₃OH in the cold gas phase. The upper limits in the warm gas phase for CH₃OH are somewhat higher, around 0.007–0.02, and hence are less constraining. The deuterium fractionation for CH₃OH in the cold gas is significantly lower than typically observed towards low-mass star-forming regions (e.g. Awad et al. 2014). Interestingly, there are, however, examples of high-mass star-forming regions with considerably higher deuteration. For example Fontani et al. (2015) detected deuterated methanol with a deuteration ratio of ~0.01.

In the warm gas phase, our upper limit of deuteration is comparable to the threshold for these forms of deuterated methanol found towards the position of source N2 in Sgr(B2) with [XD]/[XH]= 0.0012 (Belloche et al. 2016) based on tentative detections of CH₂DOH lines. This corresponds to a deuterium fractionation (D/H ratio) of 0.0004 considering the different arrangements between H and D atoms. For the deuterium fractionation of CH₃OD towards SgrB2(N), an upper limit of ≤7 × 10⁻⁴ was estimated, hence the non-detection of a deuterated form of methanol towards G328.2551-0.5321 is not unexpected. This low deuterium fractionation is a clear difference to the hot corino IRAS 16293–2422 (Jørgensen et al. 2018).

We also estimate an upper limit for the deuterium fractiona-tion of H₂S, which is found to be about 0.001–0.002. The models show a low deuteration of H₂S (Hatchell et al. 1999; Awad et al. 2014), where the molecule is formed at slightly warmer temperatures on the grain surfaces. This formation route and low deuterium fractionation therefore agrees with our measurements.

Figure 15 shows a comparison of the D/H ratio of deuter-ated molecules for G328.2551-0.5321 and the hot corino, IRAS 16293–2422, and the TMC-1 low-mass core. We find that the deuterium fractionation of several molecules (HDS, HDCO, N₂D⁺, and CH₂DOH) in the cold envelope of G328.2551-0.5321 is lower than in IRAS 16293-2422, but it is similar to what is observed towards TMC-1 (Butner et al. 1995; Minowa et al. 1997; Turner 2001; Markwick et al. 2002). Other high-mass star-forming regions show similarly low deuterium fractionation (Fontani et al. 2011; Rivilla et al. 2020). High deuteration ratios can be achieved by ion-molecule reactions with the deuterated form of ${\text{H}}_3^{+}$ and ${\text{CH}}_3^{+}$, which efficiently forms in exothermic reactions at low temperatures (Roberts et al. 2003; Parise et al. 2009).

The deuterium chemistry therefore is a good indicator of the conditions at the formation time of deuterated species and reflects the temperature and the time the gas spent at a low temperature. A longer timescale for the cold collapse phase allows the reaction between HD and H₂D⁺ to lead to a higher deuterium fractionation in the cold gas. The gas temperature during the cold collapse phase also plays a key role because it sets the efficiency of the deuterium exchange reactions. Therefore, the generally low deuterium fractionation may suggest that either the temperature during the collapse is higher than that of low-mass protostars, or that the timescales for the cold collapse phase are shorter.

Fig. 14 Abundance ratios of OCS and H2S relative to SO2 for each source for the layers at Τ = 100 Κ for G328.2551-0.5321 and the sample from Herpin et al. (2009) ordered according to their assumed age sequence, G328.2551-0.5321 being the youngest and IRAS 18162-2048 the most evolved.

Fig. 15 Comparison of the deuteration in G328.2551-0.5321 (blue circles) and the dark cloud TMC1 (yellow triangles) and the hot corino IRAS 16293–2422 (pink diamonds). See the caption of Table 5 for the details about the references. The grey line represents the maximum deuteration in G328.2551-0.5321.

5.3 Isotopie ratios

Based on the LTE modelling in Sect. 3, we estimated the temperature, size of the emitting region, and molecular column density, as well as their abundance for each isotopologue. We did this individually for each molecule as the large bandwidth covered and our sensitivity allows a large enough number of transitions to be detected and fitted, as listed in Table 3. This allowed us to compute the isotopic ratios for many of these species in the cold gas phase. We discuss here, in particular, the isotopic frac-tionation for ¹²C/¹³C, ³²S/³⁴S, and ³²S/³³S. In Table 6 we list the molecules we used to estimate the isotopic fractionation, the value of the isotopic ratio, and compare it to the literature values as well as to isotopic ratios towards the Galactic centre.

In Fig. 16, we show two approaches to evaluate the isotopic ratio for the C and S atoms, and compare it to the expected values from the literature (left and right panels, respectively). We first investigated the line area ratios of the two detected transitions; (J = 2–1, J =3–2) of C¹⁸O and ¹³C¹⁸O, while for C¹⁷O and ¹³C¹⁷O, we only used the J=2–1 transition. While lines from isotopologues including ¹³C are likely to be optically thin, the transitions from C¹⁸O could be affected by a higher optical depth. Therefore, these line area ratios, measured over the entire extent of the velocity profile, place lower limits on the ¹²C/¹³C ratio, which is about 30 for the ¹⁸O, and 40 for the ¹⁷O containing isotopologues. We expect the latter to be rather close to the real ¹²C/¹³C value given the low expected optical depth of the C¹⁷O J = 2–1 line.

We also used C³⁴S and ¹³C³⁴S to further investigate the ¹²C/¹³C isotopic ratios towards G328.2551-0.5321. Using our Weeds modelling, we selected transitions with τ ~0.1–0.6 and computed line area ratios as shown in the left panel of Fig. 16, where we find a value between 19 and 26 for ¹²C/¹³C. Another approach is to use the molecular abundance estimates obtained from the Weeds modelling. We found similar temperatures and source sizes for C³⁴S and ¹³C³⁴S in our modelling, suggesting that they originate from the same gas. Hence, we used their abundance ratios to infer a ¹²C/¹³C isotopie ratio of 23, which agrees with our estimates based on line area ratios.

Because these values obtained from CO and CS isotopologues are below the expected value of 60 (Milam et al. 2005), we also checked CH₃OH and its ¹³C isotopologue. The molecular abundance ratio based on our Weeds modelling for CH₃OH and its ¹³C isotopologue gives a ¹²C/¹³C ratio of around 25, which is very close to the values obtained from the CS isotopologues. While there is scatter in the obtained ¹²C/¹³C values between 19 and 26, considering CS and CH₃OH, an uncertainty of about 40% on our Weeds column density estimates is needed to make these estimates consistent with the value of 40 obtained from the C¹⁷0/¹³C¹⁷O measurements. However, at least an order of magnitude uncertainty on the column densities obtained from Weeds would be required to make these low isotopie ratios consistent with the expected literature value of 60. Therefore, we conclude that the ¹²C/¹³C ratio is lower than the literature value of 60 towards G328.2551-0.5321. To estimate column densities from highly optically thick isotopologues, such as HCN, HNC, HCO⁺, and H₂S, we used a ¹²C/¹³C of 30, corresponding to the rough average of our estimates. In parentheses, we list the values corresponding to the canonical ¹²C/¹³C of 60 for completeness.

Similarly to the above, to estimate ³²S/³⁴S, we used optically thin transitions of ¹³CS and ¹³C³⁴S as well as SO and ³⁴SO. These molecules consistently give a value of ³²S/³⁴S = 12 (left panel of Fig. 16). Similarly, the ratio of ³⁴S/³³S could be determined by using C³⁴S and C³³S and is found to be around 3.7.

In the right panel of Fig. 16 we show the ratio of the literature values and our estimates and find that our isotopie ratios for ¹²C/¹³C, ³²S/³⁴S, and ³²S/³³S are all about a factor of 2–3 lower than expected at the galactocentric distance of G328.2551-0.5321. Isotopie ratios may exhibit variations as a function of galactocentric radius due to the variations of nucleosynthetic processing of the interstellar medium over the disk of the Milky Way (e.g. Penzias 1981a,b; Wilson & Rood 1994; Milam et al. 2005; Wouterloot et al. 2008; Giannetti et al. 2014; Jacob et al. 2020). The galactocentric variations for the ¹²C/¹³C and ³²S/³⁴S isotopie ratios are given by ¹²C/¹³C = (8 ± 2)D_GC + (8 ± 13) and ³²S/³⁴S = (3.3 ± 0.5D_GC + (3.1 ± 4.1) (Milam et al. 2005; Chin et al. 1996). Following Roman-Duval et al. (2009), we computed the galactocentric distance as $D_{\text{DG}}=\sqrt{{\left(R_0\cos \left(l\right)-d\right)}^2+R_0^2{\sin }^2\left(l\right)}$, where R₀ is the distance to the Galactic centre, l is the longitude of the source, and d is its heliocentric distance. With an estimated distance of 2.5 kpc, G328.2551-0.5321 is located at a distance of 6.4 kpc from the Galactic centre, hence its expected isotopie ratio is about 60 ± 25 for ¹²C/¹³C (Milam et al. 2005) and 25 ± 6 for ³²S/³⁴S (Chin et al. 1996), at least two times higher than suggested by several ¹³C and ³⁴S bearing molecules. Instead, the carbon and sulphur isotopie ratios resemble that of the Galactic centre. Assuming a galactocentric variation, we find that the observed values resemble the expected ratios around a galactocentric distance of 2 kpc, which is inconsistent with the location of our source, suggesting that these low values are intrinsic to our source.

While a scatter in isotopic abundances is often observed (cf. Giannetti et al. 2014), this deviation from the expected values suggests an anomalous isotopic fractionation in G328.2551-0.5321. Recent observations of hot cores reported such anomalies of unusually low ¹²C/¹³C ratio towards IRAS 20126+4104 and G31.41+0.31, where ¹²C/¹³C as low as 10 and 30 was found based on CH₃OH, and CH₃CN (and also CH₃OCHO) molecules in the warm gas, respectively (Palau et al. 2017; Beltrán et al. 2018).

Because we estimate these isotopic ratios for the cold component of the envelope, it is unlikely that internal effects due to the emerging protostar play a role. Fractionation effects caused by UV photons have been taken into account in chemical models for example in van Dishoeck & Black (1988), who showed that UV photons can significantly impact the isotopic ratios. However, the UV photons from the protostar are unlikely to be strong enough to lead to the decrease in isotopic ratios. Cosmic rays are also found to influence the isotopic fractionation chemistry. The recent study of Colzi et al. (2020) indicates that a cosmic ionisation rate higher than 1.3 × 10⁻¹⁷ s⁻¹ tends to decrease the ¹²C/¹³ C ratio. Additionally, low temperature isotopic exchange reactions have been considered in chemical networks (Roueff et al. 2015; Colzi et al. 2020; Loison et al. 2020). These reactions lead to an enhancement in the ¹³C isotopologues in molecules formed from CO and C₃ (when it does not react with atomic oxygen) and a decrement in molecules forming from atomic carbon (Loison et al. 2020). The models from Loison et al. (2020) show that isotopic fractionation for CO, CS, methanol, but also CCH, HC₃N, CN, for example, is time dependent. Additionally, their models present different isotopic fractionation for each molecules, in particular for CO, CS, and CH₃OH. The carbon isotopic fractionation appears to be stronger in CO than in CS and CH₃OH.

Isotopic fractionation of S-bearing molecules through low-temperature reactions is found to be low (Loison et al. 2019), suggesting that ³⁴S-bearing molecules are good tracers of their ³²S main isotopologue. However, our results towards G328.2551-0.5321 suggest an elevated sulphur fractionation based on CS, similarly to Tercero et al. (2010), and especially on SO molecules, which is an intriguing result in light of these models. Similarly as above, for molecules with highly optically thick lines, such as CS and H2S, we used an isotopic ratio for ³²S/³⁴S of 12, and give in parentheses the column densities using the canonical value of 25 as well.

Fig. 16 Line area ratio for the pairs of isotopologues (left panel). The ratios expected at the galactocentric distance of G328.2551-0.5321 from galactocentric variation measurements are indicated (see Table 6 for references). The error bars represent the errors on the line area given by the Gaussian fitting procedure from CLASS. Right panel: ratio of the measured isotopie ratio obtained from the LTE modelling of the column densities and the ratios expected at the galactocentric distance of G328.2551-0.5321 from galactocentric variation measurements. The error bars represent a 20% error on the column densities.

6 Molecular composition: COMs

The complete 206 GHz frequency coverage of this study offers a large range of upper-level energies for the different COMs. We detect nine COMs towards this object: methanol (CH₃OH), methyl formate (CH₃OCHO), dimethyl ether (CH₃OCH₃), formamide (HC(O)NH₂), acetaldehyde (CH₃CHO) (O-bearing COMs), methyl cyanide (CH₃CN), ethyl cyanide (C₂H₅CN), vinyl cyanide (C₂H₃CN) (N-bearing COMs), and methyl mer-captan (CH₃SH) (S-bearing COM). All these molecules except for acetaldehyde are detected in the warm component of the envelope, four (CH₃OH, CH₃OCHO, CH₃OCH₃, and HC(O)NH₂) are also found to be present in the accretion shocks based on their spectral profile (Fig. 5), and four are also identified in the cold component of the envelope (CH₃OH, CH₃CHO, CH₃OCH₃, and CH₃CN). Compared to the study of Csengeri et al. (2019) using ALMA, we detect here two additional COMs, CH₃OCH= and CH₃SH. This is because the limited spectral coverage of these ALMA observations of 7.5 GHz around 345 GHz was not suitable for the detection of these two molecules. We do not detect ethylene glycol ((CH₂OH)₂), formic acid (HCOOH), acetone (CH₃COCH₃), and ethanol (C₂H₅OH), however, which are the species with the most compact emission found in Csengeri et al. (2019), and are thus heavily diluted in the beam of the APEX telescope. The column densities, abundances relative to H₂, methanol, and dimethyl ether are reported in Table 4 and Fig. 17. We also list upper limits for several COMs in the cold gas component, where we assumed a size of 26″ and a temperature of 16 K, which corresponds to an average of the physical properties of COMs detected in the cold gas component.

6.1 Accretion shocks in G328.2551-0.5321

The accretion shocks reported in Csengeri et al. (2018) using methanol, and then typically O-bearing COMs (Csengeri et al. 2019) show a velocity offset of ~± 4.2 kms⁻¹ relative to the source v_lsr (Fig. 5). Interestingly, in high-excitation transitions of several COMs, such as CH₃OH, CH₃OCHO, CH₃OCH₃, HC(O)NH₂, and HDO, we also see these velocity components even with a single-dish telescope. In addition, we detect several rotational transitions of methanol from within its torsion-ally excited state, v_t = 1, which only trace the accretion shocks (Csengeri et al. 2018).

Because we have more transitions in our spectral survey than Csengeri et al. (2019), we attempted to place more stringent constraints on the gas temperature based on our LTE modelling. Our models provide parameters for the two velocity components that have been spatially resolved in Csengeri et al. (2019). These velocity components are blended with the envelope component in several of these cases, hence we find a larger dispersion for the peak velocity of the accretion shocks. We estimate a gas temperature between 170 and 190K for the accretion shocks (Fig. B.8) based on all the molecules that show these velocity components. This is consistent with the range of values around 180 K found by Csengeri et al. (2018) based on CH₃OH alone. The column densities and hence the abundances are about an order of magnitude lower compared to the spatially resolved measurements based on ALMA, however, which is expected because our observations do not spatially resolve these regions.

6.2 COMs in the eold component of the envelope

We were able to detect four COMs, CH₃OH, CH₃OCH₃, CH₃CHO, and CH₃CN, in the cold component of the envelope with a gas excitation temperature between 20 and 34 K. Towards low-mass cores, COMs were found at low temperatures (Öberg et al. 2010; Bacmann et al. 2012; Jaber et al. 2014; Vastel et al. 2014), but in the vicinity of hot cores, most studies reveal COMs in the warm component of the envelope (cf. Bisschop et al. 2007; Widicus Weaver et al. 2017). The derived column densities in the cold phase (2 × 10¹⁵ cm⁻² for methanol and 5 × 10¹³–10¹⁵ cm⁻² for the other molecules) are about two orders of magnitude lower than in the warm phase, while the abundances relative to H₂ (1.9 × 10⁻⁸ for methanol and 3.4 × 10⁻¹⁰–1.5 × 10⁻⁹ for the other molecules) are about one order of magnitude lower in the cold component of the envelope (Fig. 17).

The formation of COMs on the surface of dust grains and ices at low temperatures was shown to be possible by Garrod et al. (2008). From the grain mantles, they are then liberated into the gas phase through non-thermal processes caused either by shocks, supersonic turbulence (Requena-Torres et al. 2006), energy release from exothermic reactions on the grains (Garrod et al. 2007), UV radiation originating from the protostar propagating through the outflow cavity (Öberg et al. 2010; Drozdovskaya et al. 2015), or by cosmic ray radiation or secondary UV photons (produced by the electron recombination after the impact of cosmic rays on H₂ molecules) (Bacmann et al. 2012). The high abundances of sulphur-bearing species suggest the presence of shocks (Sect. 5.1.3), and the presence of COMs in the cold component of the envelope could be also explained by the sputtering of grain mantles due to such shocks. No strong UV radiation can be expected from the protostar itself because it does not show a strong radio continuum emission, and energetic photons from the protostar could escape through the cavity formed by a collimated outflow. Chemical desorption or low-velocity shocks could also give rise to non-thermal desorption, leading to the observed cold-phase COMs. It then remains to understand to which extent this composition of COMs reflects that of the grain surfaces. The cold component of the envelope exhibits a diverse molecular composition of simple molecules including a large number of molecular ions. This suggests rather active chemical processes in the cold gas phase, which is expected to represent the pristine, pre-collapse gas not affected by the heating of the central object. Low-temperature gas-phase chemistry could therefore have an impact on the molecular composition of COMs.

Fig. 17 COM column densities (top left panel) and abundances relative to H2 (top right panel), methanol (bottom left panel), and dimethyl ether (bottom right panel). Only the main isotopologues are represented. The cold and warm components of the envelope are described by the blue and red histograms. The two shocks are represented in yellow and green. The upper limits are indicated by the black arrows.

6.3 COMs in the warm envelope

After methanol, the second most abundant COM in the warm component of the envelope is dimethyl ether. Both molecules are detected in the cold and warm components of the envelope as well as in the accretion shocks. Because they are chemically linked (Garrod et al. 2008), they are expected to be found in the same gas components. We also detect methyl formate at an abundance relative to H₂ of 1.5 × 10⁻⁸ and temperature of 86 K. We detect it in the accretion shocks and in the warm component of the envelope, but the frequency range of our survey is not suitable to detect it in the cold component of the envelope. Similarly, CH3CHO is chemically linked to these species, which is present in the warm component seen by ALMA (Csengeri et al. 2019), but our sensitivity is not good enough to detect it because it originates form a very compact region in the vicinity of the protostar. Methyl formate and dimethyl ether are only an order of magnitude less abundant than methanol.

In G328.2551-0.5321, COMs with a CN functional group (vinyl and ethyl cyanides) are detected in the warm gas (90–100 K) with high abundances of 2.6 × 10⁻⁸ for C₂H₅CN and 3.5 × 10⁻⁹ for C₂H₃CN. They exhibit a slightly higher temperature than O-bearing COMs, but this difference would be even more pronounced if we had used a two-component model for ethyl cyanide (see Sect. 3.4.2), suggesting that the gas temperatures in the immediate vicinity of the protostar are still moderate, again pointing to a very compact radiatively heated inner region. Gas temperatures of these cyanides are found to be typically higher than for O-bearing COMs (Widicus Weaver et al. 2017), and towards hot cores, they often trace the hot gas component above 150 K (Belloche et al. 2016). We find that the case of CH₃CN is somewhat particular. We detect several high-energy transitions from the K = 7–9 lines. Our model with a temperature of 200 K, a size of 1″, and a high column density of 1.5 × 10¹⁶ cm⁻² provides a relatively good fit to the K = 0–6 lines, but it underestimates the K = 7–9 components, suggesting that a hot compact component traced by CH3CN could be present. A hot envelope component with temperatures up to 280K has already been observed with single-dish telescopes towards hot cores using the CH3CN line (e.g. Widicus Weaver et al. 2017; Bisschop et al. 2007). Similar energy levels are available in our frequency coverage for methanol as well, but they do not show this hot component.

We also detect the sulphur-bearing COM methanethiol (CH3SH) towards the warm component of the envelope at a temperature of 69 K and with an abundance relative to H2 of 3.4 × 10⁻⁹. CH₃SH represents only a small fraction of the total sulphur reservoir.

6.4 Origin of COMs in the envelope of G328.2551-0.5321: Direct sublimation of grain-surface products?

We have shown that based on their typical temperature, line profile, and size of the emitting region, three different physical components can be distinguished in the envelope of G328.2551-0.5321, corresponding to the cold and warm gas components and the accretion shocks. The overall molecular composition of the cold and warm envelope components is very diverse, likely pointing to a different chemical history of the gas. Based on COMs detected both in the cold and warm gas phases, we investigate here their potential formation based on their abundance ratios in the cold and warm gas listed in Table 7. Because only a few COMs are detected both in the cold and in the warm gas, we are limited to perform only a qualitative comparison to models, for which we chose the gas-grain chemical models by Garrod & Herbst (2006); Garrod et al. (2008, 2017).

In the cold gas, we detect methanol, dimethyl ether, and acetaldehyde at temperatures between 20 and 27 K with molecular abundances with respect to H₂ of 1.9 × 10⁻⁸, 1.5 × 10⁻⁹, and 3.4 × 10⁻¹⁰, respectively (Sect. 6.2). Pure gas-phase chemistry in the cold gas phase is not sufficient to explain such high abundances for COMs, although its role has been emphasized in recent, revised chemical models of low-mass prestellar cores (Vasyunin et al. 2017). Garrod & Herbst (2006) demonstrated that methanol, methyl formate, and dimethyl ether may form efficiently on the grain surfaces, which then sublimate to the gas phase when the temperature increases. Therefore, the detection of COMs in the cold gas towards G328.2551-0.5321 suggests grain surface formation and subsequent desorption to the gas phase. In the cold envelope, this invokes non-thermal desorp-tion processes that need to be efficient to lead to the observed molecular abundances of these COMs.

The warm gas typically exhibits molecular abundances in methanol, dimethyl ether, and methyl formate that are an order of magnitude higher than the cold gas, demonstrating that the inner, warmer regions become richer in saturated oxygen-bearing COMs. The difference of abundances relative to H₂ between the cold and warm components of the envelope could be explained by a considerably higher efficiency of the thermal desorption in the warm gas phase. Observations suggest a tight chemical connection between methanol, methyl formate, and dimethyl ether (e.g. Garrod & Herbst 2006; Brouillet et al. 2013; Jaber et al. 2014), and indeed we detect these molecules in high abundance in the warm gas and in the accretion shocks. Their abundances in the warm gas are qualitatively consistent with the models of Garrod & Herbst (2006) and Garrod et al. (2008, 2017), who predicted similar abundances around 10⁻⁸−10⁻⁷ for methyl formate and dimethyl ether, while the peak abundance of methyl formate decreases in their models as a function of warm-up time. Our observations suggest that the abundances of methanol are an order of magnitude lower than in these models, therefore the abundance ratios of methyl formate and dimethyl ether to methanol are about an order of magnitude lower only. This is consistent with observational results of other hot cores (e.g. Taquet et al. 2015), but it suggests some discrepancies to these chemical models.

The abundance ratio of CH₃OCH₃ to CH₃OH in the cold and warm components of the envelope are different by a factor of only ~4 (Table 7). These ratios are representative of the difference of chemistry and complexity between the cold and warm components of the envelope. The higher CH₃OCH₃/CH₃OH ratio in the warm component of the envelope may suggest a higher degree of molecular complexity.

Formamide is a potential precursor to heavier COMs, such as amino acids (cf. López-Sepulcre et al. 2015). While obser-vationally it is detected in various environments, its formation conditions remain debated. On the grains, the successive hydro-genation of HNCO is thought to be efficient (Charnley 1997; Mendoza et al. 2014; López-Sepulcre et al. 2015), but the link between HNCO and HC(O)NH2 is still debated in experiments (e.g. Noble et al. 2015; Haupa et al. 2019). During the cold phase on the grains, further hydrogenation can lead to more complex molecules (Raunier et al. 2004). At ~100K, formamide desorbs and is present in the gas phase. Towards G328.2551-0.5321, the ratio of HC(O)NH2 to HNCO is 0.04 in the warm gas phase at a temperature of ~80K, slightly below its desorption temperature. This ratio is somewhat lower than what was found towards high-mass protostars (López-Sepulcre et al. 2015), and it is considerably lower than the extreme star-forming region, SgrB2(N2), for which Belloche et al. (2017) reported a ratio of HC(O)NH₂/HNCO of 1.8. When comparing the molecular abundance of HC(O)NH₂ with respect to H2 in the sample of Bisschop et al. (2007), we see that the abundance of G328.2551-0.5321 is a factor of a few lower than that of these hot cores. Considering the lower temperatures in the warm component of the envelope, this could suggest that the thermal desorption has not been complete and only part of the formamide has been liberated from the grains, if it formed on the grain surfaces.

Another formation route of formamide is in the gas phase from H₂CO and NH₂ (Skouteris et al. 2017). This reaction has an energy barrier of 26.9 K (Barone et al. 2015), which makes it efficient even at the low temperatures of the cold gas towards G328.2551-0.5321, where the bulk of the gas is at 20 K. The models of Vasyunin et al. (2017) show that including this reaction can lead to abundances relative to H2 of 5 × 10⁻¹² in the cold component of the envelope. However, Belloche et al. (2017) showed that when this reaction is included in Sgr B2(N2), it leads to an overestimation of the abundance of formamide in the warm component of the envelope, favouring a grain surface production. Such an abundance is in any case below our detection limit and is therefore consistent with the non-detection of formamide in the cold gas phase.

The chemical network leading to the formation of cyanides has been modelled for example by Belloche et al. (2009) and Garrod et al. (2017), who suggested that C₂H₅CN can be formed on the grains by the hydrogenation of C₂H₃CN formed from HC₃N during the collapse and the warm-up phase. C₂H₅CN then sublimates to the gas phase at a temperature of ~ 100 K (Garrod et al. 2017). Subsequent gas-phase reactions lead to the destruction of most of the C₂H₅CN formed on the grains, hence the ratio of ethyl to vinyl cyanide has been discussed as a potential chemical clock (cf. Allen et al. 2018, but see also Caselli et al. 1993) reflecting the minimum age of the protostar. Towards G328.2551-0.5321, we find here a ratio of ethyl to vinyl cyanide of ~6, which is broadly consistent with the ratio of 2.4–3.4 found by Csengeri et al. (2019), who based their abundance estimates on adopting the kinetic temperature of methanol for all COMs. Our results confirm that saturated nitriles are more abundant than unsaturated nitriles in this source, which is expected as they are chemically more stable. On the other hand, comparing the ratio of ethyl to vinyl cyanide to the values of other hot cores from Fontani et al. (2007), we find here a somewhat higher value, implying a larger amount of C₂H₅CN, suggesting that the source could be considered as younger than the hot cores in their sample.

We do not detect any of these complex nitriles in the cold envelope, although Garrod et al. (2017) showed that vinyl cyanide can be formed efficiently through gas-phase reactions between CN and C₂H₄ up to 100 K, leading to gas-phase abundances of ~10⁻¹¹, which is below our detection limit. Therefore, our observations within the limitation of the assumptions for the upper limits given in Table 4 are not sensitive enough to distinguish between the models of Garrod et al. (2017).

We saw in Sect. 5.1.1 that the ratio of methanol and methanethiol is lower than what is observed in SgrB2(N2), OrionKL, or G327.3 or by chemical modelling, suggesting an enhancement of the S-bearing molecules in this object. CH3OH and CH₃SH are thought to be formed by successive hydrogena-tion of H₂CS and H₂CO on the ices and then released to the gas phase (Gibb et al. 2000; Majumdar et al. 2016). In the warm component of the envelope, the ratio CH₃SH/H₂CS(o+p) is equal to 0.2, while CH₃OH(A+E)/H₂CO(o+p) is close to two orders of magnitude higher, which is in agreement with the chemical modelling from Müller et al. (2016). This can be explained by subsequent reactions in parallel specific to the sulphur and oxygen chemistry. For instance CH₃SH can be formed through reactions between the two radicals CH₃ and SH, and H₂CS through reactions between CH₂ and S (Müller et al. 2016). CH₃SH/H₂CS is lower by an order of magnitude than what is observed in IRAS 16293-2422 (~3) and in comet 67P (~3) (Drozdovskaya et al. 2018).

In the cold component of the envelope, the only molecule that is well reproduced by the models from Garrod et al. (2008) is acetaldehyde, for which we estimated an abundance of 3.4 × 10⁻¹⁰. Acetaldehyde forms on the grains like other complex molecules, but can desorb at low temperatures.

To summarise, the cold component of the envelope of G328.2551-0.5321 is richer in small molecules than the warm gas phase. In contrast, the warm component of the envelope shows a richness in COMs that suggests an increasing chemical complexification at higher temperatures, although several heavier COMs that are expected to form on the grains remain undetected. This suggests chemical differences in the cold and warm components of the envelope, favouring a scenario in which the molecular composition of the gas in the cold component of the envelope is chemically evolved, but it is poor in COMs due to its low temperature. The warm component of the envelope is expected to inherit the molecular composition of the grain mantles, which have reached the sublimation temperature only recently as the protostar is young. At the estimated age of the protostar, no physical and chemical equilibrium could be reached.

7 Comparison of G328.2551-0.5321 to hot corinos and hot cores

Csengeri et al. (2019) proposed G328.2551-0.5321 to be a precursor of hot molecular cores because its radiatively heated region is more compact and is at a moderate temperature compared to hot cores. Its richness in COMs has been first attributed to the impact of accretion shocks, which has been shown to lead to an increased abundance of O-bearing COMs. Here we can disentangle the contribution of the accretion shocks from the molecular composition of the warm inner region, which in addition to being among the most compact of hot cores (i.e. Allen et al. 2018; Ginsburg et al. 2018), also shows a typically colder temperature of about 100 K compared to the classical hot cores, which reach higher temperatures (Widicus Weaver et al. 2017), typically >150 K within a size scale above ~0.01 pc (e.g. Sánchez-Monge et al. 2014). Here we aim to place the molecular composition of the cold and warm gas phase of G328.2551-0.5321 into context by comparing it to the archetypical hot corino IRAS 16293–2422 and hot cores from Bisschop et al. (2007).

To do this, we compared our detections with the findings of the ASAI (Lefloch et al. 2018) and TIMASS (Caux et al. 2011) surveys. In Sect. 3 we have classified the 2935 identified lines into the categories: O-, S-, and N-bearing and carbon chain molecules. The larger number of O-bearing lines compared to that of carbon chains confirms that the molecular content of G328.2551-0.5321 is more similar to hot cores and hot corinos than to WCCC objects (Lefloch et al. 2018). The ratio of O-bearing and N-bearing lines is similar to that of IRAS 16293–2422 (Caux et al. 2011). The proportion between the different types of molecules also highlights a large number of lines from S-bearing molecules that is more similar to shocked regions, such as L1157-B1 (Lefloch et al. 2018).

We show a quantitative comparison of the molecular abundances of G328.2551-0.5321 and IRAS 16293–2422 in Fig. 18, based on single-dish observations of the O-bearing COMs and the light molecules (Schöier et al. 2002; Maret et al. 2005; Jaber et al. 2014; Taquet et al. 2014; Majumdar et al. 2016), where the cold and warm components of the envelope have been distinguished. In the warm envelope, we were able to compare the O-bearing molecules as well as HDO, HC₃N, and S-bearing molecules. For the comparison of the cold component of the envelope, we could select S-, N-, O-molecules, carbon chains, but also the four COMs detected in the cold gas, which are detected in both objects. As Fig. 18 illustrates, the cold gas shows higher molecular abundances for all species in G328.2551-0.5321, and we find abundances typically higher by an order of magnitude than IRAS 16293-2422. This enhancement is not specific to one type of molecules, but to the ensemble of the cold component of the envelope. However, molecules in the warm gas exhibit similar abundances in both objects, with values between 10⁻⁹ and 10⁻⁷.

Finally, we also compare our object to hot cores from Bisschop et al. (2007), where molecules originating from the cold and warm components of the envelope could be distinguished. We find that the warm gas towards G328.2551-0.5321 has an abundance of COMs lower by an order of magnitude relative to H₂ compared to their sample of hot cores.

Overall, we detect in G328.2551-0.5321 similar families of molecules compared to hot corinos and hot cores, although both hot cores and hot corinos themselves show some diversity in their molecular composition. Quantitatively, however, using single-dish molecular abundance measurements with distinction between the warm and the cold components of the envelope for these samples, we find that the molecular abundances of COMs with respect to H₂ resemble that of the hot corino associated with IRAS 16293-2422 rather than hot cores, which typically exhibit higher molecular abundances of COMs in their warm gas phase than G328.2551-0.5321.

Fig. 18 Molecular abundances relative to H2 towards IRAS 16293-2422 (Schöier et al. 2002; Maret et al. 2005; Jaber et al. 2014; Taquet et al. 2014; Majumdar et al. 2016) and hot cores from Bisschop et al. (2007) vs. that of G328.2551-0.5321. The diamond represents the COMs in the hot cores. The small symbols are used to show each hot core. The larger symbol shows the averaged value. The blue data points represent the cold component of the envelope, and red represents the warm component of the envelope. The grey shaded region corresponds to a factor of 5 scatter from the one-to-one relation.

8 Summary and conclusions

We performed a spectral survey with the APEX telecope between 159 GHz and 374 GHz towards the high-mass protostar, G328.2551-0.5321. We detected 39 species (61 including isotopologues), and report tentative detections for the NS⁺, PN,

DNC, and HCNH⁺ species. The detected species include carbon chains, O-, N-, S-bearing, deuterated molecules, and COMs. Using an iterative process relying on rotational diagrams and LTE modelling with Weeds, we estimated the physical conditions (e.g. column density, excitation temperature, and size of the emitting region), and based on a simple source model, we inferred the molecular abundances of these species. Our main results are listed below

• We identify 2935 lines from O-, S-, and N-bearing and carbon chain molecules. The proportion between the different types of molecules suggests that the molecular content of G328.2551-0.5321 is more similar to hot cores and hot corinos than to WCCC objects, but we find high abundances relative to H₂ from sulphur-bearing species.

• While most lines show a Gaussian profile, we also find kinematic signatures of outflow and infall motions. Several transitions require a two-component Gaussian fitting, corresponding to a narrow and broad velocity component characterised by Δv < 6kms⁻¹, and Δv > 6 km s⁻¹, where the position of the line peak coincides with the source v_lsr. The quiescent gas component is found to show T_ex < 50 Κ (for Δv < 6kms⁻¹), while the broader velocity component exhibits T_ex > 50Κ (for Δv > 6kms⁻¹). These two components correspond to a warm and cold gas phase in the envelope.

• Based on their line profiles, we also identify emission from the accretion shocks that show a velocity offset of ~ ± 4 km s⁻¹ with respect to the source v_lsr. Rotational transitions from the first torsionally excited state of CH₃OH exclusively pinpoint the accretion shocks, while we find that these velocity components are also present in the rotational transitions from the torsional ground state of CH₃OH, and rotational transitions from CH₃OCH₃, CH₃OCHO, HC(0)NH₂, and HDO. The accretion shocks exhibit the narrowest line profiles with Δv ~3 km s⁻¹, and are associated with gas at T_ex = 170–190 K.

• The molecular composition of the cold gas is characterised by small molecules and four COMs, CH₃OH, CH₃OCH₃, CH₃CHO, and CH₃CN, while among the detected molecules, the COMs (CH₃OH, CH₃OCH₃, CH_зOCHO, CH₃CHO, HC(0)NH₂, CH₃CN, C₂H₅CN, C₂H₃CN, CH₃SH) dominate the chemical composition of the warm envelope.

• The cold component of the envelope is particularly rich in sulphur-bearing molecules with abundances as high as 3.9 × 10⁻⁸ for H₂S, which is even higher than the CH₃OH abundance of 1.9 × 10⁻⁸. Sulphur-bearing ions such as HCS⁺ and SO⁺ have abundances of 1.5–1.9 × 10⁻¹⁰, respectively. Overall, we find that the abundances of sulphur-bearing species are higher than typically observed towards other star-forming regions, suggesting that sulphur is either not as strongly depleted as towards other objects with a depletion greater than 1%, or that efficient mechanisms liberate sulphur-bearing species into the gas phase.

• We detect DCN, DCO⁺ HDCS, and HDCO towards the cold envelope component, while HDO is found in the warm component of the envelope as well as towards the accretion shocks. Similarly to other hot cores, in these molecules, G328.2551-0.5321 exhibits a deuterium fractionation of 0.2%–3%, which is lower than for hot corinos. However, the deuteration values in G328.2551-0.5321 are found to be close to what is observed in dark clouds, suggesting that the prestellar stage might be shorter or warmer than the low-mass star formation scenario.

• We investigate the isotopic ratios within the cold gas and find that G328.2551-0.5321 exhibits an unusually low isotopic ratio of ¹²C/¹³C of ~30 based on column density estimates of C³⁴S and ¹³C³⁴S. Based on ¹³CS and ¹³C³⁴S, we also find a lower than expected value for the ³²S/³⁴S isotopic ratio of 12, which is similarly about half of its expected value. The ³⁴S/³³S isotopic ratio is about half of its expected value.

• Consistent with the conclusion of Csengeri et al. (2019), we identify here emission from the warm component of the envelope, which is found to be compact (typically <2′′) and exhibits a moderate temperature of ~63–166 K that is lower than that of the accretion shocks. Despite the presence of warm gas, the global molecular emission is dominated by the cold component of the envelope, with a temperature around 20K.

• Four COMs are identified in the cold gas phase of the envelope (CH₃OH, CH₃OCH₃, CH₃CHO, CH₃CN) showing molecular abundances similar to that of chemical models of grain surface production and non-thermal desorption. – The richness in COMs within the warm gas suggests an increase in chemical complexity. The differences in molecular abundances between the cold and warm gas favour different chemical processes in action. The observed temperatures in the warm gas are found to be close to the thermal desorption temperatures of several COMs, suggesting that the molecular composition of the warm gas may reflect that of the grains.

• Comparing the molecular richness to the archetypical hot corino, IRAS 16293-2422, we find that it exhibits a similar richness in the number of molecules detected in its cold gas. A quantitative comparison between the cold and warm gas phases of the archetypical hot corino IRAS 16293-2422 and G328.2551-0.5321 shows that the cold component of the envelope has molecular abundances that are an order of magnitude higher in G328.2551-0.5321, while the warm gas exhibits similar molecular abundances in both types of objects. While the molecular emission of the warm gas is similar to that of hot cores, the molecular abundances of COMs relative to H₂ towards typical hot cores are found to be an order of magnitude higher than that of G328.2551-0.5321 and IRAS 16293-2422. The moderate temperature of the warm gas phase together with abundances of COMs that are an order of magnitude lower than expected may suggest that thermal desorption has not yet liberated all grain surface products to the gas phase.

This unbiased spectral survey allowed us to investigate the molecular composition of G328.2551-0.5321 in great detail. The number of detected molecules suggests an emerging molecular complexity in this envelope, in which a cold and warm gas phase can be distinguished. Altogether, this points to a picture in which the overall molecular diversity of the source is similar to that of hot cores and hot corinos, but the gas in the vicinity of the high-mass protostellar embryo resembles that of a hot corino in the warm gas phase, rather than a hot core. A statistical study of high-mass protostars at a similarly early evolutionary stage should reveal whether this deeply embedded hot corino stage is a common phenomenon in the early evolutionary phase of high-mass star formation.

Acknowledgements

We thank the referee for comments that helped improve the manuscript. This publication is based on data acquired with the Atacama Pathfinder Experiment (APEX). APEX is a collaboration between the Max-Planck-Institut für Radioastronomie, the European Southern Observatory, and the Onsala Space Observatory. L.B. thanks the following institutions for support: the International Max-Planck-Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne, the Bonn-Cologne Graduate School (BCGS) and the Collaborative Research Centre 956 (sub-project B6) funded by the Deutsche Forschungsgemeinschaft (DFG). T. Cs. has received financial support from the French State in the framework of the IdEx Université de Bordeaux Investments for the future Program.

Appendix A Simple molecules

A.1 Carbon chain molecules

Carbon chains are unsaturated molecules, for example C_nH, c-C₃H₂, CH₃CCH, and HC_2n+1N. Because we do not cover the 3 mm atmospheric window, we do not have access to the heavier carbon chain molecules in the cold phase, such as HC_2n+1N, where n > 1, because their spectral lines in the range of our survey are too weak to be detected at low temperature. While prevalent in quiescent molecular clouds, they have also been studied in detail towards nearby protostellar objects (Turner 1991; Beuther et al. 2008; Sakai et al. 2008, 2010), which led to the definition of protostellar envelopes with warm carbon chain chemistry (WCCC), which are characterised by being particularly rich in unsaturated hydrocarbons (e.g. Sakai et al. 2008). This molecular composition is different from hot cores and hot corinos (e.g. Ceccarelli et al. 2000, 2007), which are instead dominated by saturated hydrocarbons and exhibit a rich emission from COMs. However, there is at least one exception, the low-mass protostar L483, which exhibits the characteristics of both chemical types (Oya et al. 2017). Because carbon chain molecules are predicted to form in larger amounts in gas heated to lukewarm temperatures of around 30 K (Sakai & Yamamoto 2013), the difference in chemistry for hot corinos and WCCC sources has been explained by a long versus short starless-core phase, implying different physical conditions in the pre-collapse gas.

Towards G328.2551-0.5321, we detect the carbon chain species CCH, c-C₃H₂, and CH₃CCH in the cold component of the envelope, typically exhibiting low excitation temperatures of 13–35 K. We also identify emission from HC₃N in the cold component of the envelope at a similar temperature. However, this molecule is also detected in the warm component of the envelope with a similar abundance. We estimate an upper limit for c-C₃H₂ and CH₃CCH in the warm envelope that is at least an order of magnitude lower than what is observed in the cold phase. Heavier unsaturated species such as C₄H, l-C₃H, and C₂S, C₂O seen towards the archetypical WCCC source, L1527, (Sakai et al. 2008) are not detected. The most promising lines of these molecules have their low upper-level energy transitions at lower frequencies than our survey, but their non-detection in our survey suggests that they are not present in the warm component of the envelope, with an upper limit on their abundances relative to H₂ of ~10⁻¹⁰ (10⁻⁹ for C₄H). These results are consistent with previous observations suggesting that the lightest carbon chain species are frequently detected towards high-mass star-forming regions (Sakai et al. 2008, 2010; Sakai & Yamamoto 2013).

Recent observations suggest that some carbon chain species, such as CCH and HC₃N, could be as abundant as CH₃OH (Taniguchi et al. 2018, 2021). The molecular abundances in the cold material resemble what is observed towards other hot cores and hot corinos (Aikawa et al. 2008; Sakai et al. 2010). The main formation path for unsaturated carbon chain molecules is mostly through reactions of CH4 and C⁺ (Aikawa et al. 2008; Sakai & Yamamoto 2013). The desorption of CH4 from the grains is expected at a temperature of ~25 K, which then reacts with C⁺ in the gas phase. This formation route is possible only for low temperature. When the temperature increases above ≥ 30 K, other radicals desorb from the grain surfaces (Garrod et al. 2008), and C⁺ reacts preferentially with them to form COMs. Therefore, WCCC develops preferentially when the temperature is just enough to desorb CH4 from the grains at a temperature of ~25K because this molecule is not strongly bound to the ices. In G328.2551-0.5321 the carbon chains in the cold component of the envelope might be formed through these reactions, which would explain their observed temperatures. HC³N is found in similar abundances in the cold and warm components of the envelope, but our upper limits show that other carbon chains have an abundance at least an order of magnitude lower in the warm component of the envelope, which suggests different chemistry in the cold and warm components of the envelope. Therefore, HC3N is more likely to form at high temperature through the reaction between H₂ and C₃N (Hassel et al. 2011), which could explain the presence of HC₃N in the hot gas phase.

A.2 N-bearing molecules

We detect several N-bearing molecules towards G328.2551-0.5321. Most of them are composed of H, N, and C atoms, CN, HNC, HCN, N₂H⁺, and CH₂NH, and two with an oxygen atom, HNCO and NO. All these molecules are important precursors of a more complex chemistry.

Methylenimine, CH₂NH (together with CH₃NH₂, CH₃CN, and HC(O)NH₂) is considered as a prebiotic molecule and potential precursor of glycine (e.g. Holtom et al. 2005). We detect it in the cold component of the envelope at a temperature of 43 K with a column density of 3.0 × 10¹⁴cm⁻² and an abundance relative to H₂ of 4.5 × 10⁻¹⁰. This temperature is similar to that of other COMs in the cold gas phase (see Sect. 6).

HNCO is thought to be a potential precursor of formamide (Mendoza et al. 2014; López-Sepulcre et al. 2015), but this link is debated based on different laboratory experiments (e.g. Noble et al. 2015; Haupa et al. 2019). We detect HNCO both in the cold and warm gas phase, but with a difference of two orders of magnitude in column density. The abundance relative to H₂ of 3.7 × 10⁻¹⁰ in the cold component of the envelope is similar to what is found in the envelope of low-mass protostars (López-Sepulcre et al. 2015), and the abundance relative to H₂ of 1.7 × 10⁻⁸ in the warm component of the envelope is similar to what is observed towards hot cores (Bisschop et al. 2007).

We also detect NO in the cold component of the envelope with a column density of 6.0 × 10¹⁵ cm⁻² and an abundance relative to H₂ of 2.3 × 10⁻⁸ at a temperature of 14 K. This value is similar to what is found in quiescent gas towards low-mass protostars (Akyilmaz et al. 2007; Codella et al. 2018) and dark clouds like TMC-1 (Gerin et al. 1993).

Appendix B Rotational diagrams

B.1 Carbon chains

Fig. B.1 Rotational diagrams of the carbon chains CH3CCH, CCH, C-C3H2, and HC3N. HC3N is described by three rotational diagrams showing the three temperatures traced by HC3N. The first rotational diagram is for the narrow component and low temperature. The two following rotational diagram are the two rotational diagrams for the two temperature components covered with the larger line width. Finally, we show the vibrationally excited state of HC3N, which is a different entry and represents the column density of the molecule present in this vibrationally excited state.

B.2 O-bearing molecules

Fig. B.2 Rotational diagrams of the O-bearing molecules HCO+ H213CO, and H2CCO.

B.3 Ν-bearing molecules

Fig. B.3 Rotational diagrams of the N-bearing molecules HCN, HN13C, CH2NH, NO, HNCO, and N2H+. The two rotational diagrams of HNCO represent the narrow and the broader components of the Gaussian fit.

B.4 S-bearing molecules

Fig. B.4 Rotational diagrams of the S-bearing molecules CS, H2 34S, ortho-H2S, OCS, SO2, and SO. The emission from the cold component of the envelope for H2S is described by the rotational diagram of o-H2 34S. The two rotational diagrams of OCS, SO2, and SO represent the narrow and the broader components of the Gaussian fit.

Fig. B.5 Rotational diagrams of the S-bearing molecules NS, HCS+, and SO+.

B.5 Deuterated molecules

Fig. B.6 Rotational diagrams of the deuterated molecules HDCO, DCO+ DCN, and HDO.

B.6 Others

Fig. B.7 Rotational diagrams of SiO.

B.7 COMs

Fig. B.8 Rotational diagrams of the COMs CH3OH, CH3OCH3, CH3CHO, HC(0)NH2, and CH3OCHO. The three rotational diagrams of methanol and dimethyl ether represent the rotational diagram of the three components (cold and warm components of the envelope and shocks) traced by these molecules.

Fig. B.9 Rotational diagrams of the COMs CH3CHO, CH3CN, C2H5CN, C2H3CN, and CH3SH.

Appendix C Spectral survey of G328.2551-0.5321

Fig. C.1 Spectrum towards G328.2551-0.5321. The grey filled histogram shows the spectra and the coloured labels the detected transitions of all the light molecules and some transitions of the COMs. The green line represents the LTE model obtained with Weeds. The upper panel for each frequency range shows the labels for the small molecules, and the lower panel corresponds to a zoom in temperature for better visibility of the COMs that are labeled.

References

All Tables

All Figures

	Fig. 1 Continuous frequency coverage of the spectral survey. The black curve shows the atmospheric transmission as computed by the ATM model written by Juan Pardo (Pardo et al. 2002) with a water vapour of 1.2 mm.
In the text

	Fig. 2 Observed spectral windows for the 1 mm band with the PI230 receiver at APEX. The solid lines represent the inner backend units, and the dashed lines show the outer backend units.
In the text

Fig. 3 Spectrum towards G328.2551-0.5321 between 252 GHz and 258 GHz. The filled grey histogram shows the spectra, and the coloured labels show the detected transitions of all the light molecules and some transitions of the COMs. The green line represents the LTE model obtained with Weeds. The upper panel for each frequency range shows the labels for the small molecules, and the lower panel corresponds to a zoom in temperature for better visibility of the COMs which are label.

In the text

	Fig. 4 Spectra, in black, of DNC(3–2), NS+ (4–3), PN(5–4), HCNH+ (3–2) shown as filled grey histograms. The dashed red line shows the systemic velocity of the source at –43.13 kms−1 (Csengeri et al. 2018). The blue line represents the best Gaussian fit for each line. Eup represents the upper-level energy of the transition in temperature units.
In the text

Fig. 5 Examples of the different line profiles observed in the spectral survey. From the top to the bottom, we show in the left column a narrow line profile for HC18O+ (3–2) and HC3N (25–24) with a two-component Gaussian fit and HC3N (33–32) with a single broad component. In the right column,HCO+ (3–2) has a line profile with redshifted self-absorption, SiO (4–3) has a broad component coming from the outflow, and the profile of CH3OH (11−5–12−4) shows two velocity components tracing the accretion shocks. All the lines are plotted in a window of 70 km s−1 centred on the rest velocity. The dotted red line shows the systemic velocity of the source at –43.13 km s−1.

In the text

Fig. 6 Hyperfine structure fitting of the transition CCH (2–1) in the 2 mm band. The result of the fitting procedure is indicated in blue, and the systemic velocity corresponds to the frequency of the most intense component hyperfine transition. The quantum numbers follow the notation from (Gottlieb et al. 1983), with Ν corresponding to the rotational levels, J to the spin doubling, and F to the hyperfine levels.

In the text

Fig. 7 Rotational diagram of HC3N for the broad velocity component towards G328.2551-0.5321. The rotational diagram when we assume that HC3N lines are optically thin is shown in blue. We represent the rotational diagram corrected for the opacity of each line in red. Rotational temperatures and column densities are indicated with the same colour code. The beam dilution was calculated for each transition for a source size of 4″. The error bars represent a lσ- statistical noise plus a 20% calibration error. To calculate the uncertainties on the rotational temperature and the column density, we used the Monte Carlo method and assumed a uniform distribution of the error for each data point.

In the text

Fig. 8 Excitation temperature vs. mean line width. A statistical study of K-Means (MacQueen 1967) was performed to determine the limit between the different components of the envelope. The blue data points represent the cold component, and the warm component of the envelope is represented by the red data points. Single transitions fitted with multiple velocity components are considered as individual measurements in this figure.

In the text

Fig. 9 Excitation temperature vs. the size of the emission region for each molecule whose emission arises from the envelope. The black curve represents the power-law fit of the data. The grey area indicates power-law temperature profiles with an exponent between –0.3 and –0.9. This profile does not include the emission originating from the accretion shocks that correspond to a local increase in temperature. It also excludes species for which no strong constraints on the size of the emitting region could be derived (see Sect. 3.4.2).

In the text

Fig. 10 Column density profile (dashed black line) and averaged column density (blue line) vs. radius (r). The blue triangle shows the averaged column density estimate from ACA array observations from Csengeri et al. (2017a). The green circle corresponds to the beam-averaged column density estimated from ATLASGAL (Csengeri et al. 2014). The yellow square shows the column density estimates roughly at the position of the accretion shocks with ALMA (Csengeri et al. 2019). The error bars indicate a 20% error on these measurements.

In the text

	Fig. 11 Sketch of G328.2551-0.5321 ignoring projection effects. The full source structure is represented on the left, and we zoom into the envelope of the protostar on the right.
In the text

	Fig. 12 Column densities (upper panel), abundances relative to H2 (middle panel), and relative column densities to methanol (lower panel) of the light molecules and methanol. Only the main isotopologues are represented.
In the text

	Fig. 13 Abundances of the S-bearing molecules relative to H2 in the cold and warm phases. Each molecule is represented by one colour.
In the text

	Fig. 14 Abundance ratios of OCS and H2S relative to SO2 for each source for the layers at Τ = 100 Κ for G328.2551-0.5321 and the sample from Herpin et al. (2009) ordered according to their assumed age sequence, G328.2551-0.5321 being the youngest and IRAS 18162-2048 the most evolved.
In the text

	Fig. 15 Comparison of the deuteration in G328.2551-0.5321 (blue circles) and the dark cloud TMC1 (yellow triangles) and the hot corino IRAS 16293–2422 (pink diamonds). See the caption of Table 5 for the details about the references. The grey line represents the maximum deuteration in G328.2551-0.5321.
In the text

Fig. 16 Line area ratio for the pairs of isotopologues (left panel). The ratios expected at the galactocentric distance of G328.2551-0.5321 from galactocentric variation measurements are indicated (see Table 6 for references). The error bars represent the errors on the line area given by the Gaussian fitting procedure from CLASS. Right panel: ratio of the measured isotopie ratio obtained from the LTE modelling of the column densities and the ratios expected at the galactocentric distance of G328.2551-0.5321 from galactocentric variation measurements. The error bars represent a 20% error on the column densities.

In the text

Fig. 17 COM column densities (top left panel) and abundances relative to H2 (top right panel), methanol (bottom left panel), and dimethyl ether (bottom right panel). Only the main isotopologues are represented. The cold and warm components of the envelope are described by the blue and red histograms. The two shocks are represented in yellow and green. The upper limits are indicated by the black arrows.

In the text

Fig. 18 Molecular abundances relative to H2 towards IRAS 16293-2422 (Schöier et al. 2002; Maret et al. 2005; Jaber et al. 2014; Taquet et al. 2014; Majumdar et al. 2016) and hot cores from Bisschop et al. (2007) vs. that of G328.2551-0.5321. The diamond represents the COMs in the hot cores. The small symbols are used to show each hot core. The larger symbol shows the averaged value. The blue data points represent the cold component of the envelope, and red represents the warm component of the envelope. The grey shaded region corresponds to a factor of 5 scatter from the one-to-one relation.

In the text

Fig. B.1 Rotational diagrams of the carbon chains CH3CCH, CCH, C-C3H2, and HC3N. HC3N is described by three rotational diagrams showing the three temperatures traced by HC3N. The first rotational diagram is for the narrow component and low temperature. The two following rotational diagram are the two rotational diagrams for the two temperature components covered with the larger line width. Finally, we show the vibrationally excited state of HC3N, which is a different entry and represents the column density of the molecule present in this vibrationally excited state.

In the text

	Fig. B.2 Rotational diagrams of the O-bearing molecules HCO+ H213CO, and H2CCO.
In the text

	Fig. B.3 Rotational diagrams of the N-bearing molecules HCN, HN13C, CH2NH, NO, HNCO, and N2H+. The two rotational diagrams of HNCO represent the narrow and the broader components of the Gaussian fit.
In the text

	Fig. B.4 Rotational diagrams of the S-bearing molecules CS, H2 34S, ortho-H2S, OCS, SO2, and SO. The emission from the cold component of the envelope for H2S is described by the rotational diagram of o-H2 34S. The two rotational diagrams of OCS, SO2, and SO represent the narrow and the broader components of the Gaussian fit.
In the text

	Fig. B.5 Rotational diagrams of the S-bearing molecules NS, HCS+, and SO+.
In the text

	Fig. B.6 Rotational diagrams of the deuterated molecules HDCO, DCO+ DCN, and HDO.
In the text

	Fig. B.7 Rotational diagrams of SiO.
In the text

	Fig. B.8 Rotational diagrams of the COMs CH3OH, CH3OCH3, CH3CHO, HC(0)NH2, and CH3OCHO. The three rotational diagrams of methanol and dimethyl ether represent the rotational diagram of the three components (cold and warm components of the envelope and shocks) traced by these molecules.
In the text

	Fig. B.9 Rotational diagrams of the COMs CH3CHO, CH3CN, C2H5CN, C2H3CN, and CH3SH.
In the text

Fig. C.1 Spectrum towards G328.2551-0.5321. The grey filled histogram shows the spectra and the coloured labels the detected transitions of all the light molecules and some transitions of the COMs. The green line represents the LTE model obtained with Weeds. The upper panel for each frequency range shows the labels for the small molecules, and the lower panel corresponds to a zoom in temperature for better visibility of the COMs that are labeled.

In the text