The complex organic molecular content in the L1517B starless core

ABSTRACT

Recent observations of the pre-stellar core L1544 and the younger starless core L1498 have revealed that complex organic molecules (COMs) are enhanced in the gas phase towards their outer and intermediate-density shells. Our goal is to determine the level of chemical complexity towards the starless core L1517B, which seems younger than L1498, and compare it with the other two previously studied cores to see if there is a chemical evolution within the cores. We have carried out 3 mm high-sensitivity observations towards two positions in the L1517B starless core: the core’s centre and the position where the methanol emission peaks (at a distance of ∼5000 au from the core’s centre). Our observations reveal that a lower number of COMs and COM precursors are detected in L1517B with respect to L1498 and L1544, and also show lower abundances. Besides methanol, we only detected CH₃O, H₂CCO, CH₃CHO, CH₃CN, CH₃NC, HCCCN, and HCCNC. Their measured abundances are ∼3 times larger towards the methanol peak than towards the core’s centre, mimicking the behaviour found towards the more evolved cores L1544 and L1498. We propose that the differences in the chemical complexity observed between the three studied starless cores are a consequence of their evolution, with L1517B being the less evolved one, followed by L1498 and L1544. Chemical complexity in these cores seems to increase over time, with N-bearing molecules forming first and O-bearing COMs forming at a later stage as a result of the catastrophic depletion of CO.

astrochemistry, ISM: abundances, ISM: molecules

1 INTRODUCTION

In astrochemistry, complex organic molecules (or COMs) are defined as carbon-bearing compounds with at least 6 atoms in their structure (Herbst & van Dishoeck 2009). It was initially believed that COMs could only form on dust grains in the presence of a source of heat via hydrogenation, atom addition and radical–radical reactions (Watanabe & Kouchi 2002; Garrod, Weaver & Herbst 2008). Indeed, COMs were firstly detected in relatively hot sources with T ≳ 100 K, such as massive hot cores (Hollis, Lovas & Jewell 2000; Hollis et al. 2006; Belloche et al. 2008, 2013) or low-mass warm cores (hot corinos; Bottinelli et al. 2004; Jørgensen et al. 2013). However, the detection of several COMs in the gas phase in starless/pre-stellar cores and dark cloud cores (e.g. Marcelino et al. 2007; Öberg et al. 2010; Bacmann et al. 2012; Cernicharo et al. 2012; Vastel et al. 2014; Jiménez-Serra et al. 2016, 2021; Scibelli & Shirley 2020) indicates that there must be another chemical pathway that allows the presence of COMs in the gas phase at temperatures as low as 10 K. In the last years, there have been different proposals to explain the formation of these COMs under these conditions. However, the chemical pathways giving rise to COMs in cold cores is not well-understood yet, and more observations are needed to constrain the proposed models (Rawlings et al. 2013; Vasyunin & Herbst 2013; Balucani, Ceccarelli & Taquet 2015; Vasyunin et al. 2017; Holdship et al. 2019; Jin & Garrod 2020; Punanova et al. 2022). For example, Scibelli et al. (2021) studied the starless core L1521E, claiming that COMs are not only formed by gas-phase reactions, bus also by surface reactions on dust grains. Other authors suggest that cosmic rays induce desorption from icy mantles on dust grains (Redaelli et al. 2021; Sipilä, Silsbee & Caselli 2021). It is remarkable that methanol (CH₃OH), which plays a central role in the formation for larger COMs, is detected systematically in starless cores (Scibelli & Shirley 2020), and thought to form on the surface of dust grains, partially returning into the gas phase upon reactive desorption (Garrod et al. 2006; Garrod, Wakelam & Herbst 2007; Vasyunin et al. 2017).

In order to understand how COMs form under the cold conditions of pre-stellar/starless cores, Jiménez-Serra et al. (2016, 2021) investigated the radial distribution of large COMs as a function of radius in the L1544 and L1498 starless cores. Methanol tends to show a ring-like morphology circumventing the dust continuum emission (see Tafalla et al. 2006; Bizzocchi et al. 2014; Spezzano et al. 2016; Punanova et al. 2022). Since mapping the emission of larger COMs in starless cores requires large amounts of telescope time, Jiménez-Serra et al. (2016, 2021) observed two positions within these cores: the centre, defined by the position of the dust peak, and the location where methanol peaks. The latter position is representative of an outer, intermediate-density shell located at radii between ∼4000 and ∼11 000 au from the core centre in L1544 and L1498, respectively. Several COM precursors have been found toward both cores such as tricarbon monoxide (CCCO) and cyanoacetylene (HCCCN), but more complex molecules like acetaldehyde (CH₃CHO), methyl formate (CH₃OCHO), or dimethylether (CH₃OCH₃) have only been detected towards L1544, which is at a more advanced stage of evolution than L1498 (see Jiménez-Serra et al. 2016, 2021). However, N-bearing COMs are detected towards L1498 with abundances close to those measured towards L1544.

In this work, we present a similar study carried out towards another starless core, L1517B, believed to be at an even earlier stage of evolution than L1498 based on the deuterium fractionation and the absence of infall motions (Crapsi et al. 2005). The goal is to compare chemical complexity measured towards L1517B with that observed in L1544 and L1498, and to establish if the observed trend of increasing chemical complexity is due to evolution (see Jiménez-Serra et al. 2016, 2021). These cores, located in the Taurus molecular cloud complex, show different observational signatures that provide information about their stage of gravitational collapse. L1544 is classified as a pre-stellar core because it shows clear evidence of gravitational collapse towards its innermost regions as probed by the emission of N₂H⁺ (Caselli et al. 2002; Redaelli et al. 2019). Its high deuterium fractionation (⁠|$N_\mathrm{N_2D^+} / N_\mathrm{N_2H^+} = 0.23 \pm 0.04$|⁠), CO depletion factor (f_CO = 14 ± 3; Crapsi et al. 2005; Redaelli et al. 2019), and central H₂ density (⁠|$n_\mathrm{H_2}\sim 10^7$| cm⁻³; Caselli et al. 2022), are also consistent with this idea. L1498 presents signatures of infall motions in the outer envelope of the core as revealed by the asymmetries observed in the line profiles of CS (Tafalla et al. 2004) but its deuterium fractionation, CO depletion factor and central H₂ density are lower than those measured towards L1544 (⁠|$N_\mathrm{N_2D^+} / N_\mathrm{N_2H^+} = 0.04 \pm 0.01$|⁠, f_CO = 7.5 ± 2.5, and |$n_\mathrm{H_2}\simeq 9.4 \times 10^4$| cm⁻³; Tafalla et al. 2004; Crapsi et al. 2005). For L1517B, the CO depletion factor and deuterium fractionation are similar to those of L1498, f_CO = 9.5 ± 2.8 and |$N_\mathrm{N_2D^+} / N_\mathrm{N_2H^+}$| = 0.06 ± 0.01 (Crapsi et al. 2005). However, although L1517B shows a slightly higher central H₂ density (⁠|$n_\mathrm{H_2}\simeq 2.2 \times 10^5$| cm⁻³), the absence of infall motions either towards the innermost regions or towards the envelope (Tafalla et al. 2004), suggests that L1517B is at a younger evolutionary stage than L1498 and L1544. Hence, L1517B, located at a distance of 159 pc (Galli et al. 2019), would be the dynamically youngest of the three cores.

The paper is organized as follows. In Section 2, we describe the observations carried out towards L1517B. In Section 3, we present the results of the analysis of the COMs and COM precursor emission and report the values of the derived excitation temperatures, column densities, and molecular abundances. Section 4 compares the abundances obtained towards L1517B with those measured towards L1498 and L1544. We also compare our results with those reported by Nagy et al. (2019) and Scibelli et al. (2021) towards the young L1521E starless core. This core has been found to be rich in COM emission despite its youth. In Section 5, we present the modelling of the COMs and COM precursor chemistry of the L1517B core, and compare the model predictions with the observations. In Section 6, we compare the column densities ratios between the N-bearing species HC₃N and CH₃CN derived towards L1517B, L1498, L1544, and L1521E with those obtained in protostellar systems, protoplanetary discs, and comets, and discuss the observed discrepancies. Finally, in Section 7, we summarize our conclusions.

2 OBSERVATIONS

The observations of the L1517B starless core were carried out from 2020 September 30th to October 4th, with the Instituto de Radioastronomía Milimétrica (IRAM) 30-m telescope (Granada, Spain). As for L1544 and L1498, we observed two positions within the L1517B core: the location of the dust peak and the position where the emission of methanol peaks. These two positions have the equatorial coordinates (in J2000 system) α = |$4^\mathrm{h}55^\mathrm{m}17\rm{.\!\!^{\rm {\mathrm{s}}}}6$|⁠, |$\delta$| = 30^○37′44″ for the dust continuum peak, and α = |$4^\mathrm{h}55^\mathrm{m}15\rm{.\!\!^{\rm {\mathrm{s}}}}7$|⁠, |$\delta$| = 30^○38′04″ for the position of the methanol peak (see Tafalla et al. 2004, 2006). The latter is located ∼32 arcsec away from the dust peak (see Fig. 1), which corresponds to ∼5000 au at a distance of 159 pc.

Figure 1.

H₂ column density map obtained from Herschel/spire data at 0.25, 0.35, and 0.50 mm, using the same procedure as the one used for L1544 in Spezzano et al. (2016). Crosses indicate the positions observed in the core: the dust peak (in black) and the methanol peak (in red). The beam size for each filter is 17.9, 24.2, and 35.4 arcsec (0.25, 0.35, and 0.50 mm, respectively), although the images of the two first filters were smoothed to the resolution of 0.50 mm, which is marked with a circle at the bottom-left corner of the image. The positions of the dust and methanol peaks are obtained from Tafalla et al. (2004, 2006). The beam sizes were retrieved from the spire Handbook: https://www.cosmos.esa.int/web/herschel/spire-overview.

The high-sensitivity 3-mm spectra were obtained in frequency-switching mode using a frequency throw of 7.14 MHz. The EMIR E090 receivers were tuned at 84.37 and 94.82 GHz with rejections of ≥10 dB. The observed frequency ranges are shown in Table 1. To identify possible weak spurious features in the observed spectra, we carried out part of the observations by shifting slightly the central frequencies by ±20 MHz (see also Jiménez-Serra et al. 2016). We used the narrow mode of the FTS spectrometer that provided a spectral resolution of 49 kHz, equivalent to 0.13–0.18 km s⁻¹ at 3 mm. Typical system temperatures ranged between 75 and 110 K and the telescope beam size was 22–31 arcsec between 78 and 111 GHz; as the beams are almost Gaussian and the dust and methanol peaks of L1517B are located ∼32 arcsec away (⁠|$\gt 2\, \sigma$|⁠), the contamination between both positions should be negligible (less than 5 percent). The spectra were calibrated in units of antenna temperature, |$T_\mathrm{A}^{*}$|⁠, and converted into main beam temperature, T_mb, by using beam efficiencies of 0.81 at 79–101 GHz and of 0.78 at 102–111 GHz.¹ The root mean square (RMS) noise level of the original observations ranged between 3 and 11 mK for both observing positions (see Table 1), having similar values for each frequency range than the ones obtained by Jiménez-Serra et al. (2016, 2021) for L1544 and L1498.

Table 1.

Frequency ranges, velocity resolution, and RMS noise level of our observations.

Frequency	Resolution	RMS noise (mK)
(GHz)	(km s−1)	Dust peak	Methanol peak
78.2–80.0	0.18	4.8	3.9
81.4–83.3	0.18	4.2	3.5
83.3–85.2	0.17	2.8	3.3
86.7–88.5	0.17	3.0	3.5
93.5–95.7	0.15	4.1	3.4
95.8–96.8	0.15	10.8	10.9
97.1–99.0	0.15	4.1	3.6
99.1–100.9	0.15	3.0	3.5
102.4–104.2	0.14	4.1	4.7
109.2–111.0	0.13	10.1	11.0

Our observations have covered the transitions of both O-bearing and N-bearing COMs and COM precursors, summarized in Table 2. For O-bearing species, we have targeted methanol (CH₃OH), methoxy (CH₃O), tricarbon monoxide (CCCO), ketene (H₂CCO), formic acid (t-HCOOH), acetaldehyde (CH₃CHO), methyl formate (CH₃OCHO), dimethylether (CH₃OCH₃), cyclopropenone (c-C₃H₂O), and propynal (HCCCHO). As N-bearing COMs and precursors, we have observed cyanoacetylene (HCCCN), isocyanoacetilene (HCCNC), vinyl cyanide (CH₂CHCN), acetonitrile (CH₃CN), and methyl isocyanide (CH₃NC).

Table 2.

Most of the COMs and COM precursors transitions covered in our L1517B observations and their derived line parameters.

			Dust peak				Methanol peak
Species	Line	Frequency	Areaa	Linewidth	LSR velocityb	S/Nc	Areaa	Linewidth	LSR velocityb	S/Nc
		(MHz)	(mK km s−1)	(km s−1)	(km s−1)		(mK km s−1)	(km s−1)	(km s−1)
CH3OH	20, 2 → 10, 1 A	96741.371	206.1 ± 1.9	0.278 ± 0.004	5.785 ± 0.002	109	241.1 ± 2.0	0.277 ± 0.005	5.788 ± 0.002	121
	21, 2 → 11, 1 A	95914.310	<7	...	...	...	>7	...	...	...
	21, 1 → 11, 0 A	97582.798	<2.4	...	...	...	<2.4	...	...	...
	21, 2 → 11, 1 E	96739.358	150.4 ± 1.9	0.262 ± 0.008	5.786 ± 0.003	80	181.6 ± 2.0	0.262 ± 0.003	5.785 ± 0.003	92
	20, 2 → 10, 1 E	96744.545	9.7 ± 1.9	0.244 ± 0.012	5.84 ± 0.04	5.1	11.2 ± 2.0	0.23 ± 0.08	5.78 ± 0.03	5.6
	21, 1 → 11, 0 E	96755.501	<7	...	...	...	<7	...	...	...
CH3Od	|$F = 1 \rightarrow 0,\, \Lambda = -1$|	82455.980	<4	...	...	...	<3	...	...	...
	|$F = 2 \rightarrow 1,\, \Lambda = -1$|	82458.252	5.4 ± 1.2	0.41 ± 0.08	5.81 ± 0.03	4.5	4.2 ± 0.9	0.39 ± 0.08	5.82 ± 0.03	4.7
	|$F = 2 \rightarrow 1,\, \Lambda = +1$|	82471.825	5.4 ± 1.2	0.41 ± 0.08	5.81 ± 0.03	4.5	4.2 ± 0.9	0.39 ± 0.08	5.82 ± 0.03	4.7
	|$F = 1 \rightarrow 0,\, \Lambda = +1$|	82524.180	<4	...	...	...	<3	...	...	...
CCCO	10 → 9	96214.813	<7	...	...	...	<7	...	...	...
t-HCOOH	11, 1 → 00, 0	87926.863	<1.9	...	...	...	<2.2	...	...	...
H2CCO	41, 3 → 31, 2 o-	81586.299	11.1 ± 1.3	0.53 ± 0.10	5.89 ± 0.04	8.5	12.6 ± 1.2	0.70 ± 0.25	5.95 ± 0.11	10
	51, 5 → 41, 4 o-	100094.510	9.9 ± 0.9	0.53 ± 0.10	5.89 ± 0.04	11	12.6 ± 1.2	0.70 ± 0.25	5.95 ± 0.11	10
	40, 4 → 30, 3 p-	80832.189	21 ± 4	0.53 ± 0.10	5.89 ± 0.04	5.2	17 ± 4	0.70 ± 0.25	5.95 ± 0.11	4.2
CH3OCHO	72, 6 → 62, 5 A	84454.754	<2.0	...	...	...	<2.5	...	...	...
CH3OCH3	31, 3 → 20, 2 EE	82650.180	<3	...	...	...	<2.4	...	...	...
	41, 4 → 30, 3 EE	99326.000	<2.0	...	...	...	<2.2	...	...	...
CH3CHO	51, 4 → 41, 3 A	98900.944	<3	...	...	...	4.0 ± 0.8	0.21 ± 0.05	5.765 ± 0.018	5.0
	50, 5 → 40, 4 A	95963.459	<7	...	...	...	<6	...	...	...
	41, 3 → 31, 2 A	79150.166	<3	...	...	...	3.9 ± 0.7	0.21 ± 0.05	5.765 ± 0.018	5.6
	41, 3 → 31, 2 E	79099.313	<3	...	...	...	3.9 ± 0.7	0.21 ± 0.05	5.765 ± 0.018	5.6
c-C3H2O	61, 6 → 51, 5	79483.519	<3	...	...	...	<3	...	...	...
HCCCHO	90, 9 → 80, 8	83775.816	<2.0	...	...	...	<2.5	...	...	...
CH3CN	60 → 50	110383.500	6.4 ± 1.9	0.13 ± 0.13	5.902 ± 0.023	3.4	<6	...	...	...
	61 → 51	110381.372	11.8 ± 1.9	0.24 ± 0.08	5.81 ± 0.03	6.2	<6	...	...	...
	62 → 52	110374.989	<6	...	...	...	<6	...	...	...
CH3NC	50 → 40	100526.541	3.6 ± 0.7	0.4 ± 0.3	5.87 ± 0.13	5.1	<2.2	...	...	...
CH2CHCN	90, 9 → 80, 8	84945.988	<2.0	...	...	...	<2.5	...	...	...
	91, 9 → 81, 8	83207.496	<3	...	...	...	<2.4	...	...	...
	100, 10 → 90, 9	94276.625	<3	...	...	...	<2.4	...	...	...
HCCCN	9 → 8	81881.461	693.6 ± 1.0	0.395 ± 0.001	5.824 ± 0.001	694	346.7 ± 0.9	0.383 ± 0.002	5.816 ± 0.001	385
	11 → 10	100076.392	286.5 ± 0.7	0.318 ± 0.001	5.842 ± 0.001	409	96.3 ± 0.8	0.290 ± 0.003	5.837 ± 0.001	120
HCCNC	88 → 78	79484.131	16.1 ± 1.1	0.31 ± 0.15	5.83 ± 0.07	15	6.9 ± 0.9	0.29 ± 0.16	5.86 ± 0.06	7.7
	1010 → 910	99354.250	6.8 ± 0.7	0.31 ± 0.15	5.83 ± 0.07	9.7	<2.3	...	...	...

Note. Line profiles were fitted using madcuba, except for methanol (CH₃OH), cyanoacetilene (HCCCN), and acetonitrile (CH₃CN), where we used class (see Section 3.2 for details).

a

Uncertainties in the line area are calculated as |$\Delta T (\Delta \nu \, \delta v)^{1/2}$|⁠, with ΔT the RMS noise nevel, Δv the line width, and |$\delta v$| the velocity resolution of the spectrum. Similarly, upper limits are calculated as 3 |$\Delta T (\Delta v \, \delta v)^{1/2}$|⁠.

b

LSR stands for local standard rest.

c

This refers to the signal to noise ratio in integrated intensity area. If a certain value has no uncertainty, this means that it had to be fixed so that madcuba could fit the LTE model

d

Hyperfine components of the |$N = 1 \rightarrow 0,\, K=0,\, J = 3/2 \rightarrow 1/2$|transition of CH₃O.

The raw spectra have been analysed and reduced with class, from the package gildas,² as well as with a Python³ pipeline written specifically for this purpose.⁴ This pipeline fits baselines to the spectra using an iterative method that first masks the strongest lines using sigma-clips and then applies rolling medians and rolling averages, interpolating the masked regions with third-order splines. Finally, we used the software madcuba (Martín et al. 2019) to search for the molecular transitions of all targeted species, and to carry out the fitting of the molecular line profiles under the assumption of local thermodynamic equilibrium (LTE). The physical parameters derived in the fitting are the molecular column density (N_obs), excitation temperature (T_ex), linewidth (Δv), and LSR radial velocity (v_LSR). For this, we used the tool slim (Spectral Line Identification and Modelling) of madcuba, employing the Cologne Database for Molecular Spectroscopy (CDMS; Endres et al. 2016) and the Jet Propulsion Laboratory (JPL) molecular catalogue (Pickett et al. 1998). In addition, for the cases of CH₃OH, CH₃CN, and HCCCN, we used the non-LTE code radex (Van der Tak et al. 2007) to do the fitting of the lines.

3 RESULTS

3.1 Detected transitions

Figs 2 and 3 show the observed spectra of some representative transitions of the COM and COM precursors detected towards L1517B, while Table 2 lists all the transitions covered in our observations with their derived line parameters. The targeted transitions are the same as in Jiménez-Serra et al. (2016, 2021), and correspond to those expected to be the brightest for an excitation temperature of ∼10 K.

Figure 2.

Sample of the COM and COM precursor lines observed towards L1517B’s dust continuum peak and their corresponding fits obtained with madcuba (LTE, in red) and radex (non-LTE, in orange). See Sections 2 and 3.1 for more details.

Figure 3.

Sample of the COM lines observed towards L1517B’s methanol peak and their corresponding fits obtained with madcuba (LTE, in red) and radex (non-LTE, in orange). See Sections 2 and 3.1 for more details. Acetaldehyde (CH₃CHO) transitions are quite noisy, but we could actually detect it because we observed several ones. It seems that acetonitrile (CH₃CN) may show a transition, but we could not properly fit the two covered transitions neither with madcuba nor radex; plus, the signal to noise is quite low (below 2 for the shown transition and below 1 for the other one).

Besides methanol, our spectra reveal the detection of other COMs and COM precursor species in L1517B: CH₃O, H₂CCO, CH₃CHO, HCCCN, HCCNC, CH₃CN, and CH₃NC although acetaldehyde (CH₃CHO) was only detected towards the methanol peak and acetonitrile (CH₃CN) and methyl isocyanide (CH₃NC) were only detected towards the dust peak. More complex molecules such as CH₃OCH₃ or CH₃OCHO were not detected within our noise levels. All detected lines lie above the |$3\, \sigma$| level in integrated intensity (area), where |$1\, \sigma$| is calculated as |$\Delta T (\Delta v \, \delta v)^{1/2}$|⁠, where ΔT is the RMS noise level, Δv is the line width, and |$\delta v$| is the velocity resolution of the spectrum (see Table 1). We are confident about the detection of the transitions since their derived radial velocities correspond to the v_LSR of the source (∼5.8 km s⁻¹). Moreover, except for CH₃OH A, CH₃NC, and HCCNC, we have measured at least two transitions above the 3|$\sigma$| level in integrated intensity for the detected species, stressing the identification of the species since the linewidths are narrow (⁠|$\sim \,$|0.3–0.5 km s⁻¹). The level of line confusion at the targeted RMS noise level is also very low, as commonly found in starless/pre-stellar cores. For the non-detections, we provide upper limits to the integrated intensities by using 3 |$\Delta T (\Delta v \, \delta v)^{1/2}$|⁠.

From Figs 2 and 3, we find that the line intensities vary for each position: for the N-bearing species (CH₃CN, HCCNC, and HCCCN), the emission is brighter towards the core’s centre with respect to the methanol peak. On the contrary, for O-bearing species (CH₃OH and H₂CCO) we find that the emission level is similar for both positions. Note that in the cases of CH₃OH, CH₃CN and HCCCN the fits were obtained with radex. In these cases the lines were also fitted with independent Gaussians using class but they are not shown in Figs 2 and 3, as they are independent fits for each transition and there are only made to obtain the parameters of the lines shown in Table 2.

3.2 Molecular column densities and excitation temperatures

3.2.1 LTE analysis

Using the tool slim from madcuba, we performed the LTE fitting to the observed line profiles for each species, obtaining the parameters shown in Table 2. As in most cases, we detected several lines of the same molecular species, this fitting allows us to obtain its column density (N_obs) and excitation temperature (T_ex). These parameters are also shown in Table 3.

However, there were cases in which the madcuba fitting algorithm did not converge. This can be due because the signal-to-noise ratio is not high enough, or also because we only have one transition detected. In these cases, we had to fix one or more parameters, so that the algorithm could converge. Similarly, for non-detections we had to fix the excitation temperature so that madcuba could fit an upper limit for the column density.

In general, when we needed to fix the temperature, we used the fitted value for H₂CCO for oxygen-bearing species, and the fitted one for HCCCN for nitrogen-bearing species, as these are molecules with several transitions which also have high signal-to-noise ratio. These values are, for H₂CCO, 7.7 and 9.7 K for the dust and methanol peaks, respectively, and, for HCCCN, 6.7 and 5.3 K. For cases where we had to fix the central velocity of the line, we used a value of v_LSR = 5.8 km s⁻¹.

For methyl isocyanide (CH₃NC) towards the core centre and isocyanoacetylene (HCCNC) towards the methanol peak, we only detected 1 transition, so we had to fix the excitation temperature to T_ex = 6.7 K for CH₃NC and to T_ex = 5.3 K for HCCNC (which are the corresponding values for HCCCN at the dust and methanol peaks, respectively).

Table 3.

Excitation/kinetic temperatures (T), column densities (N_obs), and abundances (|$\chi _{\mathrm{obs}}\,$|) of COMs and COM precursors towards the dust and methanol peaks in L1517B.

	Dust peak			Methanol peak
Molecule	|$\boldsymbol{T}\, (\mathrm{K})$|	|$\boldsymbol{N}_{\mathrm{\mathbf {obs}}}\, (\mathrm{cm}{}^{-2})$|	|$\boldsymbol{\chi }_{\mathbf {\mathrm{\mathbf {obs}}}}\,$|	|$\boldsymbol{T}\, (\mathrm{K})$|	|$\boldsymbol{N}_{\mathrm{\mathbf {obs}}}\, (\mathrm{cm}{}^{-2})$|	|$\boldsymbol{\chi }_{\mathrm{\mathbf {obs}}}\,$|
CH3OH A	10.0	(4.29 ± 0.11) × 1012	|$1.23_{-0.16}^{+0.21}\times 10^{-10}$|	10.0	(4.96 ± 0.12) × 1012	|$5.2_{-0.5}^{+0.6}\times 10^{-10}$|
CH3OH E	10.0	|$4.59_{-0.25}^{+2.63}\times 10^{12}$|	|$1.5_{-0.3}^{+0.7}\times 10^{-10}$|	10.0	(5.62 ± 0.20) × 1012	|$5.9_{-0.6}^{+0.7}\times 10^{-10}$|
CH3OH	10.0	|$8.9_{-0.3}^{+2.6}\times 10^{12}$|	|$2.7_{-0.5}^{+0.7}\times 10^{-10}$|	10.0	(1.058 ± 0.023) × 1013	|$1.10_{-0.10}^{+0.12}\times 10^{-9}$|
CH3O	10 ± 4	(2.8 ± 0.9) × 1011	|$8.0_{-2.5}^{+2.8}\times 10^{-12}$|	7 ± 4	(1.8 ± 0.6) × 1011	|$1.9_{-0.6}^{+0.7}\times 10^{-11}$|
CH3OCHO	7.7	<9 × 1011	<4 × 10−11	9.7	<1.0 × 1012	<1.5 × 10−10
CH3OCH3	7.7	<1.0 × 1012	<5 × 10−11	9.7	<5 × 1011	<7 × 10−11
CH3CHO	7.7	<2.5 × 1011	<1.2 × 10−11	9.7	(2.1 ± 0.4) × 1011	|$2.2_{-0.4}^{+0.5}\times 10^{-11}$|
t-HCOOH	7.7	<3 × 1012	<1.5 × 10−10	9.7	<1.9 × 1012	<3 × 10−10
c-C3H2O	7.7	<1.5 × 1010	<7 × 10−13	9.7	<5 × 1010	<7 × 10−12
H2CCO	7.7 ± 1.0	(8.4 ± 1.3) × 1011	|$2.4_{-0.4}^{+0.6}\times 10^{-11}$|	10 ± 3	(7.2 ± 1.9) × 1011	|$7.6_{-2.0}^{+2.2}\times 10^{-11}$|
CCCO	7.7	<2.4 × 1011	<1.1 × 10−11	9.7	<1.2 × 1012	<1.8 × 10−10
HCCCHO	7.7	<1.6 × 1011	<8 × 10−12	9.7	<1.9 × 1011	<3 × 10−11
HCCNC	6.7	(2.4 ± 1.0) × 1011	(7 ± 3) × 10−12	5.3	(1.9 ± 0.8) × 1011	|$1.9_{-0.8}^{+0.9}\times 10^{-11}$|
CH2CHCN	6.7	<4 × 1010	<1.8 × 10−12	5.3	<3 × 1010	<5 × 10−12
CH3NC	6.7	(1.6 ± 0.9) × 1010	(5 ± 3) × 10−13	5.3	<2.3 × 1010	<3 × 10−12
CH3CN	10.0	(2.1 ± 0.6) × 1011	|$6.0_{-1.9}^{+2.1}\times 10^{-12}$|	10.0	<2.0 × 1011	<3 × 10−11
HCCCN	10.0	(5.16 ± 0.05) × 1012	|$1.48_{-0.18}^{+0.25}\times 10^{-10}$|	10.0	|$3.72_{-0.13}^{+0.14}\times 10^{12}$|	(3.9 ± 0.4) × 10−10

Note. Temperatures (T) refer to excitation temperatures (T_ex) for all the species except for methanol (CH₃OH), cyanoacetilene (HCCCN), and acetonitrile (CH₃CN), where they refer to kinetic temperatures (T_kin). We used madcuba to derive the molecular parameters from the observations except for methanol, cyanoacetilene, and acetonitrile, where we used radex. Molecular abundances were calculated using an H₂ column density of (3.5 ± 0.5) × 10²² cm⁻² for the dust continuum peak and of (9.6 ± 1.0) × 10²¹ cm⁻² for the position of the methanol peak. For the non-detections and also for some detections, we had to fix the excitation temperature (T_ex) so that madcuba could fit the column density (N_obs).

For ketene (H₂CCO), we calculated the ortho-to-para ratio (we have detected two ortho transitions and one para transition), by carrying out a separate fit of the transitions and by dividing the resulting column densities of the ortho and para species. We obtain an ortho-to-para ratio of |$1.0_{-0.3}^{+0.4}$| for the position of the dust peak and |$0.9_{-0.4}^{+0.5}$| for the methanol peak. This ratio is different to what has been previously observed towards the molecular cloud TMC-1 and the pre-stellar core L1689B, where the ratio was ∼3 (Ohishi et al. 1991; Bacmann et al. 2012). An exact value of 3 would correspond to the statistical ratio due to the nuclear spin degeneracy. However, lower values are also possible for the lowest temperatures, like in the case of formaldehyde (H₂CO), where an ortho-to-para ratio of 1 would indicate that this molecule is formed under thermal equilibrium at a temperature of ∼10 K (Kahane et al. 1984). We note, however, that we have only detected two ortho lines and one para line of ketene, and therefore the derived ortho-to-para ratios can be subject to significant uncertainties.

3.2.2 Non-LTE analysis

For some cases, in particular for methanol (CH₃OH), cyanoacetilene (HCCCN), and acetonitrile (CH₃CN), we have decided to use the non-LTE code radex (Van der Tak et al. 2007) to infer the column density N_obs of the species instead of madcuba. In the case of methanol and acetonitrile, we did it because madcuba could not fit properly all the transitions, and the errors were greater than |$3\, \sigma$|⁠. For the case of cyanoacetilene (HCCCN), we opted for radex because it allows us to estimate the H₂ number density, as we have two transitions. As input parameters for these calculations, we assumed the linewidth obtained with the madcuba fit (∼0.28 km s⁻¹) and a kinetic temperature of 10 K [as inferred for all radii across the L1517B starless core; see Tafalla et al. (2004) and Section 5 below]. Appendix  A describes the procedure used to derive the physical properties and their uncertainties from the modelling.

In the case of CH₃OH, we have to distinguish between the A and E species, making a separate fit for each molecule. We used the collisional coefficients with H₂ given by Rabli & Flower (2010). In addition to deriving the column densities reported in Table 3, we also derived the H₂ number density from the observations of methanol E, as this parameter can be determined by the intensity ratio of the detected lines. Unfortunately, for the case of methanol A we only detected one transition so the column density was estimated by assuming the H₂ number density obtained by Tafalla et al. (2004): 2.20 × 10⁵ cm⁻³ for L1517B core centre and 1.23 × 10⁵ cm⁻³ for the location of the methanol peak. For methanol E, we could derive the H₂ number density, although the uncertainties are high, obtaining values of |$5_{-4}^{+8} \times 10^4$| cm⁻³ for the core’s centre and |$2.0_{-0.6}^{+0.9} \times 10^5$| cm⁻³ for the methanol peak. These values are consistent, within the uncertainties, with those derived by Tafalla et al. (2004), as we have to take into account the possible uncertainties in the determination of the density profiles from the absorption coefficient of the dust. As for the excitation temperatures, they are close to the assumed kinetic temperature of 10 K, consistent with those expected for the derived H₂ densities, although this is biased by the fact that the lowest temperature for which we have collisional rates is 10 K.

For CH₃CN, we fixed the H₂ number density to the values from Tafalla et al. (2004) (as the signal-to-noise ratio was not high enough to get a good fit). For the core’s centre, one of the lines could be fitted within |$2\, \sigma$|⁠, providing a better fit than the one by madcuba. For the methanol peak, we did not detect any transition, so we obtained an upper limit to its column density from radex.

For HCCCN, we have used the two transitions to derive the column density and the H₂ number density towards both positions. We used the collisional coefficients with H₂ given by Faure, Lique & Wiesenfeld (2016). Fig. A3 shows the results obtained from the non-LTE analysis for the core’s centre. As expected the derived H₂ densities and column densities are related. We obtained the HCCCN column densities shown in Table 3 and H₂ densities of |$1.80_{-0.04}^{+0.04} \times 10^5$| and |$5.56_{-0.18}^{+0.19} \times 10^4$| cm⁻³ for the dust and methanol peaks, respectively. Although these results present lower uncertainties than those of methanol, our H₂ densities are still consistent with the prediction by Tafalla et al. (2004) taking into account the beam size of our observations and the uncertainties in the determination of the density from dust emission, as we previously discussed. In the following sections, we use the values from Tafalla et al. (2004) because they are more direct than the derivation through radex, and because the authors assume the same dust emissivity (κ = 0.005 cm² g⁻¹ for 1.2 mm) as that considered by Crapsi et al. (2005) for inferring the CO depletion factor of L1517B, which has been employed in our astrochemical simulations (Section 5).

3.3 Molecular abundances

Once we know the column densities of COMs and COM precursors detected toward L1517B, we can compute the molecular abundances of these species by using the H₂ column densities measured towards the positions of the dust and methanol peaks (Table 3). For the position of the dust peak, the H₂ column density of (3.5 ± 0.5) × 10²² cm⁻² was obtained using the 1.2-mm continuum emission observed by Tafalla et al. (2004) with the mambo 1-mm bolometer array at IRAM 30-m telescope. For the position of the methanol peak, however, we employed the data from the Herschel space telescope at 0.25, 0.35, and 0.5 mm assuming a dust optical depth index of β = 1.5 (Spezzano et al. 2016); the derived H₂ column density is (9.6 ± 1.0) × 10²¹ cm⁻². As discussed in Jiménez-Serra et al. (2021), this is the best procedure to probe the total H₂ column density, respectively, towards the innermost regions (with bolometers at 1.2 mm) and outer shells (with Herschel) in starless cores.

From Table 3, we can see that the abundance of the detected species is enhanced towards the position of the methanol peak with respect to the core’s centre. For methanol (CH₃OH), its abundance is enhanced by a factor of ∼4, while for the rest of detected species in both positions (CH₃O, H₂CCO, HCCNC, and HCCCN) the factor of enhancement is ∼3. Acetaldehyde (CH₃CHO) is only detected in the methanol peak, the ratio would be ≳2. As for acetonitrile (CH₃CN) and methyl isocianide (CH₃NC) we only find lines towards the dust peak, and the ratio would be ≲6. Therefore, for these three molecules, the lower/upper limits for the factor of enhancement toward the methanol peak are consistent with those of the molecules detected in both positions. Additionally, the ratio of the abundances of methanol A and E species toward the dust and methanol peak is, respectively, |$0.91_{-0.32}^{+0.08}$| and |$0.88_{-0.04}^{+0.04}$|⁠, which is similar to what has been observed in the L1544 and L1498 cores (Jiménez-Serra et al. 2016, 2021) and to the expected values of |$\sim \,$|0.7–1.0 (Wirström et al. 2011).

4 COMPARISON WITH OTHER CORES

In this Section, we compare our results obtained for the starless core L1517B with those of the L1544 and L1498 cores (see Jiménez-Serra et al. 2016, 2021) and also with those of the L1521E core (see Nagy et al. 2019; Scibelli et al. 2021).

4.1 Comparison with L1498 and L1544

For L1517B, L1498, and L1544, the same molecular species have been observed toward two positions (i.e. the dust peak and the methanol peak for each of them), which allows us to obtain information about the radial distribution of the abundances of COM and COM precursors within the cores.

Fig. 4 presents the abundances of the COM and COM precursor molecules observed toward these three cores and for the two measured positions. The abundances measured towards the dust and methanol peaks are plotted, respectively, in yellow and red for the L1544 pre-stellar core (Jiménez-Serra et al. 2016), in light and dark blue for the L1498 starless core (Jiménez-Serra et al. 2021), and in light and dark purple for the L1517B starless core (this work).

Figure 4.

Bar plot of the abundances of different COMs and COM precursors measured for the cores L1517B, L1498, and L1544, toward the dust peak (lighter colors) and the methanol peak (darker colours) for each of them. Upper limits are indicated both with arrows at the top of the corresponding bar and with stripes. Abundances for L1544 are taken from Jiménez-Serra et al. (2016, 2021), while for L1498 they are taken from Jiménez-Serra et al. (2021). Note that this is not a cumulative bar plot, so the upper edge of each bar, for both colours (that is, for both the centre and the methanol peak, for each core), indicates the abundance of the corresponding species.

Consistently with L1498 and L1544, the abundances of the species detected towards L1517B tend to be higher toward the position of the methanol peak than towards the dust peak (t-HCOOH towards L1544 is an exception). Fig. 4 also shows that the L1517B core presents a lower number of detections than L1544 and a similar number to L1498. For the N-bearing species, we have detected only HCCCN and HCCNC towards both positions of L1517B from a total of five targeted species, although we have found CH₃CN and CH₃NC towards the dust peak. A significant difference with respect to L1498 and L1544 is that vinyl cyanide (CH₂CHCN) is not detected toward L1517B, showing upper limits to its abundance that are factors ≳10 lower than the abundances measured toward L1498 and L1544. Following the same trend, for O-bearing species only four COMs and COM precursor species (CH₃OH, CH₃O, CH₃CHO, and H₂CCO) out of 10 have been detected toward L1517B, and in the case of CH₃CHO it is only detected towards the methanol peak.

If we compare in detail L1517B and L1498, we can see that the first one has a level of chemical complexity somewhat lower than L1498. First, the abundances of the COMs and COM precursors detected towards L1517B tend to show lower abundances than (or are comparable to) those measured toward L1498. Secondly, if we count the total number of detections for any position, we obtain 13 detections for L1517B and 14 for L1498, which may seem a small difference. However, if we consider the detections and non-detections for each position, the lower level of chemical complexity for the L1517B starless core becomes more apparent.

To quantitatively measure this we have calculated, for each pair of cores, the fraction of molecules with the same state of detection (detected or not detected). That is, for each source, we build a vector with one entry for each molecule and position (15 molecules × 2 positions, 30 entries), whose value is 1 if there is a detection and 0 if not. Then, we build another vector with the same number of entries, that are 1 if the two input entries (of the two cores) have the same value and 0 otherwise. Finally, we sum up all the elements of this vector and divide it by the number of entries, in order to normalize it. We call the resulting number the fractional similarity, whose values lie between 0 and 1. In this way, if for each position and molecule both entries have the same value towards both sources (either a detection or a non-detection), the fractional similarity is 1. On the contrary, if for each position and molecule both entries are different, the resulting similarity is 0. If we compute the fractional similarity between L1517B and L1498, we obtain 0.70 ± 0.03.⁵ This indicates that, although the two cores are similar in terms of number of detections, they are truly different regarding their level of chemical complexity since the similarity, clearly less than 1, is produced by 9 different entries in the vectors (that is, in the detections for each molecule and position).

For the sake of completeness, we have also computed the fractional similarity between L1544 and L1517B, which is 0.57 ± 0.03, and between L1544 and L1498, which is 0.43 ± 0.03. This indicates that the three cores show a significantly different level of chemical complexity in terms of detections and non-detections, and the two more similar cores would be L1498 and L1517B.

We have used two other methods to evaluate the level of chemical complexity towards the three cores. First, we have computed for each source and position the mean molecular mass of the species weighted by the abundances in a logarithmic scale.⁶ This can be viewed as an estimate of the level of complexity of each source and position, as COMs tend to have more atoms with increasing complexity, and thus a greater molecular mass. From Table 4, we find that the lowest value of the mean molecular mass is obtained for L1517B, while the greatest one is for L1544, with L1498 being close to it. This would indicate that L1517B presents the lowest level of chemical complexity of the three cores, followed by L1498 and L1544. Additionally, we also computed for each source and position the geometric mean of the abundance of the species weighted by the molecular mass.⁷ The resulting value represents a measure of the amount of the more complex species in the core, as we are giving larger weights in the mean to more complex species. As shown by Table 4, L1517B shows the lowest mean abundance, followed by L1498 and L1544, which present barely the same quantity. All of this indicates that the level of chemical complexity in the L1517B starless core seems to be lower than that of L1498, and of L1544.

Table 4.

Weighted geometric mean of the abundance with respect to H₂, 〈χ_obs〉, and weighted mean molecular mass, 〈m_molec〉 for the targeted species in L1517B, L1498, and L1544.

Source		|$\boldsymbol{\left\langle \chi _{\mathrm{obs}}\right\rangle }\, \left(\times 10^{-12}\right)$|		|$\boldsymbol{\left\langle m_{\mathrm{molec}}\right\rangle }\, \left(\mathrm{g/mol}\right)$|
L1517B	dust peak	|$7.1_{-1.5}^{+1.5}$|	|$15_{-3}^{+3}$|	|$45.5_{-0.3}^{+0.3}$|	|$45.61_{-0.27}^{+0.20}$|
	meth. peak	|$23_{-5}^{+5}$|		|$45.7_{-0.4}^{+0.3}$|
L1498	dust peak	|$10.5_{-2.1}^{+2.1}$|	|$32_{-5}^{+5}$|	|$46.18_{-0.24}^{+0.18}$|	|$46.25_{-0.17}^{+0.15}$|
	meth. peak	|$53_{-10}^{+10}$|		|$46.32_{-0.24}^{+0.23}$|
L1544	dust peak	|$19_{-4}^{+4}$|	|$30_{-4}^{+4}$|	|$46.8_{-0.4}^{+0.3}$|	|$46.46_{-0.21}^{+0.15}$|
	meth. peak	|$42_{-7}^{+7}$|		|$46.13_{-0.19}^{+0.16}$|

Note. The weights for the mean abundance are the molecular masses for each molecule. The weights for the mean molecular mas are the abundance for each molecule in logarithmic scale. The values on the second columns for each magnitude are the arithmetic means of the magnitudes for both positions. To deal with uncertainty propagation and upper limits we made simulations with statistical distributions (see Appendix  B).

From all these results, we suggest that L1517B is in a less chemically evolved stage than L1498 and L1544. This is supported by the deuterium fractionation of L1517B, similar to that of L1498, but L1517B shows no evidence of infall motions (Tafalla et al. 2004; Crapsi et al. 2005). Additionally, L1517B shows the lowest level of chemical complexity within the three cores, presenting only the simplest COMs and COM precursors and with lower abundances in most of the cases.

If we take into account the possible time evolution between the three cores (with L1517B being the youngest and L1544 the oldest), we propose a scenario in which the N-bearing molecules would form first in starless cores. Then, as the cores evolve, they would accrete gas from the surrounding molecular cloud, becoming denser and yielding a strong depletion of carbon monoxide (CO). At the pre-stellar core stage, the catastrophic depletion of CO takes place, triggering the formation of O-bearing COMs and COM precursors as observed in L1544 (Vasyunin et al. 2017). Therefore, chemical complexity would increase over time. In the next section, we explore this scenario by modelling the chemistry of N-bearing and O-bearing COMs and COM precursors towards L1517B.

4.2 Comparison with L1521E

L1521E is another well-studied starless core located in the Taurus molecular cloud. Like L1517B, L1521E is also considered to be young. Indeed, L1521E shows no evidence of infall motions, and has a central density of 2.7 × 10⁵ cm⁻³ (Tafalla & Santiago 2004). Its CO depletion factor is f_CO = 4.3 ± 1.6 (Nagy et al. 2019), lower than those of L1517B and L1498. Although its deuterium fractionation |$N_\mathrm{N_2D^+} / N_\mathrm{N_2H^+}$| has not been measured yet, several authors agree in the fact that this core is young (e.g. Tafalla & Santiago 2004, Kong et al. 2015, Scibelli et al. 2021). Nagy et al. (2019) observed the molecular line emission towards the dust peak of this core (including some COMs and COM precursors) using the IRAM 30-m telescope. More recently, Scibelli et al. (2021) detected four COMs towards L1521E (CH₃OCHO, CH₃OCH₃, CH₃CHO, and CH₂CHCN) using the 12-m ARO (Arizona Radio Observatory) telescope, with a much larger beam size (∼70 arcsec versus 22–31 arcsec for the IRAM 30-m telescope). This implies that while the observations of Nagy et al. (2019) would cover the location of just the dust peak, the observations of Scibelli et al. (2021) would cover both the dust peak and the methanol ring, as its methanol peak is found 22–29 arcsec away from the core centre (Nagy et al. 2019). We thus refer to the observations of Scibelli et al. (2021) as those performed toward the bulk of the core.

Among the 15 targeted species measured towards L1517B, there are 9 of them that have also been targeted toward L1521E: six towards the dust peak and four for the bulk of the core (with just one molecule in common for both positions). Fig. 5 reports the abundances of the common species observed towards L1517B and L1521E. For the data under the bulk label, we use the values derived by Scibelli et al. (2021) assuming two source size cases: the lower limits of the error bars refer to the values obtained assuming a beam filling factor of 1, while the error bars upper limits refer to the abundances obtained with the best-fitting source size of 35 arcsec. Note, however, that vinyl cyanide (CH₂CHCN) is an exception since its abundance was only derived assuming a beam filling factor of 1 (the source size for this molecule could not be estimated; Scibelli et al. 2021).

Figure 5.

Bar plot of the abundances of different COMs and COM precursors measured for the young starless cores L1517B (in light and dark purple) and L1521E (in light and dark green). For L1517B, the values are as in Fig. 4. For L1521E, the abundances shown in dark green refer to the bulk COM emission observed by Scibelli et al. (2021) towards L1521E with a beam size of ∼70 arcsec. The large beam includes not only the dust peak position but also the methanol ring (peak located at 22–29 arcsec from L1521E’s dust peak). Abundances from L1521’s dust peak are shown in light green and are taken from Nagy et al. (2019). Upper limits are indicated with both arrows and stripes.

In general, there is a good agreement between the abundances measured toward the dust peak position of the two cores (see light purple and light green bars in Fig. 5), with the sole exception of cyclopropenone (c-C₃H₂O). In contrast, the differences are bigger between the bulk abundances of L1521E and the abundances measured towards the dust and methanol peaks of L1517B. However, if we take into account the uncertainties with the error bars upper/lower limits, the abundances of L1521E’s bulk are largely consistent with those measured towards both positions of L1517B, with the exception of vinyl cyanide (CH₂CHCN). A possible explanation for this would be that L1521E is more chemically evolved than L1517B, lying closer to L1498 than to L1517B. Alternatively, although these cores are located in the same molecular cloud complex, the local physical and chemical properties of their environment could differ, causing significant differences in their final chemical composition. Studies of the COM chemical content toward a larger sample of starless cores are needed to establish whether differences in the environment affect the COM chemical evolution in starless cores.

We finally note that calculations of the fractional similarity and of the geometric means of the abundance and of the molecular mass for L1521E, cannot be performed since the number of data available for this core is rather small.

5 MODELLING THE FORMATION OF COMS AND COM PRECURSORS IN L1517B

5.1 The model

We have modelled the chemistry of COMs and COM precursors towards the starless core L1517B by using the 0D gas–grain chemical code monaco (Vasyunin & Herbst 2013; Vasyunin et al. 2017), which has been successfully applied previously to the starless cores L1544 and L1498 (Jiménez-Serra et al. 2016, 2021). Our goal is to compare the observed molecular abundances with those predicted by the model, and to infer the radial distribution of COMs and COM precursors towards L1517B and its age.

The monaco chemical code is a rate equations-based, three-phase (gas, ice surface, and ice bulk) numerical model that provides the evolution of the fractional abundances of atomic and molecular species with time. To obtain the radial distribution of the abundances of COMs and COM precursors in L1517B, the code is run for a grid of distances on the density and temperature profile of L1571B, assuming physical parameters to be constant in time. The physical structure of the L1517B starless core is obtained from the H₂ gas density profile by Tafalla et al. (2004), and the gas kinetic temperature distribution inferred by the same authors using NH₃ observations. The assumed H₂ gas density profile for L1517B is:

$$\begin{eqnarray} n(r) = \frac{2.2 \times 10^5 \mathrm{cm}^{-3}}{1 + \left(\frac{r}{35\,\mathrm{arcsec}} \right)^{2.5}} \mathrm{,} \end{eqnarray}$$

(1)

where r is the angular distance (in arcsec). As for the temperature radial profile, we have used a constant temperature of T = 10 K (see Tafalla et al. 2004). Fig. 6 shows both profiles as a function of distance to the centre of the core. In our model, we assume that the dust and gas temperatures are equal.

Figure 6.

H₂ gas density and temperature profiles used for modelling the chemistry of COMs and COM precursors towards L1517B. These radial profiles are taken from Tafalla et al. (2004). The vertical grey dashed line indicates the distance of the observed methanol peak with respect to the centre of the core.

Given the low temperatures in starless cores (T = 10 K), the precursors of COMs are formed on the surface of dust grains via hydrogenation processes. Once formed, a small fraction of these compounds are non-thermally desorbed from dust grains into the gas phase, where they undergo gas-phase chemical reactions that yield the observed COMs (Vasyunin & Herbst 2013; Vasyunin et al. 2017). As non-thermal desorption processes, the code includes UV photodesorption, cosmic ray-induced desorption (Hasegawa & Herbst 1993) and reactive (chemical) desorption (Minissale et al. 2016). The chemical network for the O-bearing and N-bearing COMs and precursors is the same as the one used for the L1498 starless core (see the details in Jiménez-Serra et al. 2021). As initial abundances, the model employs the results of a simulation of a diffuse cloud model with constant gas density of 100 cm⁻³ and gas temperature of 20 K for 10⁷ yr, using the low metals initial abundances EA1 from Wakelam & Herbst (2008).

5.2 Comparing simulations and observations

Using the initial conditions explained in Section 5.1, we have simulated the chemical evolution of L1517B by running the monaco code for 10⁶ yr. To constrain the age of the core, we have used the CO depletion factor of f_CO = 9.5 ± 2.8 measured by Crapsi et al. (2005) towards L1517B, and the location of the methanol peak at a distance of 5000 au. Interestingly, our model cannot reproduce simultaneously the location of the methanol peak at 5000 au and the CO depletion factor of |$\sim \,$|7–12. In Fig. 7, we show the evolution of the predicted CO depletion factor over time and its observed value. Note that the monaco code considers the telescope beam size when calculating the CO depletion factor. By comparing the model results with the CO and COMs and COM precursor observations, the best agreement is reached after 1.5 × 10⁵ yr. The CO depletion factor at this time is ∼7.3, which is on the lower edge of the estimated interval for the CO depletion factor (see Fig. 7).

Figure 7.

Evolution of the simulated CO depletion factor over time. The blue line and shaded region show the value f_CO = 9.8 ± 2.8, derived by Crapsi et al. (2005), while the grey line indicates the best-fitting age for L1517B.

At this time of chemical evolution, however, the location of the methanol peak in the model does not match the observed one towards L1517B. This is clearly shown in Fig. 8, where we present the radial distribution of the abundance of COMs and COM precursors modelled for L1517B, together with the measured values and upper limits (see dots and vertical arrows). As shown in Fig. 8, the location of the methanol peak in the model is at ∼15 000 au away from the dust peak (vertical dashed line), which is a factor of 3 further away that observed (at |$\sim 5000\,$| au; see dashed–dotted vertical line).

Figure 8.

Radial distribution of the modelled abundances of the COMs and COM precursors observed for L1517B. Colour curves represent the modelled abundances, while dots and vertical arrows represent the observed abundances, and their upper limits, towards the two observed positions in L1517B, the dust and methanol peaks. The position of some of the dots (observed COMs and COM precursor abundances) has been slightly shifted to enhance visibility. The dotted–dashed grey line indicates the position of the observed methanol peak, while the grey dashed line indicates the position of the methanol peak according to the chemical simulations.

The origin of this discrepancy likely arises from the fact that the central H₂ gas density towards L1517B is 2.2 × 10⁵ cm⁻³, i.e. significantly lower than measured towards the L1544 pre-stellar core. If this density were higher, as in L1544 (∼10⁷ cm⁻³; Caselli et al. 2022), the CO depletion factor (9.5 ± 2.8) would be reached earlier in the simulations, enabling a better match between the modelled and the observed location of the methanol peak at 5000 au. This would also help reconciling the age of the L1517B starless core predicted by the model (of ∼1.5 × 10⁵ yr) with the one expected from observations. As already mentioned, L1517B does not show evidence for gas accretion either towards the innermost regions or from the outer envelope, as observed towards L1498 or L1544, which suggests that L1517B is at an earlier stage of evolution, or at least, at a similar evolutionary stage as L1498 (note that they both have similar central H₂ gas densities and deuterium fractionation values; Tafalla et al. 2004; Crapsi et al. 2005). Another possible explanation of this discrepancy would be a prominent asphericity of the core, although from Fig. 1 it seems that L1517B is rather spherical.

For the abundances of COMs and COM precursors, Table 5 compares the values between the modelled and observed abundances towards the positions of the dust and methanol peaks in L1517B. All the modelled and observed abundances agree within a factor of 10 except for CH₃CN in the core’s centre and HC₃N in both positions. Actually, at earlier times of the simulations, where the simulated methanol peak is closer to the observed one, abundances of those two species are closer to the observed ones. In any case, we should take into account that astrochemical models entail intrinsic uncertainties derived from uncertainties from some of the constants used in the simulation, such as the chemical reaction rates (Vasyunin et al. 2004; Wakelam et al. 2005, 2010).

Comparing the predicted age for L1517B (1.5 × 10⁵ yr) with those predicted for L1498 and L1544, we find that, contrary to what we thought, L1517B would not be the youngest core, as the predicted ages for L1498 and L1544 are, respectively, 9.8 × 10⁴ and 1.6 × 10⁵ yr (Jiménez-Serra et al. 2016, 2021). Therefore, according to the chemical ages based on similar initial conditions, L1517B would be more evolved than L1498 and less than L1544. However, although L1517B shows some concentration of H₂ at its centre, there is no clear sign of current contraction motions towards this core. This, together with the observed level of chemical complexity, seems to suggest that L1517B is less evolved that L1544 and L1498. Note also that our chemical modelling does not include the dynamical evolution of the core and, thus, the chemical age derived by our model does not necessarily coincide with the actual dynamical age of the core. In fact, if the initial conditions of the cores were different, the final COM chemical composition of the cores would be significantly altered as compared to our chemical modelling.

For L1521E, Scibelli et al. (2021) also used the monaco code to estimate the age of the core (of ∼6 × 10⁴ yr). However, they could only reproduce the observed abundances of acetaldehyde (CH₃CHO), out of a total of five modelled COMs (with differences in the predicted abundances by factors of 20–160). Therefore, additional aspects in the chemical modelling of these young starless cores such as different initial conditions and/or dynamical evolution of the cores, may be required.

6 EVOLUTION OF THE HC₃N / CH₃CN ABUNDANCE RATIO ACROSS LOW-MASS STAR FORMATION

Recent observational campaigns have been devoted to investigate the chemical composition in COMs (both O-bearing and N-bearing) towards all stages in the process of star formation. In this section, we compare the information available for several sources for two of these COMs and COM precursors, the N-bearing species cyanoacetylene (HC₃N) and acetonitrile (CH₃CN), with the abundances obtained towards our limited sample of starless cores (L1544, L1498, L1517B, and L1521E). We have selected these two molecular species because they have been measured systematically in the past few years across a variety of objects, from Class 0/I protostars (Bergner et al. 2017) to Class II protoplanetary discs (see the maps ALMA large program; e.g., Ilee et al. 2021), and comets (Biver et al. 2021).

Table 5.

Observed and modelled abundances (χ_obs, χ_mod) for several COMs and COM precursors in L1517B.

	Dust peak			Methanol peak
Species	|${\chi }_{\mathrm{\mathbf {obs}}}$|	|${\chi }_{\mathrm{\mathbf {mod}}}$|	Agreement	|${\chi }_{\mathrm{\mathbf {obs}}}$|	|${\chi }_{\mathrm{\mathbf {mod}}}$|	Agreement
CH3OH	|$2.7_{-0.5}^{+0.7}\times 10^{-10}$|	1.21 × 10−9	+	|$1.11_{-0.10}^{+0.12}\times 10^{-9}$|	2.23 × 10−9	+
CH3O	|$8.0_{-2.5}^{+2.9}\times 10^{-12}$|	3.0 × 10−12	+	|$1.9_{-0.6}^{+0.7}\times 10^{-11}$|	6.6 × 10−12	+
CCCO	|$\lt \, 8\times 10^{-12}$|	1.12 × 10−11	+	|$\lt \, 1.3\times 10^{-10}$|	2.05 × 10−11	+
t-HCOOH	|$\lt \, 1.0\times 10^{-10}$|	8.6 × 10−11	+	|$\lt \, 2.0\times 10^{-10}$|	1.63 × 10−10	+
H2CCO	|$2.4_{-0.5}^{+0.6}\times 10^{-11}$|	1.12 × 10−10	+	|$7.6_{-2.0}^{+2.3}\times 10^{-11}$|	2.17 × 10−10	+
CH3OCH3	|$\lt \, 3\times 10^{-11}$|	2.32 × 10−13	+	|$\lt \, 5\times 10^{-11}$|	5.2 × 10−13	+
CH3CHO	|$\lt \, 8\times 10^{-12}$|	2.09 × 10−12	+	|$2.2_{-0.4}^{+0.5}\times 10^{-11}$|	4.0 × 10−12	+
HCCCHO	|$\lt \, 5\times 10^{-12}$|	1.62 × 10−13	+	|$\lt \, 2.1\times 10^{-11}$|	3.2 × 10−13	+
CH3CN	|$6.0_{-1.9}^{+2.1}\times 10^{-12}$|	4.0 × 10−13	–	|$\lt \, 3\times 10^{-11}$|	7.0 × 10−13	+
CH2CHCN	|$\lt \, 1.2\times 10^{-12}$|	1.46 × 10−14	+	|$\lt \, 3\times 10^{-12}$|	2.9 × 10−14	+
HCCCN	|$1.48_{-0.18}^{+0.24}\times 10^{-10}$|	1.04 × 10−11	–	|$3.9_{-0.4}^{+0.4}\times 10^{-10}$|	1.85 × 10−11	–

Note. The modelled and measured abundances agree (+) or not (−) within a factor of 10.

Fig. 9 reports the values of the column density ratio of HC₃N/CH₃CN for the starless cores L1544, L1498, and L1517B (for both positions measured in these cores), the starless core L1521E in the dust peak, 12 Class 0/I protostellar systems (studied by Bergner et al. 2017), 4 Class II protoplanetary discs (studied by Ilee et al. 2021), and 2 comets (studied by Biver et al. 2021, and references therein). To obtain this ratio, we used the column densities for both molecules, except for the comets, for which we employed the abundances of both molecules with respect to water. From Fig. 9, it seems that there is a trend for the starless cores to present much higher HC₃N/CH₃CN ratios (of |$\sim \,$|200–700) than those measured towards Class 0/I protostars and Class II protoplanetary discs, and towards comets (in this case, their HC₃N/CH₃CN abundance ratios reach values even lower than 1). The core L1521E would be an exception, as its ratio is similar to those of the protostars. In any case, the observed general trend suggests that the observed HC₃N/CH₃CN abundance ratio decreases when starless cores enter into the Class 0 phase, remaining approximately constant until the Class II protoplanetary phase, and droping by two orders of magnitude for comets.

Figure 9.

Abundance ratio between cyanoacetylene (HC₃N) and acetonitrile (CH₃CN) for different sources: the starless cores L1517B, L1498, and L1544 (with both positions, the dust and methanol peaks), the starless core L1521E in the dust peak, 12 Class 0/I protostars, 4 Class II proto-planetary discs, and 2 comets (46P and 67P). The references for each source are: L1517B, this work; L1498, Jiménez-Serra et al. (2021); L1544, Jiménez-Serra et al. (2016); L1521E, Nagy et al. (2019); Class 0/I protostellar systems, Bergner et al. (2017); Class II protoplanetary discs, Ilee et al. (2021); comets, Biver et al. (2021). The points of the plot as well as the original data used to obtain the ratios can be found in Appendix  C.

The observed range of values for the HC₃N/CH₃CN ratio could be explained by the fact that it is difficult to form HC₃N at late evolutionary stages in the process of low-mass star formation. Indeed, HC₃N is mainly a gas-phase product whose formation is favoured by the low-densities and cold temperatures found in the outer envelopes of pre-stellar systems and protostars (see section 5.1 in Bergner et al. 2017 for more details). In the same work, the authors report similar CH₃CN abundances for different kinds of sources (young stellar objects, comets, hot cores, and hot corinos), which would explain the decreasing trend in the HC₃N/CH₃CN abundance ratio with time.

7 CONCLUSIONS

The high-sensitivity spectra obtained for the starless core L1517B reveal a chemical complexity a bit poorer than that observed in the core L1498, and quite lower than in L1544, as it only presents a few and simple O-bearing and N-bearing species, such as CH₃O, H₂CCO, CH₃CN, and HCCCN, and in general with lower abundances than in L1498 and L1544. Quantitatively, the geometric mean of the abundance of the targeted species in L1517B is about half than in L1498 and L1544. Similarly to these latter cores, the targeted molecules are more abundant by a factor of ∼3 towards the methanol peak of L1517B, located at ∼5000 au from the core’s centre.

We have also modelled the chemical evolution of L1517B and, except for HCCCN and CH₃CN, our simulations agree with the observed abundances and upper limits within a factor of 10. The model also predicts the observed enhancement of the abundances as we move further away from the core’s centre, although we overestimate the distance of the methanol peak by a factor of 3. Our model suggests that L1517B is chemically older than L1498, but the absence of infall motions and its relatively poor chemical complexity makes us think that this core is at an earlier evolutionary stage than L1498. Actually, our model would predict a younger age for L1517B than L1498 if the central density of the core were higher.

We propose a scenario in which N-bearing molecules are formed first, followed by O-bearing molecules once the catastrophic depletion of CO takes place. Furthermore, complexity increases with time, with bigger molecules being formed at later stages of the core, as observed in L1544. We note, however, that the starless core L1521E does not seem to follow this trend, although the sample of molecules studied is smaller than for the other three cores and the beam of the ARO observations is larger. More observational studies of starless and pre-stellar cores are needed in order to clearly determine the influence of the environment, initial conditions and dynamics on the chemical evolution of the cores.

Finally, we have also studied the HC₃N/CH₃CN abundance ratio measured in sources at different stages in the formation of a low-mass stellar system, observing a decreasing trend over time, which could be explained by the adverse conditions of the late stages of the low-mass star formation to form HC₃N.

ACKNOWLEDGEMENTS

We thank an anonymous referee for their valuable comments and suggestions, which improved the quality of the paper. We also thank J. Aguirre and M. Fernández-Ruz for their input about ways to measure complexity. AM aknowledge support from grant PRE2019-091471 under project MDM-2017-0737-19-2 (Unidad de Excelencia ‘María de Maeztu’ – Centro de Astrobiología) funded by the Spanish Ministry of Science and Innovation/State Agency of Research, MCIN/AEI/10.13039/501100011033, and by ‘ERDF, A way of making Europe’. IJS and JMP aknowledge support from grant PID2019-105552RB-C41 by the Spanish Ministry of Science and Innovation/State Agency of Research, MCIN/AEI/10.13039/501100011033, and by ‘ERDF, A way of making Europe’. The work by AV is supported by the Russian Ministry of Science and Higher Education via the State Assignment contract FEUZ-2020-0038. The Instituto de Radioastronomía Milimétrica (IRAM) is supported by the Institut National des Sciences de l'Univers (INSU) of the Centre National de la Recherche Scientifique (CNRS, France), the Max-Planck-Gesellschaft (MPG, Germany), and the Instituto Geográfico Nacional (IGN, Spain).

DATA AVAILABILITY

The original spectra used in this work will be shared on reasonable request to the corresponding author. Most of the codes used to process and visualize the data can be downloaded from GitHub.⁸

Footnotes

1

https://publicwiki.iram.es/Iram30mEfficiencies

2

https://www.iram.fr/IRAMFR/GILDAS

3

https://www.python.org

4

https://github.com/andresmegias/gildas-class-processing

5

The value used for the uncertainty is 1/30 ≃ 0.03, which is the difference in the fractional similarity produced by one entry that is different in the two vectors.

6

We used logarithmic weights due to the range of the values of our abundances, which encloses several orders of magnitude.

7

We used the geometric mean due to the range of the values of our abundances.

8

https://github.com/andresmegias/gildas-class-python

https://github.com/andresmegias/madcuba-slim-scripts

https://github.com/andresmegias/radex-python

https://github.com/andresmegias/richvalues

9

https://www.python.org

10

https://scipy.org

11

https://github.com/andresmegias/radex-python

12

We prefer to use the median, as it is more robust to outliers.

13

https://pypi.org/project/richvalues

APPENDIX A: DERIVATION OF THE METHANOL COLUMN DENSITY AND ITS UNCERTAINTY WITH radex

In this section, we explain the procedure we have followed to determine the column densities of methanol, acetonitrile and cyanoacetilene, and their uncertainties, using the non-LTE molecular excitation code, radex (Van der Tak et al. 2007). It is a radiative transfer model that allows us to predict the line intensity of the selected molecule from the following parameters:

• Molecule. It has to be included in the radex molecule list, therefore having collisional cross-section derived, experimentally or theoretically.

• Spectral range. Frequency range for the transitions to be predicted for the selected molecule.

• Excitation conditions.

  • Molecular hydrogen number density, |$n_\mathrm{H_2}$|⁠.

  • Kinetic temperature of the molecule, T_kin.

  • Background temperature, T_bg.

• Radiative transfer parameters.

  • Column density, N.

  • Line width, Δv.

For each species, the results of the calculation would be a set of intensity lines, |$\lbrace \tilde{I_i}\rbrace$|⁠, where i is the index of the line. In our case, the molecules are methanol (CH₃OH) A and E, acetonitrile (CH₃CN), and cyanoacetilene (HC₃N), and we have observations for their lines in a specific spectral range (see Table 2), so we have a set of observed intensities, {I_i}. We have an uncertainty for each line (the RMS noise), so we will also have a set of uncertainties {ΔI_i}. We use the frequency range of the observed lines, plus a little margin, for the radex calculations. The rest of the parameters except the column density (and optionally the H₂ number density) are known (see Section 3.2 for more details):

• Molecular hydrogen number density (not fixed for the case of methanol E and cyanoacetilene):

  • Core’s centre: |$n_\mathrm{H_2}$| = 2.20 × 10⁵ cm⁻³ (Tafalla et al. 2004).

  • Methanol peak: |$n_\mathrm{H_2}$| = 1.24 × 10⁵ cm⁻³ (Tafalla et al. 2004).

• Background temperature: T_bg = 2.73 K (temperature of the cosmic microwave background).

• Kinetic temperature: T_kin = 10 K (Tafalla et al. 2004).

• Line width: Δv ∼ 0.3 km s⁻¹ (we choose for each case the result of the madcuba fit).

Our goal is to optimize the column density and optionally also the H₂ number density (when more than one transition is available) to minimize the difference between the observed lines, {I_i}, and the predicted ones, |$\lbrace \hat{I_i}\rbrace$|⁠, taking into account the uncertainties, {ΔI_i}. To measure that difference, we define a loss function, |$\mathcal {L}$|⁠, which we choose to be the chi-square (⁠|$\chi^2$|⁠) error:

$$\begin{eqnarray} \mathcal {L} = \sum _i {\left(\frac{\hat{I}_i - I_i}{\Delta I_i} \right) ^2} . \end{eqnarray}$$

(A1)

In this way, we are measuring the quadratic error weighted with the squared inverse of the uncertainty, so that a difference between the observations and the model contributes less to the loss if the uncertainty on the line intensity is greater.

Once we have properly defined our loss function, we can minimize it with respect to the column density for each of our 8 cases (four species, methanol A and E, cyanoacetilene, and acetonitrile; and two positions, dust and methanol peaks), and also with respect to the H₂ number density for methanol E and cyanoacetilene. We did this with Python,⁹ using the methods cobyla and Nelder–Mead within the function minimize from the library SciPy.¹⁰

Finally, we have to estimate the uncertainty for the optimized value/s of our parameter/s. To do so, we take a threshold for the loss, greater than the minimized value, which defines a lower and an upper uncertainty. In order to find a proper threshold, we create 10⁵ sets of variations of the observed intensities, {{I_i}_j}, so that the value of each intensity, I_i,j, comes from a normal distribution with mean I_i and standard deviation ΔI_i. Then, we calculate the loss for the optimized model with equation (A1) for each set of intensities {I_i}_j, obtaining a distribution of losses. We choose the threshold so that the loss values lower than such a threshold are the 68.27 percent of the total values (similarly to the definition of |$1\, \sigma$| confidence interval, but for an asymmetric distribution with a minimum value). The threshold thus corresponds to the 68.27 percentile (see Fig. A1).

Figure A1.

Histogram of the distribution of losses obtained from the intensities of the observed lines and their uncertainties, for the case of methanol E in the dust peak of L1517B. The scale in the horizontal axis is symmetrical-logarithmic with threshold 1, that is, it is linear from 0 to 1 and logarithmic beyond 1. The blue line on the left shows the value of the loss without taking into account the uncertainties of the observed lines (equation A1), which is almost 0; and the blue dashed line indicates the 68.27 percentile.

Then, we have to explore our parameter space around the minimized values in order to find when the loss function reaches our threshold value. If we only fit the column density, we only have to calculate the loss function in the surroundings of the minimized value until we reach the threshold value. However, if we also want to minimize the H₂ number density, we have to explore a bidimensional parameter space. We opted to make an adaptive grid that starts calculating the loss values at the surroundings of the obtained minimum and continues enlarging its size until it encloses all the loss values minor to the threshold with a reasonable margin. This way, we obtained Figs A2 and A3.

Figure A2.

Plot of the loss (⁠|$\mathcal {L}$|⁠) versus the column density (N) and the H₂ number density (⁠|$n_\mathrm{H_2}$|⁠) for a range of values around the optimized values that minimize the loss function, for the case of methanol E in the dust peak of L1517B. The gray lines indicate the optimized values for N and |$n_\mathrm{H_2}$|⁠. The colour map is composed by two ranges of colours in order to highlight the points in which the loss (⁠|$\mathcal {L}$|⁠) is lower than the threshold value (marked in blue in the colour bar), which will define the uncertainties for both optimized parameters.

Figure A3.

Same as Fig. A2, but for cyanoacetilene (HC₃N) for the dust peak of L1517B.

Following this procedure, we have written a Python script that uses radex and analyses our data in an automated way.¹¹ In this way, we obtained the column density values for methanol, acetonitrile, and cyanoacetilene in L1517B shown in Table 3.

APPENDIX B: PROPAGATION OF UNCERTAINTIES AND TREATMENT OF UPPER LIMITS

Let’s consider a group of n variables {x_i}, with central values μ_i and uncertainties σ_i. This means that the probability density function (PDF) of the variable x_i is centred around μ_i with a width of the order of σ_i, so that the |$1\, \sigma$| confidence interval (68.27 percent) is (μ_i − σ_i, μ_i + σ_i). To propagate the uncertainties through a function f applied to the variables {x_i}, we can draw a sample of a large number of values |$(\gtrsim \, 10^4\sqrt{n})$| of each variable x_i, apply the function to each of the elements of the samples, and then obtain a central value and an uncertainty for the resulting distribution.

To do so, we need two things: an appropriate PDF for converting each variable x_i to a distribution of values, and a proper method to obtain a central value and an uncertainty from the resulting distribution. For the last task, we can use the mean or the median as the central value,¹² and the |$1\, \sigma$| confidence interval with respect to the central value (that includes 68.27 per cent of the distribution) to obtain the lower and upper uncertainties. As for the PDF, if the domain of the variable x is (−∞, ∞), a proper function is the normal distribution:

$$\begin{eqnarray} f(x) = \frac{1}{\sqrt{\tau} \,\, \sigma } \exp \left(- \frac{1}{2} \left(\frac{x-\mu }{\sigma } \right)^2 \right) \,\,\,\, , \end{eqnarray}$$

(B1)

with |$\tau\equiv 2\pi$|⁠. However, if the domain of the variable is not (−∞, ∞), this function would be incorrect. Therefore, we have to use another function as the PDF.

Let’s suppose a domain (b₁, b₂). We define the left and right amplitudes, a₁ and a₂, as the distances between the limits of the domain and the median, that is, a_j = |b_j − μ| for j = 1, 2. Now, as these amplitudes can be different, we will split our desired PDF in two halfs, one for x ≤ μ and other for x > μ. Then, we will use an amplitude a as a reference, which must be greater than the uncertainty, a > σ. If the amplitude is quite greater than the uncertainty, a ≫ σ, a good PDF would be just the normal distribution truncated to the domain (b₁, b₂). However, for amplitudes closer to the uncertainty, it would be considerably incorrect, as the truncation shifts the median of the distribution and modifies the confidence intervals, and thus the uncertainties. To fix this, we make the following variable change:

$$\begin{eqnarray} x - \mu \,\,\,\, \rightarrow \,\,\,\, \tilde{x} - \mu \,\, \equiv \,\, \frac{4}{\tau} \, a \, \tan \left(\frac{\tau}{4} \,\, \frac{x - \mu }{a} \right) . \end{eqnarray}$$

(B2)

Using this new variable |$\tilde{x}$| with a normal distribution, we are able to compress the original domain of (−∞, ∞) to (−a, a). However, we get two disadvantages: first, the normalization constant is now different, and secondly, the relation between the parameter of the standard deviation of the original normal distribution and the |$1\, \sigma$| confidence interval (from which we define the uncertainty, σ) is now different; therefore, we should rename the standard deviation of the original normal distribution to s, which we will call width. The first change is not a problem, as the PDF will be normalized computationally for each case. As for the second one, we have characterized computationally the relation between σ and s. It happens that s/σ decreases with a/σ, having s/σ ≃ 2.65 when a/σ = 2 and s/σ → 1 when a/σ → ∞. We have then a relation between the width, the uncertainty and the amplitude, that is, s = s(σ, a).

Therefore, the resulting PDF would be the following:

$$\begin{eqnarray} f(x) \,\, \propto \,\, \exp \left(- \frac{1}{2} \left(\frac{\frac{4}{\tau} a \tan \left(\frac{\tau}{4} \frac{x - \mu }{a} \right)}{s(\sigma ,a)} \right)^2 \right). \end{eqnarray}$$

(B3)

We observed that s/σ increases quasi-exponentially as a/σ approaches a minimum value of |$\sim \,$|1.47. Moreover, we discovered that if a/σ ≲ 1.7, the estimated relation s(σ, a) starts to be incorrect due of the increasing dispersion in s. Therefore, we must use another PDF for the cases in which a/σ ≲ 1.7. Actually, for that limit case the shape of the PDF is almost a uniform distribution between −a and a. We can model this PDF as a trapezoidal function with a rectangular core and triangles in the edges. By doing so, one can easily demonstrate that this kind of PDF would only work if a/σ is greater than the inverse of the integrated area corresponding to the |$1\, \sigma$| confidence interval (⁠|$\sim \,$|0.683); that is, a/σ ≳ 1.46, which is consistent with the found asymptote of the calculated relationship s(σ, a).

We then choose to only use this distribution with a ≥ 2σ, for the sake of simplicity and following a conservative approach. For a ≫ σ, this PDF tends to a normal PDF, so we use this general PDF instead of a truncated Gaussian for this regime (see Fig. B1). We call this distribution a bounded normal distribution.

Figure B1.

Probability density functions for a normal distribution and two bounded normal distributions with different domains.

Finally, for σ < a < 2σ, we use a shifted and mirrored lognormal distribution whose PDF would be:

$$\begin{eqnarray} f(x) &=& \frac{1}{\sqrt{\tau} \,\, s^{\prime }} \frac{1}{|x-(\mu -a)|}\\ &&\times \,\,\exp \left(- \frac{1}{2} \left(\frac{|x-(\mu -a)|-m^{\prime }}{s^{\prime }} \right)^2 \right), \end{eqnarray}$$

(B4)

with m′ = ln (a) and s′ = ln (a − σ) − ln (a). This PDF meets our requirements: the median is equal to μ and the uncertainty associated with the |$1\, \sigma$| confidence interval is σ. We call this distribution a mirrored lognormal distribution.

We have defined our PDFs for the case of a variable x with a central value μ and an uncertainty σ, building the final PDF with two halfs with amplitude a = a_j for j = 1, 2. In case we had lower and upper uncertainties, σ₁ and σ₂, we should just replace σ by σ₁ for the left half of the PDF and by σ₂ for the right half.

Lastly, we should also address the case of a variable with an upper/lower limit or even a finite interval. Let’s consider an interval (x₁, x₂), which may be finite or infinite. If it is finite, we choose a uniform distribution between x₁ and x₂ as the corresponding PDF. But if it is infinite, we choose a log-uniform distribution with finite thresholds for 0 and ±∞, which we set to ±10⁻⁹⁰ and ±10⁹⁰. For example, for an interval of (−100, ∞), we would build a sample {x₋} from a uniform distribution between −90 and 2 and a sample {x₊} from a uniform distribution between −90 and 90. Our final distribution would be the joining of the samples of |$\lbrace -10^{\lbrace x_{-}\rbrace }\rbrace$| and |$\lbrace 10^{\lbrace x_{+}\rbrace }\rbrace$|⁠.

To sum up, if we have a set of variables with central values μ_i and uncertainties σ₁, σ₂, we first build distributions {x_i} using the mentioned PDFs. Then, we apply the function to the distributions, f({x_i}), obtaining a new distribution. Finally, we use an algorithm to detect if the distribution corresponds to an interval (that can be an upper/lower limit) or a defined value with uncertainties, and derive the corresponding parameters.

Using this approach, we have written a Python library that allows to deal with values with uncertainties and upper/lower limits, performing the uncertainty propagation automatically.¹³ We have used it throughout the calculations of this paper, e.g. for obtaining the values shown in Fig. 9.

APPENDIX C: ABUNDANCE RATIO OF CIANOACETYLENE AND ACETONITRILE

Table C1 shows the abundance ratio of cyanoacetylene (HC₃N) and acetonitrile (CH₃CN) through different sources, from starless cores to comets. The propagation of the uncertainties and the upper limits in the column density is done through simulations using statistical distributions (see Appendix  B).

Table C1.

Column densities of HC₃N and CH₃CN (N) and abundance ratio of this molecules (⁠|$\chi _{\mathrm{HC_3N}} \, / \, \chi _{\mathrm{CH_{3}CN}}$|⁠) for the starless cores L1517B, L1498, L1544, and L1521E, 12 protostellar systems, 4 protoplanetary discs, and 2 comets.

Source			|${N_{\mathrm{HC_{3}N}}}\, (\mathrm{cm}{}^{-2})$|	|${N_{\mathrm{CH_{3}CN}}}\, (\mathrm{cm}{}^{-2})$|	|${\chi _{\mathrm{HC_{3}N}}\, /\, \chi _{\mathrm{CH_{3}CN}}}$|
Starless cores	L1517B	Dust peak	|$(5.16\, \pm \, 0.05)\times 10^{13}$|	|$(2.1\, \pm \, 0.6)\times 10^{11}$|	|$250_{-50}^{+100}$|
		Methanol peak	|$3.72_{-0.13}^{+0.14}\times 10^{12}$|	|$\lt \, 2.0\times 10^{11}$|	|$\gt \, 17$|
	L1498	Dust peak	|$(1.6\, \pm \, 0.3)\times 10^{13}$|	|$\lt \, 8\times 10^{10}$|	|$\gt \, 100$|
		Methanol peak	|$(2.2\, \pm \, 1.0)\times 10^{13}$|	|$(1.0\, \pm \, 0.1)\times 10^{11}$|	|$220_{-100}^{+100}$|
	L1544	Dust peak	|$(1.0\, \pm \, 0.3)\times 10^{14}$|	|$(1.5\, \pm \, 0.2)\times 10^{11}$|	|$670_{-200}^{+230}$|
		Methanol peak	|$(4.2\, \pm \, 0.4)\times 10^{13}$|	|$\lt \, 9.1\times 10^{10}$|	|$\gt \, 310$|
	L1521E	Dust peak	|$(8.4\, \pm \, 2.5)\times 10^{12}$|	|$(4.8\, \pm \, 1.4)\times 10^{11}$|	|$18_{-6}^{+9}$|
Protostellar systems	B1-a		|$(4.2\, \pm \, 1.2)\times 10^{12}$|	|$(4.9\, \pm \, 1.1)\times 10^{11}$|	|$9_{-3}^{+3}$|
	B1-c		|$(4\, \pm \, 3)\times 10^{12}$|	|$(3.5\, \pm \, 0.6)\times 10^{11}$|	|$12_{-9}^{+11}$|
	B5 IRS1		|$(3.2\, \pm \, 1.2)\times 10^{12}$|	|$\lt \, 1.7\times 10^{11}$|	|$\gt \, 5$|
	IRAS 03235		|$(3.9\, \pm \, 1.0)\times 10^{12}$|	|$(1.6\, \pm \, 0.9)\times 10^{11}$|	|$25_{-11}^{+31}$|
	IRAS 03245		|$(4\, \pm \, 3)\times 10^{12}$|	|$(1.9\, \pm \, 1.5)\times 10^{11}$|	|$19_{-15}^{+69}$|
	IRAS 03271		|$(1.7\, \pm \, 1.4)\times 10^{12}$|	|$(1.9\, \pm \, 1.1)\times 10^{11}$|	|$9_{-7}^{+14}$|
	IRAS 23238		|$(3.2\, \pm \, 2.6)\times 10^{12}$|	|$(1.4\, \pm \, 0.3)\times 10^{11}$|	|$22_{-18}^{+21}$|
	L1014 IRS		|$(4.4\, \pm \, 3.6)\times 10^{11}$|	|$\lt \, 6\times 10^{10}$|	|$\gt \, 0.7$|
	L1455 IRS3		|$(6\, \pm \, 5)\times 10^{11}$|	|$\lt \, 9\times 10^{10}$|	|$\gt \, 0.7$|
	L1455 SMM1		|$(2.4\, \pm \, 1.9)\times 10^{12}$|	|$\lt \, 1.1\times 10^{11}$|	|$\gt \, 2.2$|
	L1489 IRS		|$(3.5\, \pm \, 2.4)\times 10^{12}$|	|$\lt \, 1.4\times 10^{11}$|	|$\gt \, 0.3$|
	SVS 4–5		|$(1.1\, \pm \, 0.3)\times 10^{13}$|	|$(5.2 \pm 0.9) \times 10^{11}$|	|$21_{-6}^{+8}$|
Protoplanetary discs	GM Aur		|$1.9_{-0.4}^{+0.4}\times 10^{13}$|	|$2.1_{-0.1}^{+0.2}\times 10^{12}$|	|$8.8_{-1.9}^{+2.0}$|
	As 209		|$2.9_{-0.5}^{+0.5}\times 10^{13}$|	|$1.7_{-0.2}^{+0.2}\times 10^{12}$|	|$17_{-3}^{+4}$|
	HD 163296		|$7_{-2}^{+3}\times 10^{13}$|	|$2.3_{-0.2}^{+0.2}\times 10^{12}$|	|$32_{-9}^{+11}$|
	MWC 480		|$8_{-3}^{+4}\times 10^{13}$|	|$3.5_{-0.2}^{+0.2}\times 10^{12}$|	|$22_{-8}^{+11}$|
Source			|$\chi _{\mathrm{HC_{3}N}}\, /\, \chi _{\mathrm{H_{2}O}}$|	|$\chi _{\mathrm{CH_{3}CN}}\, /\, \chi _{\mathrm{H_{2}O}}$|	|$\chi _{\mathrm{HC_{3}N}}\, /\, \chi _{\mathrm{CH_{3}CN}}$|
Comets	46P		|$\lt \, 3\times 10^{-5}$|	|$(1.7\, \pm \, 0.01)\times 10^{-4}$|	|$\lt \, 0.21$|
	67P		4 × 10−6	5.9 × 10−5	0.068

Note. References for each source: L1517B: this work; L1498: Jiménez-Serra et al. (2021); L1544: Jiménez-Serra et al. (2016); L1521E: Nagy et al. (2019); protostellar systems: Bergner et al. (2017); protoplanetary discs: Ilee et al. (2021); comets: Biver et al. (2021).