# The GUAPOS project: G31.41+0.31 Unbiased ALMA sPectral Observational Survey. IV. Phosphorus-bearing molecules and their relation to shock tracers

F. Fontani, C. Mininni, M. T. Beltrán, V. M. Rivilla, L. Colzi, I.

Jiménez-Serra, Á. López-Gallifa, Á. Sánchez-Monge, S. Viti

► To cite this version:

F. Fontani, C. Mininni, M. T. Beltrán, V. M. Rivilla, L. Colzi, et al. The GUAPOS project: G31.41+0.31 Unbiased ALMA sPectral Observational Survey. IV. Phosphorus-bearing molecules and their relation to shock tracers. Astronomy & Astrophysics - A&A, 2024, 682, 10.1051/0004-6361/202348219. insu-04822491

HAL Id: insu-04822491

https://insu.hal.science/insu-04822491v1

Submitted on 6 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The GUAPOS project: G31.41+0.31 Unbiased ALMA sPectral Observational Survey

## IV. Phosphorus-bearing molecules and their relation to shock tracers

F. Fontani<sup>1,2,3</sup> ◦, C. Mininni<sup>4</sup> ◦, M. T. Beltrán<sup>1</sup> ◦, V. M. Rivilla<sup>5</sup> ◦, L. Colzi<sup>5</sup> ◦, I. Jiménez-Serra<sup>5</sup>, Á. López-Gallifa<sup>5</sup>, Á. Sánchez-Monge<sup>6,7</sup>, and S. Viti<sup>8,9</sup>

Received 10 October 2023 / Accepted 21 November 2023

## ABSTRACT

**Context.** The astrochemistry of the important biogenic element phosphorus (P) is still poorly understood, but observational evidence indicates that P-bearing molecules are likely associated with shocks.

**Aims.** We study P-bearing molecules and some shock tracers towards one of the chemically richest hot molecular cores, G31.41+0.31, in the framework of the project "G31.41+0.31 Unbiased ALMA sPectral Observational Survey" (GUAPOS), which is being carried out with the Atacama Large Millimeter Array (ALMA).

**Methods.** We observed the molecules PN, PO, SO, SO₂, SiO, and SiS through their rotational lines in the spectral range 84.05–115.91 GHz covered by the GUAPOS project.

**Results.** PN is clearly detected, while PO is tentatively detected. The PN emission arises from two regions southwest of the hot core peak, named regions 1 and 2 here, and is undetected or tentatively detected towards the hot core peak. The PN and SiO lines are very similar both in spatial emission morphology and spectral shape. Region 1 is partly overlapping with the hot core and is warmer than region 2, which is well separated from the hot core and located along the outflows identified in previous studies. The SO, SO₂, and SiS emissions are also detected towards the PN-emitting regions 1 and 2, but arise mostly from the hot core. Moreover, the column density ratio SiO/PN remains constant in regions 1 and 2, while SO/PN, SiS/PN, and SO₂/PN decrease by about an order of magnitude from region 1 to region 2, indicating that SiO and PN have a common origin even in regions with different physical conditions. The PO/PN ratio in region 2, where PO is tentatively detected, is ~0.6–0.9, which is in line with the predictions of pure shock models.

**Conclusions.** Our study provides robust confirmation of previous observational evidence that PN emission is tightly associated with SiO and is likely a product of shock chemistry, as the lack of a clear detection of PN towards the hot core allows us to rule out relevant formation pathways in hot gas. We propose the PN-emitting region 2 as a new astrophysical laboratory for shock-chemistry studies.

**Key words.** stars: formation – ISM: molecules – ISM: individual objects: G31.41+0.31

## 1. Introduction

Hot molecular cores (HMCs) are the ‘cradles’ of future high-mass stars ($M \ge 8 M_\odot$) and rich stellar clusters, and represent the chemically richest environment in the local interstellar medium (ISM; Bisschop et al. 2007; Fontani et al. 2007; Belloche et al. 2013; Rivilla et al. 2017). Indeed, many complex (i.e. at least six atoms) astronomical molecules have been detected for the first time towards the well-known HMCs Orion-KL (Genzel & Stutzki 1989) and Sgr B2 (e.g. Bonfand et al. 2019), including complex organic ones (COMs) such as dimethyl ether, ethanol, methyl cyanide, vinyl cyanide, and many others (McGuire 2018, 2022).

Such rich chemistry is triggered by several physical processes. First, HMCs typically harbour one or more deeply embedded protostars, which heat up their surrounding molecular gas up to temperatures of $T \ge 100$ K. In this warm environment, the molecules that formed in the early cold phase in the ice mantles of dust grains sublimate back into the gas phase via thermal and non-thermal desorption (e.g. Garrod et al. 2022), and gas reactions not efficient at low temperature start to proceed and form new, more complex molecules (Garrod & Herbst 2006; Bonfand et al. 2019; Gieser et al. 2019; Barger & Garrod 2020). Second, collimated jets and molecular outflows from protostars can trigger local chemistry typical of shocked gas, such as grain sputtering and non-equilibrium chemistry

<footnote n="1">¹ INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Florence, Italy
e-mail: francesco fontani@inaf.it</footnote>

<footnote n="2">² Centre for Astrochemical Studies, Max-Planck-Institute for Extraterrestrial Physics, Giessenbachstrasse 1, 85748 Garching, Germany</footnote>

<footnote n="3">³ LERMA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, 92190 Meudon, France</footnote>

<footnote n="4">⁴ INAF – Istituto di Astrofisica e Planetologia Spaziali, Via Fosso del Cavaliere 100, 00133 Roma, Italy</footnote>

<footnote n="5">⁵ Centro de Astrobiología (CSIC-INTA), Ctra Ajalvir km 4, 28850 Torrejón de Ardoz, Madrid, Spain</footnote>

<footnote n="6">⁶ Institut de Ciències de l'Espai (ICE, CSIC), Can Magrans s/n, 08193, Bellaterra, Barcelona, Spain</footnote>

<footnote n="7">⁷ Institut d'Estudis Espacials de Catalunya (IEEC), Barcelona, Spain</footnote>

<footnote n="8">⁸ Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="9">⁹ Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>
(Hollenback & McKee 1989; Bachiller 1996; Tychoniec et al. 2021). Therefore, HMCs with a well-known physical structure are excellent astronomical laboratories for studying astrochemistry in a variety of conditions.

In this framework, the aim of the project 'G31.41+0.31 Unbiased ALMA sPpectral Observational Survey' (GUAPOS; Mininni et al. 2020) is to observe the full 3mm spectral window accessible with the Atacama Large Millimeter Array (ALMA) towards the well-known HMC G31.41+0.31 (G31 hereafter), which is located at a heliocentric distance of 3.75 kpc (Immer et al. 2019, see Mininni et al. 2020 for a detailed description of the source). The unbiased survey allows us to observe many transitions of different species and to confirm (or discard) their detection in a robust way. The proper identification of the species and of all their transitions is also important for the robust derivation of physical parameters, such as abundances and excitation temperatures, which are needed in order to constrain chemical models. The GUAPOS project has already provided important constraints on the formation and/or destruction routes of potential prebiotic species such as the three isomers of C₂H₄O₂ (Mininni et al. 2020), peptide-like bond molecules (Colzi et al. 2021), and oxygen- and nitrogen-bearing COMs (Mininni et al. 2023; López-Gallifa et al. 2023). In this work, we focus our analysis on phosphorus-bearing molecules.

Phosphorus (P) is a basic ingredient for life as we know it. P-compounds are unique in their involvement in the formation of large biomolecules thanks to their extreme structural stability and functional reactivity. For these reasons, P is a crucial component – in the form of phosphate – of nucleic acids (RNA and DNA), cellular membranes (phospholipids), and adenosine triphosphate (ATP), the key molecule for energy transfer in cells (see e.g. Macía 2005; Pasek et al. 2017). Despite its low solar abundance relative to hydrogen (P/H ~3 × 10⁻⁷, Asplund et al. 2009) with respect to other biogenic elements, such as C, O, and N, its abundance in living organisms is several orders of magnitude higher than the solar abundance; for example, up to P/H ~10⁻³ in bacteria (e.g. Fagerbakke et al. 1996). Therefore, identifying the main source and reservoir of P in space, and understanding how its compounds form and evolve, how they are transformed and/or conserved in star-forming regions, and are finally delivered to planetary systems are of huge importance in astrophysics.

P is believed to be mainly formed in massive stars by neutron capture on silicon (Si) in hydrostatic neon-burning shells in the presupernova stage (e.g. Koo et al. 2013 and references therein), and also in explosive carbon- and neon-burning layers during supernova explosions (Woosley & Weaver 1995). In the ISM, it is less depleted than previously thought based on observations of P-bearing molecules obtained in massive star-forming regions and evolved stars (Rivilla et al. 2016; Ziurys et al. 2018). Nevertheless, a handful of P-bearing molecular species (eight in total, including the tentative detection of SiP) have been detected in the ISM, and only three of them in star-forming regions, namely PO (Rivilla et al. 2016; Lefloch et al. 2016; Wurmser & Bergner 2022), PN (Ziurys 1987; Fontani et al. 2016; Rivilla et al. 2018; Bergner et al. 2022; Wurmser & Bergner 2022), and PO⁺ (Rivilla et al. 2022). Single-dish observations and interferometric maps of PN and PO, the P-bearing molecules that are the easiest to detect in star-forming regions (Fontani et al. 2016; Mininni et al. 2018; Rivilla et al. 2016, 2018; Bergner et al. 2022; Wurmser & Bergner 2022) and evolved stars (Agúndez et al. 2007; De Beck et al. 2013; Ziurys et al. 2018), agree with the hypothesis that PN and PO are derived from a solid phosphorus carrier, which is based on their spatial association with SiO emission, a tracer

of protostellar outflows and shocks (Lefloch et al. 2016; Mininni et al. 2018; Rivilla et al. 2018; Fontani et al. 2019; Bernal et al. 2021).

Although some theoretical works have begun to determine the key reactions of the chemical network surrounding phosphorus (Fernández-Ruz et al. 2023), the reactions leading to the formation of PN and (especially) PO are far from clear. For example, Rivilla et al. (2020) proposed that to justify the observed abundances, the molecules should not be directly sputtered from the grains but formed through gas-phase photochemistry induced by ultraviolet (UV) photons from the protostar in post-shocked gas. García de la Concepción et al. (2021) proposed that PO can be formed from atomic P reacting with OH, a reaction that is particularly efficient in warm environments, and Jiménez-Serra et al. (2018) also proposed mechanisms involving energetic processing. Moreover, Rivilla et al. (2020) and Bergner et al. (2022) found that PN and PO emission is spatial with low-velocity and not with high-velocity SiO and SO emission, further complicating the simple scenario of pure grain sputtering.

Here, we study the P-bearing molecules detected in GUAPOS together with other shock tracers: the silicon-bearing molecules SiO and SiS, and the sulphur-bearing species SO and SO₂. Our analysis of the aforementioned sulphur- and silicon-bearing species is limited to the regions in the core where emission from phosphorus-bearing molecules is significant. The paper is structured as follows: in Sect. 2 we describe our observations and data reduction. In Sect. 3, we present our observational results, which we then discuss in Sect. 4. Conclusions and future perspectives are given in Sect. 5.

## 2. Observations and data reduction

Observations towards HMC G31 were made with ALMA during Cycle 5 (project 2017.1.00501.S, P.I.: M.T. Beltrán), whereby we obtained an unbiased spectral survey in Band 3 from 84.05 GHz to 115.91 GHz. The frequency resolution is 0.49 MHz, corresponding to a velocity resolution of ~1.6 km s⁻<fnref n="1" /> at 90 GHz. The final angular resolution is 1.′′2 (~4500 au). The primary beam is ~68″ at 84 GHz and ~50″ at 115 GHz. The pointing centre of the observations is RA(J2000) = 18h47m34.′′312 and Dec(J2000) = -01°12′45.′′9. The uncertainty on the flux calibration is ~5%. For more details on the data reduction (calibration, baseline subtraction, cleaning, and line identification), we refer to Mininni et al. (2020). We performed our spectral analysis, as described in Sect. 3.3, with the MAdrid Data CUBe Analysis (MADCUBA<fnref n="1" />; Martín et al. 2019) software.

## 3. Results

### 3.1. Emission morphology

#### 3.1.1. PN, PO, and SiO

Figure 1 shows the map of PN $J = 2-1$ (panel a) integrated in the velocity range with emission above $3\sigma$ rms, namely 85.7–105.8 km s⁻<fnref n="1" />. As reference, the systemic velocity of G31 is 96.5 km s⁻<fnref n="1" />. The emission morphology is extended and located mostly towards two regions, both south of the hot core peak. The emission peak of the region showing the most intense emission

<footnote n="1">¹ MADCUBA is software developed in the Madrid Center of Astrobiology (INTA-CSIC); it enables the user to visualise and analyse single spectra and data cubes: https://cab.inta-csic.es/mADCuba/</footnote>
<figure><img src="image_2.png" /></figure>

**Fig. 1.** Intensity maps of PN and SiO integrated in velocity. (a) PN *J* = 2–1 integrated in the range 85.7–105.8 km s⁻¹. The white contour is the 3σ rms level of the integrated map (σ = 1.14 × 10⁻² Jy km s⁻¹), while the black contours are in steps of 1σ rms. The numbered regions 1 and 2 are those used to extract the spectra of all species analysed. The synthesised beam is in the bottom-left corner. The yellow contour is the 3 mm continuum emission at 16 mJy beam⁻¹, corresponding to 20σ rms. The yellow stars indicate the continuum sources identified by Beltrán et al. (2021). (b) Map of the intensity of SiO *J* = 2–1 integrated in the same velocity range as PN (colour scale). The PN emission regions identified in panel a are highlighted in white. Contours start at the 3σ rms level of the integrated emission (3 × 10⁻² Jy km s⁻¹), and correspond to 3, 15, 30, 50, 80, and 120σ. (c) Spectra of PN *J* = 2–1 extracted from the emission peak in regions 1 (upper panel) and 2 (lower panel). The red dashed horizontal line indicates the 3σ rms level, and the vertical lines show the position in velocity of the line hyperfine components. (d) SiO *J* = 2–1 emission integrated in the velocity ranges 63.5–85.7 km s⁻¹ (blue contours) and 105.8–130.0 km s⁻¹ (red contours). In both cases, the starting contour is the 3σ rms level of the integrated map (1.95 × 10⁻² Jy km s⁻¹ for the red contours, 1.62 × 10⁻² Jy km s⁻¹ for the blue contours), and the step is 25σ. The dashed lines correspond to the six outflows identified by Beltrán et al. (2022) from SiO *J* = 5–4. The grey filled areas correspond to the PN emission regions 1 and 2 identified in panel a. (e) Same as panel d but for SiO *J* = 5–4, obtained at an angular resolution of ~0.22” with ALMA, and already published in Beltrán et al. (2018). The integration velocity intervals are the same as in panel d.

is offset by −1.′′6,−0.′′9 from the phase centre, corresponding to
~6900 au, and is labelled 1 in Fig. 1. The second region has
a complex shape with multiple peaks. The main intensity peak
is offset by −3.′′9,−3.′′3 (corresponding to ~19 000 au) from the
phase centre, and is labelled region 2 in Fig. 1. Towards both

peaks, we extracted the PN *J* = 2–1 spectra, which are shown
in Fig. 1. The intensity of both lines is clearly above the 3σ rms
level. PO is only tentatively detected towards G31, and so we do
not show the integrated map here; we refer to Sect. 3.2.2 for the
analysis of this tentative detection.
Figure 1 also shows the emission of SiO *J* = 2–1 integrated in the same velocity range as PN *J* = 2–1 (85.7–105.8 km s⁻¹, panel b). The emission morphology is in very good agreement with that of PN. In particular, SiO *J* = 2–1 does not peak on the hot core but towards the two PN emission regions. Evidence that PN and SiO are spatially associated was already found from both single-dish studies (Mininni et al. 2018; Rivilla et al. 2018; Fontani et al. 2019; Lefloch et al. 2016) and interferometric studies (Bergner et al. 2022). The maps in Fig. 1 confirm this association very clearly. The emission from PN is less extended than that from SiO, which is perhaps due to insufficient sensitivity. To quantify this difference, we first derived the contours at half maximum of the emission of PN and SiO in region 1. The contours subtend solid angles with equivalent diameters of ~3" and ~4", respectively, and are therefore comparable. Second, we can quantify the expected emission of PN in the envelope detected in SiO by scaling the SiO integrated intensity by the PN/SiO abundance ratio. The average SiO integrated intensity in the envelope surrounding the emission peak 1 is ~0.4 Jy km s⁻¹. Multiplying this value by the lowest expected PN/SiO relative abundance, which is ~0.1 as we show in Sect. 4.3, the expected integrated emission of PN is ~4 × 10⁻² Jy km s⁻¹, which is marginally higher than the 3σ rms level in the PN map, that is 3.5 × 10⁻² Jy km s⁻¹, and is therefore consistent with a non-detection in PN.

Moreover, SiO is also detected at higher blueshifted and redshifted velocities with respect to the velocities detected in PN. We show the integrated emission in the blue and red wings (63.5–87 km s⁻¹ and 102–130 km s⁻¹, respectively) in panel c of Fig. 1. PN is undetected at these high velocities. However, the remarkable agreement between PN and SiO in the velocity range in which PN emission is significant – and where SiO emission is most intense – suggests that the two molecules have a very similar emission structure.

Panel d of Fig. 1 shows maps of the high-velocity wings of the SiO *J* = 5–4 line published in Beltrán et al. (2018). These maps were obtained at an angular resolution of ~0.22", which is approximately five times higher than that of the GUPOS maps. The orientation of the six outflows in the region and the positions of their driving sources were recently improved with ALMA (~0.09"; Beltrán et al. 2022). We can see that the PN emission is entirely to the SW of the sources driving the outflows, and can be associated with four of the six flows identified in Beltrán et al. (2022) at higher angular resolution. Thus, the complex PN emission we see in both regions 1 and 2 is probably due to the superposition of several outflow lobes. The superposition also suggests that PN is mainly associated with the blueshifted SiO outflow lobes (panel d in Fig. 1). Moreover, although region 1 is associated with both redshifted and blueshifted emission, region 2 is only detected in the blueshifted emission. This could be due to the blueshifted lobes being brighter than the red ones, although there could be other reasons for this, such as inhomogeneity of the dense gas in the clump where the hot core and outflows are embedded.

The PN emission is also clearly away from the location of the hot core, and this result is consistent with what has been observed with ALMA in the young stellar objects AFGL5142 and B1 (Rivilla et al. 2020; Bergner et al. 2022). We checked whether or not some PN emission also arises from the hot core by analysing the spectrum extracted towards the 3 mm continuum (Mininni et al. 2020). The spectrum extracted from the 3 mm continuum contour surrounding the hot core (Fig. 1) is shown in Fig. 2. We chose the 3 mm continuum contour at 20σ rms to disentangle the emission of the hot core from that of the

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Integrated PN *J* = 2–1 spectrum extracted from the hot core. The upper panel shows the spectrum extracted from the 3 mm continuum contour shown in Fig. 1. The lower panel shows the spectrum extracted from a beam around the intensity peak of the 3 mm continuum. The red vertical lines indicate the expected velocities of the hyperfine components, and their length is proportional to the relative intensity of the components. The horizontal lines illustrate the 3σ rms level.</figcaption></figure>

ultracompact HII region placed NE of the hot core (Mininni et al. 2020). We also show the spectrum extracted from a beam centered on the continuum intensity peak of the hot core. We detect faint emission at ~100 km s⁻¹ and ~80 km s⁻¹, which is away from the systemic velocity of G31, of namely 96.5 km s⁻¹. The feature at ~100 km s⁻¹ is unlikely to be a nearby line because no other species have been found to emit lines at the frequency of this feature in previous GUPOS works (Mininni et al. 2020; Colzi et al. 2021; López-Gallifa et al. 2023). Moreover, this latter feature is narrower (full width at half maximum (FWHM) of ~6 km s⁻¹) than the lines detected towards the hot core, which have a FWHM broader than 7 km s⁻¹ (Mininni et al. 2020, 2023; Colzi et al. 2021). The hint of an absorption feature nearby (at ~90 km s⁻¹) suggests that this could be partially self-absorbed PN. The fact that these absorption and emission features are both more prominent towards the total hot core than towards the peak could indicate that it is PN emission from the envelope, possibly associated with narrower features, and/or partially self-absorbed in the blue part of the line. However, the candidate absorption feature is too close to the noise level in both spectra, and we therefore conclude that the emission from the hot core is lacking or negligible, as also suggested in López-Gallifa et al. (2023). The feature at ~80 km s⁻¹ could be due to n-C₃H₇CN 32<sub>5,28</sub> – 32<sub>4,23</sub> at 93.9839 GHz and/or gGg'-(CH₂OH)₂ 15<sub>1,14</sub> – 14<sub>2,12</sub> at 93.9818 GHz.

We stress that we identified all the lines in the full GUPOS spectrum towards region 2, and found no lines of other identified species that could contaminate the PN *J* = 2–1 line. Moreover,
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Emission of SO, SiS, and SO₂ integrated in velocity. Upper panels: velocity-integrated emission of, from left to right, SO *N* = 3–2, *N* = 2–1(2), *N* = 2–1(3), and *N* = 4–4 (see Table 1 for the spectral parameters). The integration velocity range is 85.7–105.8 km s⁻¹ in all images to match the velocity interval where the PN emission is detected. The black contours start from the 3σ rms level of the integrated maps, which is, from left to right: 4.8 × 10⁻² Jy beam⁻¹ km s⁻¹; 2.5 × 10⁻² Jy beam⁻¹ km s⁻¹; 3.6 × 10⁻² Jy beam⁻¹ km s⁻¹; 2.7 × 10⁻² Jy beam⁻¹ km s⁻¹, and the step is of 10σ rms. The white contour corresponds to the PN integrated emission, and the numbered regions 1 and 2 are those defined in Fig. 1. The synthesised beam is illustrated in the lower-left corner. To the right of the quantum numbers, we indicate the *E*<sub>up</sub> of the transition. Lower panels: same as the upper panels for the SiS *J* = 5–4 and *J* = 6–5 lines, and the SO₂ *J*<sub>K<sub>a</sub>,K<sub>b</sub></sub> = 3<sub>1,3</sub> – 2<sub>0,2</sub> and *J*<sub>K<sub>a</sub>,K<sub>b</sub></sub> = 10<sub>1,9</sub> – 10<sub>0,10</sub> lines. The 3σ rms level is: 2.9 × 10⁻² Jy beam⁻¹ km s⁻¹ for SiS *J* = 5–4; 2.3 × 10⁻² Jy beam⁻¹ km s⁻¹ for SiS *J* = 6–5; 2.95 × 10⁻³ Jy beam⁻¹ km s⁻¹ for SO₂ *J*<sub>K<sub>a</sub>,K<sub>b</sub></sub> = 3<sub>1,3</sub> – 3<sub>0,2</sub>; 1 × 10⁻² Jy beam⁻¹ km s⁻¹ for SO₂ *J*<sub>K<sub>a</sub>,K<sub>b</sub></sub> = 10<sub>1,9</sub> – 10<sub>0,10</sub>.</figcaption></figure>

all lines detected towards this region were identified. Towards region 1, the identification of all species is more difficult because the full spectrum is very rich in lines. However, the clean detection towards region 2 and the lack of candidate contaminating lines from other species towards the hot core (Fig. 2) makes a contamination from other lines in region 1 also unlikely. The identification of all species in regions 1 and 2, and their comparison, goes beyond the scope of this paper, and will be the subject of a forthcoming work.

3.1.2. SO, SiS, and SO₂

Figure 3 shows the integrated intensity maps of the SO, SiS, and SO₂ transitions observed in GUAPOS and detected towards the PN-emitting regions. An excitation analysis of all transitions of these species detected towards the hot core is presented in López-Gallifa et al. (2023). The SO emission arises mostly from the hot core, and shows some extended emission to the southwest, with a secondary peak towards the PN-emitting region 2 clearly visible in transitions *N* = 3–2 and *N* = 2–1(3). Therefore, overall, the bulk of the SO emission has a different morphology from that of both PN and SiO. The emission from both SiS and SO₂ is also dominated by the hot core, but there is a clear secondary peak coincident with region 2. The increase in the SO, SiS, and SO₂ integrated intensity towards the hot core is naturally explained if warm or hot gas-phase chemistry is boosting the formation of sulphur-bearing species. However, the SiS lines are also likely contaminated. Indeed, SiS *J* = 5–4 is partly blended in the hot core with a line at 90.7697 GHz, which is unidentified so far, and the SiS *J* = 6–5 is blended with CH₃COCH₃ at 108.9237 GHz

and with an unidentified line detected in the hot core at 108.9216 GHz. In Sect. 3.2, we describe how such contaminations are clear closer to the hot core in region 1, and disappear in region 2.

In summary, from the integrated intensity maps of all species, the SiO emission is the most similar to the PN emission. This is further illustrated in Fig. 4, where we plot the pixel-by-pixel ratio between the PN integrated intensity map and those of SiO *J* = 2–1, SiS *J* = 5–4, SO *N* = 2–1(2), and SO₂ *J*<sub>K<sub>a</sub>,K<sub>b</sub></sub> = 3<sub>1,3</sub> – 2<sub>0,2</sub>. The ratio is nearly constant (within a factor two) for SiO/PN in both regions 1 and 2, while region 1 shows an increase in the ratios SiS/PN, SO/PN, and SO₂/PN by up to a factor of 10 towards the centre of the hot core. This shows that, furthermore, both PN and SiO are formed predominantly in the outflows or in their cavities, while the sulphur-bearing species are produced mostly in the hot core.

3.1.3. ²⁹SiO and ³⁴SO

Figure 5 shows the integrated intensity maps of the lines of the less abundant but clearly detected isotopologues. Because these latter are expected to be optically thin, they should better illustrate the morphology of the molecular emission compared to the main isotopologues. Contrary to the main isotopologue, the ²⁹SiO emission peaks on the hot core. However, this line is likely contaminated by two transitions: CH₃OCHO 21<sub>5,16</sub> – 21<sub>4,17</sub> at 85.761876 GHz, and C₂H₅CN *J*<sub>(K<sub>1</sub>,K<sub>2</sub>)</sub> = 11<sub>2,10</sub> – 11<sub>1,11</sub> at 85.760502 GHz. The centroid velocity of this line is expected to be displaced by about –4.5 km s⁻¹ from ²⁹SiO *J* = 2–1. Because both CH₃OCHO and C₂H₅CN peak towards the hot core, as is
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Maps of integrated intensity ratios. The plots show the ratio between the PN integrated intensity and that of, from left to right, SiO *J* = 2–1, SiS *J* = 5–4, SO *J* = 2–1(2), and SO₂ *J*<sub>K<sub>d</sub>,K<sub>b</sub></sub> = 3<sub>1,3</sub> – 2<sub>0,2</sub>. We show only the emission included in or close to the two PN-emitting regions 1 and 2, as the PN integrated emission goes rapidly to zero elsewhere. The star in the top left corner indicates the phase centre, and the ellipse in the bottom left corner in the first panel is the GUAPOS synthesised beam. The black and white contours correspond to the 3 mm continuum emission and the PN integrated emission, respectively, as shown in Fig. 1.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Velocity-integrated emission of less abundant isotopologues. We show, from left to right, <sup>29</sup>SiO *J* = 2–1, <sup>34</sup>SO *N* = 3–2, *N* = 2–1(2), and *N* = 2–1(3). The integration velocity range is 85–101 km s<sup>-1</sup> for <sup>29</sup>SiO, which is the range in which PN and <sup>29</sup>SiO are both detected. For <sup>34</sup>SO, we used: 81–105 for the *N* = 3–2 line, 88–104 km s<sup>-1</sup> for the *N* = 2–1(2) line, and 85.7–105.8 km s<sup>-1</sup> for the *N* = 2–1(3) line. The black contours start at the 3σ rms level of the integrated intensity maps, which is, from left to right, 1.6 × 10<sup>-2</sup> Jy beam<sup>-1</sup> km s<sup>-1</sup>; 1.8 × 10<sup>-2</sup> Jy beam<sup>-1</sup> km s<sup>-1</sup>; 1.7 × 10<sup>-2</sup> Jy beam<sup>-1</sup> km s<sup>-1</sup>; and 8.7 × 10<sup>-3</sup> Jy beam<sup>-1</sup> km s<sup>-1</sup>, and the step is 10σ rms. The white contour corresponds to the 3σ rms level of the integrated map of PN, and the numbered regions 1 and 2 are those defined in Fig. 1. The synthesised beam is illustrated in the bottom-left corner.</figcaption></figure>

the case for the other COMs detected in the region, the fact that the emission peak in the map is on the hot core is likely due to this contamination.

The emission morphologies of the <sup>34</sup>SO lines are very sim-
ilar to those of their main isotopologues: the *N* = 3–2 and *N* =
2–1(3) transitions both peak towards the hot core but have a clear
secondary peak towards the PN region 2. The *N* = 2–1(2) line
is mostly concentrated on the hot core and the secondary peak
is barely visible, perhaps due to insufficient sensitivity. There-
fore, the maps of the less abundant isotopologues of SO confirm
the morphology of the main one, while that of SiO is different
due to blending with a line arising from the hot core. However,
we stress that a possible partial contamination in region 1 is
also possible in the <sup>34</sup>SO lines. In particular, the *N* = 2–1(2)
line could be contaminated by nearby lines of CH<sub>3</sub>COOH and
CH<sub>3</sub>COCH<sub>3</sub> detected in the hot core. This would explain why
the map of <sup>34</sup>SO in this line is more concentrated towards the hot
core than the others. The *N* = 3–2 and *N* = 2–1(3) could also
be contaminated by nearby lines detected, but unidentified, in
the hot core. In summary, the analysis of all these lines towards
region 1 should therefore be regarded with caution given the
strong contaminations.

3.2. Spectra

3.2.1. Main isotopologues

Figures 6 and 7 show the spectra of the transitions of the main isotopologues reported in Table 1, extracted at the two emitting regions of PN *J* = 2–1 illustrated in Fig. 1. The spectra are in brightness temperature units. The conversion from flux-density units to temperature units was performed using Eq. (1) of Fontani et al. (2021), which provides an average brightness temperature over the angular surface of each region.

The spectral profiles of SiO and PN are compared in Fig. 8,
and appear similar in both regions 1 and 2. This similarity is
especially apparent in the central velocity channels, while it is
less clear at high velocities (≥102 km s⁻¹, and ≤87 km s⁻¹),
where the wings detected in SiO disappear in PN. This is likely
due to insufficient sensitivity in this velocity range for PN. The
detection of PN in a narrower velocity range with respect to SiO
has also been found through ALMA observations in the high-
mass protostar AFGL5142 (Rivilla et al. 2020) and the low-mass
protostar B1-a (Bergner et al. 2022), as well as in single-dish
studies (Lefloch et al. 2016; Mininni et al. 2018; Rivilla et al.
2018; Fontani et al. 2019).
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Spectra extracted from the regions detected in PN *J* = 2–1. We show the transitions (from top to bottom) of PN (*J* = 2–1), SiO (*J* = 2–1 and *J* = 5–4) and SiS (*J* = 5–4 and *J* = 6–5). At the top of each column, we give the number of the extraction region according to the labels in Fig. 1. The blue dashed vertical lines correspond to the reference LSR velocity (96.5 km s⁻¹; Mininni et al. 2020). The black vertical lines in the panels of the SiS *J* = 5–4 line correspond to an unidentified line detected in the hot core at 90.7697 GHz, and those of the *J* = 6–5 transition correspond to CH₃COCH₃ at 108.9237 GHz and to an unidentified line at 108.9216 GHz, both detected in the hot core (Mininni et al. 2023). The red curves superimposed on the lines represent the best LTE fit of the analysed molecules obtained with MADCUBA. For PN, we show the fit obtained fixing *T*ex to that of SiO, that is 26 K (see Sect. 3.3.2).</figcaption></figure>

Both SiS lines appear displaced (redshifted) in velocity
towards region 1 with respect to PN, SiO, and SO (Fig. 6). Such
a velocity shift is not seen in region 2, where the SiS velocity
peak in both lines is placed very close to that of PN, SiO, and
SO. We think that, in region 1, the SiS *J* = 5–4 line is contami-
nated by an unidentified line clearly detected towards the hot core
at 90.7697 GHz, and the SiS *J* = 6–5 line is contaminated by
emission of CH₃COCH₃ at 108.9237 GHz arising from the hot
core, and by an unidentified line at 108.9216 GHz, both clearly
detected towards the hot core (Mininni et al. 2023; Colzi et al.
2021).

The profiles of the SO transitions *N* = 2–1(2) and *N* = 3–2
are, overall, similar to PN in region 1 (Fig. 7), although the inten-
sity peaks are slightly different: PN *J* = 2–1 and SO *N* = 3–2
peak at ~94 km s⁻¹, while SO *N* = 2–1(2) peaks at the G31 sys-
temic velocity, that is 96.5 km s⁻¹. In region 2, the SO lines are
narrower than those of PN and SiO, suggesting that SO could
be dominated by more quiescent material, maybe surrounding
the outflow lobes, towards this position rather than by shocked
material.

Finally, the SO₂ lines have a profile similar to that of SO
in region 2, while in region 1 there are some differences: the
Jₖₐ,ₖ = 101,9 – 100,10 transition has a profile similar to that

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** As in Fig. 6, but showing the spectra of SO₂ (*J* = 3–2 and *J* = 10–10) and SO (*J* = 2–1, *J* = 3–2, and *J* = 4–4). We keep the spectra of PN for reference in the top panels, as in Fig. 6. The black vertical lines in the spectrum of SO₂ *J* = 3–2 are SO₂ *J* = 16–15 (Table 1) and ¹³CH₃OH *J* = 13–14 at 104.0266 GHz.</figcaption></figure>

of SO, while the 3<sub>1,3</sub> – 2<sub>0,2</sub> transition shows a blended pro-
file with the SO<sub>2</sub> 16<sub>2,14</sub> – 15<sub>3,13</sub> line (undetected in region 2)
and with <sup>13</sup>CH<sub>3</sub>OH 13<sub>–3,11</sub> – 14<sub>–2,13</sub> at 104.02655 GHz, which
is detected towards the hot core (Mininni et al. 2023) and not
towards region 2.

3.2.2. Tentative detection of PO

The four lines of PO listed in Table 1 extracted from regions 1
and 2 are shown in the spectra of Fig. 9. We used MADCUBA to
fit the four transitions simultaneously. The fitting method, which
assumes local thermodynamic equilibrium (LTE) conditions,
is fully described in Sect. 3.3.1. We have used the spectro-
scopic entry 047507 of the Cologne Database for Molecular
Spectroscopy (CDMS), which is based on the laboratory mea-
surements of Bailleux et al. (2002) and Kawaguchi et al. (1983),
and the dipole moment determined by Kanata et al. (1988).

In region 1 (upper panel in Fig. 9), a detection of PO can-
not be ruled out but it is difficult to claim due to the spectral
complexity and the contamination by emission lines from other
molecules (Fig. 9). In particular, the two PO transitions in the
right panel of Fig. 9, namely the transitions at 109.2062 GHz
and 109.2812 GHz, are blended with transitions of ethyl formate
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Spectral comparison between PN, SiO, and SO. The PN *J* = 2–1 spectra (grey histogram) extracted from regions 1 and 2 in Fig. 1 are superimposed on SiO *J* = 2–1 (red histograms), SO *N* = 2–1(2) (magenta histograms), SO *N* = 3–2 (blue histograms), <sup>29</sup>SiO *J* = 2–1 (orange histograms), <sup>34</sup>SO *N* = 2–1(2) (purple histograms), and <sup>34</sup>SO *N* = 3–2 (green histograms) extracted from the same regions. The PN line intensity scale has been appropriately multiplied to perform a consistent spectral comparison with SiO and SO. In each plot, the white vertical bars are the theoretical positions in velocity of the PN *J* = 2–1 hyperfine components (Table 1), with the strongest one being at the systemic velocity of 96.5 km s<sup>-1</sup>.</figcaption></figure>

(C<sub>2</sub>H<sub>5</sub>OCHO) and ethylene glycol (aGg-CH<sub>2</sub>OH)<sub>2</sub>, respectively, as already found in source W51 (Rivilla et al. 2016). The lines contaminating the PO transitions in the left panel are not easily identifiable, as a proper determination of all contaminating species would require the study of the full GUAPOS spectrum towards this position. This goes beyond the scope of this paper and will be the subject of a forthcoming work.

In region 2 (lower panel in Fig. 9), which is less affected by nearby lines, the synthetic PO spectrum is consistent with a tentative detection. If we consider the peak intensity, the four transitions are below the 3σ rms noise, but if we consider the S/N of the line integrated intensity, the two strongest transitions, those at 108.9984 GHz and 109.2062 GHz, are consistent with a detection. The S/N of the integrated intensity can be computed from the expression $\int T_B dV / [\sigma \times \sqrt{\Delta V / \text{FWHM}} \times \text{FWHM}]$, where $\Delta V$ is the spectral resolution in velocity. Plugging the line

integrated intensity (0.51 K km s<sup>-1</sup>), the 1σ rms (0.02 K), and the FWHM obtained with MADCUBA (12 km s<sup>-1</sup>, see Table 2) into this expression, we obtain an S/N of ~6. The spectrum towards this position was smoothed in velocity by a factor 2 to increase the sensitivity. The tentative detection is supported by the fact that the emission is most likely detected towards the two lines expected to be the strongest in the quadruplet (i.e. at 108.9984 GHz and 109.2062 GHz).

### 3.2.3. Less abundant isotopologues

The spectra of the less abundant isotopologues, namely <sup>29</sup>SiO, <sup>34</sup>SO, and <sup>33</sup>SO, extracted from regions 1 and 2 are shown in Fig. 10. We clearly detected <sup>29</sup>SiO *J* = 2–1 and <sup>34</sup>SO *N* = 2–1(1), 3–2, and 4–4 in both regions. The detection of <sup>34</sup>SO *N* = 2–1(2) and 4–4 in region 2 is likely but difficult to
<figure><img src="image_10.png" /><figcaption>Fig. 9. Spectra containing the PO lines in the GUAPOS observations. Upper panel: spectrum extracted towards region 1. The red curve indicates the PO lines fitted according to the best-fit parameters in Table 2 (assuming $T_{ex}=26$ K). The horizontal red dashed line corresponds to the $3\sigma$ rms level. Lower panel: same as top panel for region 2, using the fit shown in Table 2, where $T_{ex}=12$ K is assumed. The observed spectra at the original frequency resolution are shown in light grey, while the spectra smoothed by a factor of 2 are shown in dark grey.</figcaption></figure>

confirm because it is very close to the noise level. The line profiles of the $^{34}$SO lines are overall similar in both FWHM and centroid velocity to the corresponding $^{32}$SO ones (compare Figs. 6 and 10).

The line profile of $^{29}$SiO $J = 2-1$ is similar to the $^{28}$SiO profile in the red wing, but this line is relatively much brighter than $^{28}$SiO in the blue wing. As already noted in Sect. 3.1.1, this excess emission is likely due to contamination from $C_2H_5CN$ emission, the expected centroid velocity of which is indicated in Fig. 10.

Figure 8 shows the PN $J = 2-1$ line profile superimposed on the detected lines of $^{29}$SiO and $^{34}$SO. For $^{34}$SO, we use the two unblended lines $N = 2-1(2)$ and $N = 3-2$. In region 2, the PN and $^{29}$SiO line profiles are very similar because $^{29}$SiO is likely less blended. The $^{34}$SO lines are very similar to PN in region 1, and are narrower in region 2. As already noted for the main isotopologue SO lines, this could be due to the fact that the SO emission is associated with more quiescent material towards this position.

We did not clearly detect any $^{33}$SO line, although towards region 1 a tentative detection of the $J = 3-2$ transition is possible, as shown in Fig. 10. Assuming that the line is optically thin in both less abundant isotopologues, we can verify whether or not the tentative detection in $^{33}$SO $N = 3-2$ is consistent with the expected line intensity: towards region 1 the $^{34}$S/$^{33}$S ratio is ~5. Because the relative isotopic ratio $^{34}$S/$^{33}$S is ~5.6 (Lodders 2003), the measured ratio is consistent with a tentative detection, but the $^{33}$SO line also suffers from severe blending with nearby transitions, and we therefore refrain from analysing it further.

### 3.3. Molecular column densities and excitation temperatures

In this section, we derive the molecular column densities, $N_{tot}$, of the analysed species. We first derive the $N_{tot}$ and excitation temperature, $T_{ex}$, of the molecules detected in multiple lines, namely SO, $^{34}$SO, $SO_2$, SiS, and SiO (Sect. 3.3.1). Then, using the excitation temperature derived from these species, we compute $N_{tot}$ of PN and PO (Sect. 3.3.2).

#### 3.3.1. SO, $^{34}$SO, $SO_2$, SiS, and SiO

The fit to the spectra of SO, $^{34}$SO, $SO_2$, SiS, and SiO shown in Figs. 6, 7, and 10, and the derivation of the physical and spectral parameters were performed with the Spectral Line Identification and LTE Modelling (SLIM) tool of MADCUBA. Through its AUTOFIT function, SLIM produces the synthetic spectrum that best matches the data assuming LTE. The input parameters are: $T_{ex}$, $N_{tot}$, radial velocity of the source (V), line FWHM, and angular size of the emission ($\theta_S$). AUTOFIT assumes that V, FWHM, $T_{ex}$, and $\theta_S$ are the same for all transitions fitted simultaneously. These input parameters have all been left free except $\theta_S$, for which we can safely assume that the emission fills the telescope beam, as can be seen in Figs. 3 and 5. The results are shown in Table 2.

Let us first present the best-fit parameters of the S-bearing species. For $^{34}$SO, the fit converged in region 1, leaving V, FWHM, and $N_{tot}$ free, but we had to fix $T_{ex}$. The best-fit $T_{ex}$ obtained by visual inspection is 58 K. The resulting column density in region 1 is ~1.1 × 10$^{15}$ cm$^{-2}$. In region 2, the fit converged, leaving all the parameters free apart from the filling factor, as explained above. The resulting $N_{tot}$ is ~8 × 10$^{13}$ cm$^{-2}$, and $T_{ex}$ is 25 K. For SO, in region 1 the fit could not converge, leaving all parameters free and using all lines. We therefore did not use the $N = 3-2$ transition in the first instance, which has a spectral shape that suggests an excessively high optical depth. We then fitted the remaining three transitions fixing FWHM to 9 km s$^{-1}$, which corresponds to the best-fit value derived from $^{34}$SO. The resulting best-fit column density and $T_{ex}$ are 1 × 10$^{16}$ cm$^{-2}$ and 60 K, respectively. In region 2, the fit converged leaving all the parameters free, and the resulting column density and $T_{ex}$ are ~1.2 × 10$^{15}$ cm$^{-2}$ and 31 K, respectively. The $SO_2$ excitation temperature in both regions is consistent within the uncertainties with the values measured from SO (~75 K in region 1 and ~27 K in region 2), and the column densities are ~9.8 × 10$^{15}$ cm$^{-2}$ in region 1 and ~1.6 × 10$^{14}$ cm$^{-2}$ in region 2.

The fit of the lines of the Si-bearing species SiS and SiO (Fig. 6) converged leaving all parameters free except the filling factor in both regions. However, as discussed in Sect. 3.2,
<table><caption>**Table 1.** Spectral parameters of the analysed transitions.</caption><thead><tr><th>Mol.</th><th>Transition</th><th>Rest freq.<br>(GHz)</th><th>log[A<sub>ij</sub>]<br>(s<sup>-1</sup>)</th><th>E<sub>up</sub><br>(K)</th></tr></thead><tbody><tr><td rowspan="6">PN</td><td>N = 2–1, J = 2–2</td><td>93.9782</td><td>-5.6144</td><td>6.8</td></tr><tr><td>N = 2–1, J = 1–0</td><td>93.9785</td><td>-5.2677</td><td>6.8</td></tr><tr><td>N = 2–1, J = 2–1</td><td>93.9798</td><td>-5.1373</td><td>6.8</td></tr><tr><td>N = 2–1, J = 3–2</td><td>93.9799</td><td>-5.0123</td><td>6.8</td></tr><tr><td>N = 2–1, J = 1–2</td><td>93.9808</td><td>-6.5687</td><td>6.8</td></tr><tr><td>N = 2–1, J = 1–1</td><td>93.9823</td><td>-5.3926</td><td>6.8</td></tr><tr><td rowspan="4">PO</td><td>J = 5/2 – 3/2,<br>Ω = 1/2,</td><td></td><td></td><td></td></tr><tr><td>F = 3–2, l = e</td><td>108.9984</td><td>-4.6712</td><td>8.4</td></tr><tr><td>F = 2–1, l = e</td><td>109.0454</td><td>-4.7160</td><td>8.4</td></tr><tr><td>F = 3–2, l = f</td><td>109.2062</td><td>-4.6689</td><td>8.4</td></tr><tr><td>F = 2–1, l = f</td><td>109.2812</td><td>-4.7143</td><td>8.4</td></tr><tr><td rowspan="2">SiO</td><td>J = 2–1</td><td>86.8470</td><td>-4.5335</td><td>6.3</td></tr><tr><td>J = 5–4<sup>(a)</sup></td><td>217.1050</td><td>-3.2843</td><td>31.3</td></tr><tr><td><sup>29</sup>SiO</td><td>J = 2–1</td><td>85.7592</td><td>-4.5500</td><td>6.2</td></tr><tr><td rowspan="2">SiS</td><td>J = 5–4</td><td>90.7716</td><td>-4.9241</td><td>13.1</td></tr><tr><td>J = 6–5</td><td>108.9243</td><td>-4.6799</td><td>18.3</td></tr><tr><td rowspan="4">SO<sup>(b)</sup></td><td>N = 2–1, J = 2–1</td><td>86.0940</td><td>-5.2799</td><td>19.3</td></tr><tr><td>N = 3–2, J = 2–1</td><td>99.2999</td><td>-4.9488</td><td>9.2</td></tr><tr><td>N = 4–4, J = 5–4</td><td>100.0296</td><td>-5.9656</td><td>38.6</td></tr><tr><td>N = 2–1, J = 3–2</td><td>109.2522</td><td>-4.9665</td><td>21.1</td></tr><tr><td rowspan="4"><sup>34</sup>SO<sup>(b)</sup></td><td>N = 2–1, J = 2–1</td><td>84.4107</td><td>-5.3056</td><td>19.2</td></tr><tr><td>N = 4–4, J = 5–4</td><td>96.7818</td><td>-6.0061</td><td>38.1</td></tr><tr><td>N = 3–2, J = 2–1</td><td>97.7153</td><td>-4.9695</td><td>9.1</td></tr><tr><td>N = 2–1, J = 3–2</td><td>106.7432</td><td>-4.9970</td><td>20.9</td></tr><tr><td rowspan="2"><sup>33</sup>SO<sup>(c)</sup></td><td>N = 2–1, J = 2–1</td><td>85.2379</td><td>-5.2928</td><td>19.3</td></tr><tr><td>N = 3–2, J = 2–1</td><td>98.4936</td><td>-4.9591</td><td>9.1</td></tr><tr><td rowspan="3">SO<sub>2</sub><sup>(d)</sup></td><td>J<sub>K<sub>a</sub>,K<sub>b</sub></sub> = 3<sub>1,3</sub> – 2<sub>0,2</sub></td><td>104.0294</td><td>-4.2264</td><td>7.7</td></tr><tr><td>J<sub>K<sub>a</sub>,K<sub>b</sub></sub> = 16<sub>2,14</sub> – 15<sub>3,13</sub></td><td>104.0336</td><td>-5.4987</td><td>138</td></tr><tr><td>J<sub>K<sub>a</sub>,K<sub>b</sub></sub> = 10<sub>1,9</sub> – 10<sub>0,10</sub></td><td>104.2393</td><td>-3.7708</td><td>55</td></tr></tbody></table>

**Notes.** All parameters are taken from the Cologne Database for Molecular Spectroscopy (CDMS: https://cdms.astro.uni-koeln.de/cdms/portal/; Endres et al. 2016) <sup>(a)</sup>From observations in Beltrán et al. (2022). <sup>(b)</sup>For simplicity, in the text and figures, we label the SO and <sup>34</sup>SO transitions only with the *N* quantum number. We add the *J* of the upper level for the two transitions with the same *N* quantum numbers. <sup>(c)</sup>We list only the two (expected-to-be) strongest hyperfine transitions, *F* = 7/2 – 5/2 and *F* = 9/2 – 7/2. <sup>(d)</sup>We list only the lines detected towards the PN-emitting regions.

the SiS lines in region 1 are both strongly blended with nearby lines, some of which are from unidentified species and therefore could not be fitted simultaneously with SiS. Therefore, we only derived $N_{tot}$ and $T_{ex}$ for region 2, where we obtain $T_{ex}$~13 K and $N_{tot}$~2.6 × 10¹³ cm⁻². For SiO, in region 1, we derive $T_{ex}$~26 K and $N_{tot}$~3.6 × 10¹⁴ cm⁻², and in region 2 $T_{ex}$~12 K and $N_{tot}$~1.8 × 10¹⁴ cm⁻². Finally, we did not estimate the parameters for <sup>29</sup>SiO because the line is too blended with a nearby line of C₂H₅CN, as discussed in Sect. 3.2.

### 3.3.2. PN and PO

For PN, we cannot derive $T_{ex}$ from the data because we have only one transition. We therefore derived the best-fit in a range of $T_{ex}$ based on the minimum and maximum values obtained from the other tracers. These are 26–75 K in region 1 and 12–30 K in region 2. The results are given in Table 2. The PN total column densities are of the order of 10¹³ cm⁻², regardless of the region or temperature assumed. The largest PN column density is measured towards region 1 assuming $T_{ex}$=75 K. We fitted the PN lines also taking the hyperfine structure into account, and derived in all cases FWHMs slightly smaller but consistent within the errors with the values obtained neglecting it.

For PO, as mentioned in Sect. 3.2.2, the PO transitions towards region 1 are significantly contaminated by brighter emission from other species. It was already noted by Rivilla et al. (2016) that the transitions at 109.2062 GHz and 109.2812 GHz are blended with transitions of ethyl formate (C₂H₅OCHO) and ethylene glycol (aGg-CH₂OH)₂, respectively. The contamination of the other two lines is less obvious and would require the analysis of the full GUAPOS spectrum towards this position, which is beyond the scope of this paper. As the presence of PO in this position cannot be confirmed, we derived an upper limit for its column density. We used the values of $T_{ex}$, FWHM, and $V_{LSR}$ derived for PN, and used MADCUBA to increase the value of $N_{tot}$ to the maximum value that is consistent with the observing spectrum. The resulting LTE model is shown in Fig. 9, and the derived upper limit is ~1.2–1.3 (Table 2).

Regarding position 2, which is much less rich in lines, the PO transitions are not contaminated, and seem to be detected, especially the two brightest transitions. We fitted them with MADCUBA, using the same $T_{ex}$, FWHM, and $V_{LSR}$ derived for PN. The resulting LTE fit and $N_{tot}$ are shown in Fig. 9 and Table 2, respectively. We find a PO/PN ratio of ≤1, which is consistent with the results found in some low-mass star-forming regions (Wurmser & Bergner 2022), but is lower than the typical PO/PN ratios measured in high-mass star-forming regions (Rivilla et al. 2016, 2018) and in the comet 67P/Churyumov-Gerasimenko, where the PO/PN ratio is at least 10 (Rivilla et al. 2020).

### 3.3.3. Opacity of SiO and SO lines

MADCUBA estimates the line opacities as explained in Martín et al. (2019), and provides column densities already corrected for optical-depth effects. The optical depth at the line centroid velocity, $\tau_{max}$, computed with MADCUBA for each molecule is listed in Table 2. In case of multiple transitions, we list $\tau_{max}$ of the line with the highest opacity. The values shown in Table 2 indicate optically thin lines in most cases except SiO, for which $\tau_{max}$ is around 1, and SO, for which $\tau_{max}$ is ~0.3.

For both SiO and SO lines, the opacity at line peak can also be computed from the less abundant isotopologues. Some lines of these isotopologues have the same quantum numbers as those of the main isotopologue (see Table 1). Their excitation temperatures should therefore also be comparable and the line opacity can be derived in an alternative way. If the emitting region is also the same for the two isotopologues, from the equation of radiative transfer, one can demonstrate that the ratio between the line brightness temperature of two isotopologues depends only on their relative abundance, and on the optical depth of the main isotopologue. Let us consider the case of SiO: the line intensity
<table><caption>**Table 2.** Best-fit spectral and physical parameters obtained with MADCUBA.</caption><thead><tr><th>Region</th><th>Molecule</th><th>FWHM<br>(km s<sup>-1</sup>)</th><th>V<br>(km s<sup>-1</sup>)</th><th>N<sub>tot</sub><br>(cm<sup>-2</sup>)</th><th>T<sub>ex</sub><br>(K)</th><th>τ<sub>max</sub></th></tr></thead><tbody><tr><td rowspan="8">1</td><td><sup>34</sup>SO</td><td>9.0(0.9)</td><td>95.7(0.4)</td><td>1.1(0.1) × 10<sup>15</sup></td><td>58</td><td>0.04</td></tr><tr><td>SO</td><td>9.0</td><td>96.0(0.5)</td><td>1.0(0.2) × 10<sup>16</sup></td><td>60(10)</td><td>0.34</td></tr><tr><td>SO<sub>2</sub></td><td>12.0(0.6)</td><td>94.8(0.3)</td><td>9.8(1.5) × 10<sup>15</sup></td><td>75(10)</td><td>0.08</td></tr><tr><td>SiS<sup>(a)</sup></td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>SiO</td><td>12.6(0.4)</td><td>95.1(0.2)</td><td>3.6(0.1) × 10<sup>14</sup></td><td>26.4(0.2)</td><td>0.9</td></tr><tr><td>PN</td><td>9.7(0.9)</td><td>94.8(0.4)</td><td>1.3(0.1) × 10<sup>13</sup></td><td>26</td><td>0.01</td></tr><tr><td>PO</td><td>9.7(0.9)</td><td>94.8(0.4)</td><td>3.1(0.3) × 10<sup>13</sup></td><td>75</td><td>0.01</td></tr><tr><td>PO</td><td>9.7</td><td>94.8</td><td>&lt; 1.7 × 10<sup>13</sup></td><td>26</td><td>—</td></tr><tr><td>PO</td><td>9.7</td><td>94.8</td><td>&lt; 3.8 × 10<sup>13</sup></td><td>75</td><td>—</td></tr><tr><td rowspan="8">2</td><td><sup>34</sup>SO</td><td>8.6(0.9)</td><td>93.8(0.4)</td><td>8(2) × 10<sup>13</sup></td><td>25(10)</td><td>0.02</td></tr><tr><td>SO</td><td>7.9(0.5)</td><td>94.8(0.4)</td><td>1.2(0.3) × 10<sup>15</sup></td><td>31(12)</td><td>0.3</td></tr><tr><td>SO<sub>2</sub></td><td>6.2(0.6)</td><td>93.8(0.3)</td><td>1.6(0.2) × 10<sup>14</sup></td><td>27(3)</td><td>0.02</td></tr><tr><td>SiS</td><td>12.3(0.9)</td><td>93.3(0.4)</td><td>2.6(0.3) × 10<sup>13</sup></td><td>13(1)</td><td>0.02</td></tr><tr><td>SiO</td><td>11.5(0.6)</td><td>93.0(0.3)</td><td>1.8(0.2) × 10<sup>14</sup></td><td>11.6(0.6)</td><td>0.7</td></tr><tr><td>PN</td><td>11.7(1.2)</td><td>94.2(0.5)</td><td>1.5(0.1) × 10<sup>13</sup></td><td>12</td><td>0.01</td></tr><tr><td>PO</td><td>12.2(1.2)</td><td>94.3(0.5)</td><td>1.6(0.1) × 10<sup>13</sup></td><td>30</td><td>0.01</td></tr><tr><td>PO</td><td>11.7</td><td>94.2</td><td>~8.5(1.4) × 10<sup>12</sup></td><td>12</td><td>0.005</td></tr><tr><td>PO</td><td>12.2</td><td>94.3</td><td>~1.3(0.2) × 10<sup>13</sup></td><td>30</td><td>&lt;0.001</td></tr></tbody></table>

**Notes.** All parameters were obtained with MADCUBA as described in Sect. 3.3.1. The numbers in parentheses are the uncertainties. The quantities without uncertainties (in parentheses) are fixed in the fit. (a) Both lines are strongly blended with nearby lines, some of which are from unidentified lines.

ratio is given by

$$ \frac{T_p^{28}}{T_p^{29}} \sim \frac{1 - \exp^{-\tau_{28}}}{1 - \exp^{-\tau_{28}/X[28/29]}}, \qquad (1) $$

where $T_p^{28}$ and $T_p^{29}$ are the brightness temperatures of the two lines, $\tau_{28}$ is the optical depth of the main isotopologue, and $X[28/29]$ is the $^{28}$Si/$^{29}$Si relative abundance ratio. If $\tau_{28} \ll 1$, Eq. (1) states that $\frac{T_p^{28}}{T_p^{29}} \sim X[28/29]$, which means that the temperature ratios should be equal to the expected isotopic ratio. For SO, an equation similar to Eq. (1) is valid for the $^{32}$S/$^{34}$S and $^{33}$S/$^{34}$S ratios.

The reference solar values for $^{28}$Si/$^{29}$Si and $^{32}$SO/$^{34}$SO are 19.7 (Anders & Grevesse 1989) and 22.5 (Lodders 2003), respectively. Inspection of Fig. 8 shows that, at the velocity where the $^{28}$SiO line peaks (i.e. ~94.8 km s$^{-1}$), the intensity ratio between the $^{28}$SiO and the $^{29}$SiO lines is $^{28}$SiO/$^{29}$SiO~8 in region 1 and $^{28}$SiO/$^{29}$SiO~20 in region 2. This indicates that in region 2, where the isotopic ratio is very close to the solar one, both lines are likely optically thin. In region 1, the computed $^{28}$SiO/$^{29}$SiO ratio gives $\tau_{28} \sim 2.35$, but this is an upper limit because of the blending of $^{29}$SiO with $C_2H_5CN$. Indeed, the value computed by MADCUBA is $\tau_{28} \sim 1$. In region 2 MADCUBA provides $\tau_{28} \sim 0.7$, which is not consistent with optically thin lines, but is at least consistent with line opacity being smaller than that in region 1.

The $^{32}$SO/$^{34}$SO ratio derived from the $N = 2-1(1)$ line is ~22 in region 2, which is consistent with optically thin emission in both lines, and $^{32}$SO/$^{34}$SO~13 in region 1, which is smaller than the reference value by a factor 1.7. This ratio provides an optical depth of $\tau_{32} \sim 1.2$, which is greater than the 0.34 value provided by MADCUBA, but is still consistent with moderately optically thick lines. As stated in Sect. 3.3, the $N = 3-2$ line is

more optically thick in region 1, but it has been excluded from the fit due to its complex spectral shape that did not allow the fit to converge.

# 4. Discussion

## 4.1. Comparison of spatial emission between PN and shock tracers

The most direct and apparent result of this study is the similar spatial distribution of PN and the SiO $J = 2-1$ bulk-velocity emission. This was already found in the intermediate-/high-mass protostellar object AFGL5142 (Rivilla et al. 2020), as well as in the low-mass protostar B1-a (Bergner et al. 2022) with ALMA observations. Our findings also agree with those of these latter two works in the fact that PN does not arise from the position of the protostar(s) embedded in the hot core. The non-detection of PN towards the hot core could be explained either by insufficient atomic P to form PN, which in turn is abundant along the outflow cavities owing to grain sputtering, or to disruption of the PN molecule in the high-temperature and high-irradiation environment of the hot core. However, other species sensitive to UV photodissociation, such as methanol, are abundant towards the hot core (Mininni et al. 2023). Jiménez-Serra et al. (2018) proposed that $PH_3$ is abundantly produced on grain surfaces via hydrogenation of P during the collapse phase. Then, upon evaporation, it is rapidly (in timescales of $10^4$ yr) converted into PN and PO. However, both molecules are not detected clearly in the hot core, and therefore $PH_3$ does not seem to be the main P-carrier here. The fact that PN is detected only in the shocked regions indicates that the main carrier of P is in the dust cores, where sputtering is needed to (partially or totally) destroy the grains.
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Same as Fig. 6, but for the lines of $^{29}$SiO J = 2–1, $^{33}$SO N = 2–1(2) and N = 3–2, and $^{34}$SO N = 2–1(2), N = 3–2, N = 4–4, and N = 2–1(3). The black solid vertical line in the spectrum of $^{29}$SiO J = 2–1 indicates the expected peak velocity of C$_2$H$_5$CN $J_{K_a,K_b}$ = 11$_{2,10}$ – 11$_{1,11}$. In each panel, the horizontal red dashed line shows the 3σ rms level.</figcaption></figure>

The comparison between PN and SiS is less clear because the SiS emission is heavily contaminated by nearby lines close to the hot core. However, the PN line studied here has an upper level energy that is lower than those of the SiS lines ($E_{up}$ ~ 7 K vs. $E_{up}$ ~ 13 and ~18 K, respectively), and hence could be associated with (slightly) different material. If, and eventually how, the emission changes going to higher excitation lines needs to be investigated through maps of higher-$J$ lines of PN. At present, the map of PN at the highest $J$ at high-angular resolution is the $J = 3-2$ map in Bergner et al. (2022), based on observations made with ALMA towards B1-a. This map appears very similar in morphology to the PN $J = 2-1$ map obtained in the same work. However, the upper energy level of PN $J = 3-2$ is ~13.5 K, which is just a factor ~2 higher than that of $J = 2-1$. It would therefore be interesting to map higher-excitation PN lines in order to check whether or not the intensity peak of the emission changes.

We also propose the PN-emitting region 2 as a new 'hot spot' for shocked material. This region, which is well separated from the hot core, would allow a study of the chemistry of shocked gas

without the influence of the hot core, similarly to other known chemically rich, shocked protostellar spots, such as L1157–B1 (Gueth et al. 1998).

## 4.2. Comparison in velocity between PN and shock tracers

Another finding from our observations is the presence of PN only at relatively low velocities if compared with the velocities attained in the wings of the SiO and SO emission (see e.g. Fig. 6). This is also in agreement with both the high-angular-resolution studies performed towards AFGL5241 and B1-a and the single-dish studies of Mininni et al. (2018) and Fontani et al. (2019). This could either be due to a lack of sensitivity in the high-velocity regime or to the fact that PN is only produced in relatively weak shocks. A similar behaviour is seen in young outflows (e.g. L1448-mm), where molecules such as H$_2$S fade away rapidly for high velocities, while SiO and SO remain very bright at all velocities. One possibility could be the destruction of the onion-shell structure of dust grains where volatile species are detected at lower velocities due to the erosion of the ices and SiO and SO are also seen at high velocities because they are released from the grain cores (Jiménez-Serra et al. 2005). However, as discussed in Sect. 3.2.1, and also in Mininni et al. (2018) and Fontani et al. (2019), the fact that the line profile of PN and SiO is similar in the low(er)-velocity channels suggests that the non-detection of PN at high velocities is more likely due to the lack of sensitivity.

To quantitatively establish this, we computed the expected intensity in the wings of the SiO lines if the SiO abundance were the same as that of PN. The column-density ratio PN/SiO is ≤1/10 (Table 2). Figure 6 indicates that the maximum intensity in the high-velocity wings of SiO J = 2–1 is ~1.5 K in region 1 and ~1 K in region 2. Scaling down these values by the abundance ratio PN/SiO~1/10, which is reasonable because the emission in the wings is optically thin, one would obtain maximum intensities of 0.15 K and 0.1 K, respectively. The 1σ rms noise in the SiO J = 2–1 spectrum in region 1 is 0.1 K, and that in region 2 is 0.033 K. Thus, the intensities in the SiO wings scaled down by the PN/SiO factor would be at most 1.5 and 3 times the rms. Considering that these are upper limits, these results are consistent with a non-detection of the high-velocity wings, as we see in PN.

Concerning the PN peak velocity, in region 1, it is identical to that of SiO, and is consistent within the uncertainties with those of SO, $^{34}$SO, and SO$_2$. In region 2, the PN peak velocity is consistent within the uncertainties with that of all the other tracers except SiO, but the difference between the two values is just ~1.2–1.3 km s$^{-1}$, which is comparable to the velocity resolution of the observations.

## 4.3. Column density comparisons

Figure 11 shows the total molecular column density ratios between PN and the other species studied in this work. We plot the ratios for regions 1 and 2 calculated assuming two $T_{ex}$ as described in Sect. 3.3.2. The SO/PN, $^{34}$SO/PN, and SO$_2$/PN ratios are higher in region 1 by about one order of magnitude than in region 2, while the SiO/PN ratios are consistent in the two regions within the uncertainties. Because $N_{tot}$ of PN remains almost constant (around 10$^{13}$ cm$^{-2}$) in both 1 and 2, the different ratios between the two regions arise mostly from the decrease in $N_{tot}$ of SO, $^{34}$SO, and SiS by an order of magnitude. In Sect. 3.3.3, we conclude that the SiO J = 2–1 line could be affected by optical depth values higher than that provided by
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Column density ratios between molecular species and PN. The different molecules are identified by different colours, and the measurements obtained in regions 1 and 2 are separated on the x-axis as indicated. For each molecule, we plot two values of $N_{\text{tot}}/N_{\text{tot}}[\text{PN}]$, which correspond to the two $N_{\text{tot}}[\text{PN}]$ obtained in the $T_{\text{ex}}$ velocity interval given in Sect. 3.3.2 (Table 2).</figcaption></figure>

MADCUBA towards region 1. The upper limit on $\tau$ estimated this way is ~2.35, instead of ~0.9 obtained with MADCUBA. By correcting $N_{\text{tot}}$ in Table 2 for this different $\tau$, we would obtain $N_{\text{tot}} \sim 6 \times 10^{14} \text{ cm}^{-2}$ towards region 1, and the SiO/PN ratios would become 66 and 24 for $T_{\text{ex}}[\text{PN}] = 12$ and 60 K, respectively. Even in this case, the SiO/PN ratio for $T_{\text{ex}}[\text{PN}] = 12$ would be marginally consistent with the value in region 2.

The SO/PN ratio was found to vary by orders of magnitude in the protostar AFGL5142 (Rivilla et al. 2020). In the latter source, a high-mass protostar is driving a bipolar high-velocity jet surrounded by a cavity, both clearly detected in SO. Several emission spots of PN and PO were detected along the cavity walls, but they were both undetected towards the protostar and the jet. Towards the protostar and the jet, the SO/PN ratio was of the order of 1000 or more, which is similar to the value obtained in G31 towards region 1. Along the outflow cavities, the SO/PN ratio instead drops down to ~70–200, which is consistent with the values we measure in region 2. As discussed in Sect. 3.3.3, the SO lines are affected by non-negligible opacities. This means that the derived total column densities of SO could be lower limits, and so the SO/PN column density ratios could be even higher. The ³⁴SO/PN should not be affected by high optical depth effects, and even in this case the column density ratio drops by about an order of magnitude from region 1 to region 2. In summary, our study confirms that PN ans SiO are very selective tracers of outflow cavities, unlike SO and SO₂.

Finally, the PO/PN ratio in region 2, where PO is tentatively detected, is ~0.6–0.9, and in region 1, the PO/PN upper limit is ~1.2–1.3, depending on the assumed temperature. These ratios, although in line with previous measurements on low- and high-mass protostars (Rivilla et al. 2016; Lefloch et al. 2016; Wurmser & Bergner 2022), are lower than those obtained in the outflow spots of AFGL5142 (Rivilla et al. 2020), where $PO/PN \ge 1$ and increases with distance from the protostar. Ratios $PO/PN \le 1$ are found in the theoretical models of Jiménez-Serra et al. (2018) in the case of pure shock models, without the need for a high cosmic-ray ionisation rate, which in turn would be needed to explain $PO/PN > 1$.

## 5. Conclusions

In the context of the GUPOS project, we studied P-bearing molecules towards the HMC G31 at high-angular resolution to investigate their connection with shock chemistry.

- We clearly detected the PN $J = 2-1$ transition and several SO, SO₂, SiO, and SiS rotational lines. PO lines are tentatively detected.

- The integrated intensity maps indicate that the emission of PN arises from two regions, labelled regions 1 and 2 here, both southwest of the hot core peak, where four of the six outflows detected in Beltrán et al. (2022) are placed. PN is not detected towards the hot core, although region 1 is partly overlapping with it. This allows us to rule out important formation pathways in hot gas.

- The PN and SiO emissions are very similar in morphology and spectral shape, both having two strong emission peaks towards regions 1 and 2, while all sulphur-bearing species emit predominantly from the hot core.

- We propose that the PN-emitting region 2 is a 'hot spot' for shocked material well separated from the hot core, and can be used to study the chemistry of shocked gas without the influence of the hot core, similarly to other known chemically rich shocked regions powered by protostellar objects.

- We derive excitation temperatures in the range of ~26–75 K in region 1, and in the range ~12–30 K in region 2. The column density ratios of all species with respect to that of PN decrease by about an order of magnitude from region 1 to region 2, except the SiO/PN ratio, which is constant within the uncertainties in both regions, further indicating a common origin of the two species.

- We derive a (tentative) column density ratio $PO/PN \sim 1$ in region 2, which is in line with a pure shock model that does not need high cosmic ray ionisation rates.

An interesting follow up of our study will be to map transitions of PN at higher excitation to test whether these lines trace different (e.g. inner) material. Moreover, observing more PN lines will allow us to derive $T_{\text{ex}}$ for PN as well, allowing us, in turn, to compute a more accurate estimate of $N_{\text{tot}}$.

*Acknowledgements.* We thank the anonymous referee for their careful reading of the article and their useful comments. C.M. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 program, through the ECOGAL Synergy grant (grant ID 855130). V.M.R. has received support from the project RYC2020-029387-I funded by MCIN/AEI/10.13039/501100011033, and from the Consejo Superior de Investigaciones Científicas (CSIC) and the Centro de Astrobiología (CAB) through the project 20225AT015 (Proyectos intramurales especiales del CSIC). I.J.-S. and L.C. acknowledge financial support through the Spanish grant PID2019-105552RB-C41 funded by MCIN/AEI/10.13039/501100011033. I.J.-S., L.C. and V.M.R. acknowledge also financial support through the Spanish grant PID2022-136814NB100 funded by MCIN/AEI/10.13039/501100011033 and by "ERDF A way of making Europe". S.V. acknowledges support from the European Research Council (ERC) Advanced grant MOPPEX 833460. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00489.S and ADS/JAO.ALMA#2017.1.00501.S.ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. TheJoint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.

## References

- Agúndez, M., Cernicharo, J., & Guélin, M. 2007, ApJ, 662, L91

- Anders, E., & Grevesse, N. 1989, GeCoA, 53, 197

- Asplund, M., Grevesse, N., Sauval, A., & Jacques, S. P. 2009, ARA&A, 47, 481

- Bachiller, R. 1996, ARA&A, 34, 111
Bailleux, S., Bogey, M., Demuynck, C., Liu, Y., & Walters, A. 2002, J. Mol. Spectrosc., 216, 465

Barger, C. J., & Garrod, R. T. 2020, ApJ, 888, 38

Belloche, A., Müller, H. S. P., Menten, K. M., et al. 2013, A&A, 559, A47

Beltrán, M. T., Cesaroni, R., Rivilla, V. M., et al. 2018, A&A, 615, A141

Beltrán, M. T., Rivilla, V. M., Cesaroni, R., et al. 2021, A&A, 648, A100

Beltrán, M. T., Rivilla, V. M., Cesaroni, R., et al. 2022, A&A, 659, A81

Bergner, J. B., Burkhardt, A. M., Öberg, K. I., et al. 2022, ApJ, 927, 7

Bernal, J. J., Koelemay, L. A., & Ziurys, L. M. 2021, ApJ, 906, 55

Bisschop, S. E., Jørgensen, J. K., van Dishoeck, E. F., et al. 2007, A&A, 465, 913

Bonfand, M., Belloche, A., Garrod, R. T., et al. 2019, A&A, 628, 27

Colzi, L., Rivilla, V. M., Beltrán, M. T., et al. 2021, A&A, 653, 129

De Beck, E., Kamiński, T., Patel, N. A., et al. 2013, A&A, 558, 132

Endres, P., Schlemmer, S., Schilke, P., Stutzki, J., & Müller, H. S. P. 2016, J. Mol. Spectrosc., 327, 95

Fagerbakke, K. M., Heldal, M., & Norland, S. 1996, Aquat. Microb. Ecol., 10, 15

Fernández-Ruz, M., Jiménez-Serra, I., & Aguirre, J. 2023, ApJ, 956, 47

Fontani, F., Pascucci, I., Caselli, P., et al. M. 2007, A&A, 470, 639

Fontani, F., Rivilla, V. M., Caselli, P., et al. 2016, ApJ, 822, L30

Fontani, F., Rivilla, V. M., van der Tak, F. F. S., et al. 2019, MNRAS, 489, 4530

Fontani, F., Barnes, A. T., Caselli, P., et al. 2021, MNRAS, 503, 4320

García de la Concepción, J., Puzzarini, C., Barone, V., Jiménez-Serra, I., & Roncero, O. 2021, ApJ, 922, 169

Garrod, R. T., & Herbst, E. 2006, A&A, 457, 927

Garrod, R. T., Jin, M., Matis, K. A., et al. 2022, ApJS, 259, 1

Genzel, R., & Stutzki, J. 1989, ARA&A, 27, 41

Gieser, C., Semenov, D., Beuther, H., et al. 2019, A&A, 631, 142

Gueth, F., Guilloteau, S., & Bachiller, R. 1998, A&A, 333, 287

Hollenback, D., & McKee, C. F. 1989, ApJ, 342, 306

Immer, K., Li, J., Quiroga-Nuñez, L. H., et al. 2019, A&A, 632, A123

Jiménez-Serra, I., Martín-Pintado, J., Rodríguez-Franco, A., & Martín, S. 2005, ApJ, 627, L121

Jiménez-Serra, I., Viti, S., Quénard, D., & Holdship, J. 2018, ApJ, 862, 128

Kanata, H., Yamamoto, S., & Saito, M. 1988, J. Mol. Spectrosc., 131, 89

Kawaguchi, K., Saito, S., & Hirota, E. 1983, J. Chem. Phys., 79, 629

Koo, B.-C., Lee, Y.-H., Moon, D.-S., et al. 2013, Science, 342, 1346

Lefloch, B., Vastel, C., Viti, S., et al. 2016, MNRAS, 462, 3937

Lodders, K. 2003, ApJ, 591, 12201

López-Gallifa, Á., Rivilla, V. M., Beltrán, M. T., et al. 2023, MNRAS, submitted Macía, E. 2005, Chem. Soc. Rev., 34, 691

Martín, S., Martín-Pintado, J., Blanco-Sánchez, C., et al. 2019, A&A, 631, 159

McGuire, B. A. 2018, ApJS, 239, 17

McGuire, B. A. 2022, ApJS, 259, 30

Mininni, C., Fontani, F., Rivilla, V. M., et al. 2018, MNRAS, 476, L39

Mininni, C., Beltrán, M. T., Rivilla, V. M., et al. 2020, A&A, 644, A84

Mininni, C., Beltrán, M. T., Colzi, L., et al. 2023, A&A, 677, A15

Pasek, M. A., Gull, M., & Herschy, B. 2017, Chem. Geol., 475, 149

Rivilla, V. M., Fontani, F., Beltrán, M. T., et al. 2016, ApJ, 826, 161

Rivilla, V. M., Beltrán, M. T., Cesaroni, R., et al. 2017, A&A, 598, A59

Rivilla, V. M., Jiménez-Serra, I., Zeng, S., et al. 2018, MNRAS, 475, L30

Rivilla, V. M., Drozdovskaya, M. N., Altwegg, K., et al. 2020, MNRAS, 492, 1180

Rivilla, V. M., García De La Concepción, J., Jiménez-Serra, I., et al. 2022, FrASS, 9, 9288

Tychoniec, L., van Dishoeck, E. F., van't Hoff, M. L. R., et al. 2021, A&A, 655, A65

Wakelam, V., & Herbst, E., 2008, ApJ, 680, 371

Woosley, S. E., & Weaver, T. A. 1995, ApJS, 101, 181

Wurmser, S., & Bergner, J. B. 2022, ApJ, 934, 153

Ziurys, L. M. 1987, ApJ, 321, L81

Ziurys, L. M., Schmidt, D. R., & Bernal, J. J. 2018, ApJ, 856, 169