© The Authors 2023

1 Introduction

Deuterated molecules are important diagnostic tools of the earliest phases of star formation, allowing astronomers to study the central regions of pre-stellar cores in detail (e.g. Caselli et al. 2002b; Pineda et al. 2022), where CO and other CO-bearing molecules are heavily frozen onto dust grains (e.g. Caselli et al. 1999; Bergin & Tafalla 2007). High levels of deuteration are also found in the cold envelopes of newly formed protostars (van Dishoeck et al. 1995; Hatchell et al. 1998; Roberts et al. 2002; Parise et al. 2004; Emprechtinger et al. 2009), bringing up the question of how much of this material can survive during the star and planet formation process. High levels of deuteration have been found in planet-forming disks (e.g. Mathews et al. 2013), comets (e.g. Altwegg et al. 2015), and carbonaceous chondrites (e.g. Robert 2003; Busemann et al. 2006; Ceccarelli et al. 2014). Our oceans are also enriched in heavy water, partially inherited from the solar pre-stellar phase (Cleeves et al. 2014). It is therefore important to understand the chemical processes that regulate deuterium fractionation in the early phases of star formation in detail and investigate possible differences between various evolutionary stages. In particular, observational constraints are needed for the chemical models to discriminate among the processes of deuteration happening in the gas phase and on the surface of dust grains.

At low temperatures, deuterium fractionation is driven by the gas phase reaction

$${\text{H}}_3^{+}+\text{HD}\rightleftharpoons {\text{H}}_2{\text{D}}^{+}+{\text{H}}_2+232\text{\hspace{0.17em}K}.$$(1)

Due to its exothermicity, the reaction proceeds more efficiently from left to right with decreasing temperature (assuming a low ortho-to-para H₂ ratio; e.g. Pagani et al. 1992). This leads to an enhancement of the abundance of H2D⁺ and an associated increase in the abundances of D₂H⁺ and ${\text{D}}_3^{+}$, and subsequently to efficient deuteration of other molecules. In the cold and dense conditions of the earliest stages of star formation, the main destructor of H₂D⁺, CO, is highly depleted from the gas phase and frozen out onto the surfaces of dust grains (e.g. Caselli et al. 1999), and this further increases the deuterium fractionation (Dalgarno & Lepp 1984). The deuterium fraction of a molecular tracer is measured by dividing the column density of the deuterated molecule by the column density of the hydrogenated molecule.

In this paper, we focus on deuterated molecules present in the well-known pre-stellar core L1544 in Taurus, and in one of the youngest (and most highly deuterated) Class 0 sources, HH211 in Perseus. L1544 is located at a distance of 170 pc (Galli et al. 2019), with clear evidence of gravitational contraction (Caselli et al. 2012). Its central temperature approaches 6 K (Crapsi et al. 2007) and its central density reaches 10⁷ cm⁻³ (Keto & Caselli 2010; Caselli et al. 2019). Towards its centre, L1544 exhibits a high level of deuteration (Crapsi et al. 2005; Redaelli et al. 2019). Within the central 2000 au, recent Atacama Large Millimeter/submillimeter Array observations have been found consistent with almost complete (99.999%) depletions of elements heavier than He (Caselli et al. 2022). L1544 is on the verge of star formation, thus shedding light on the initial conditions in the process of star formation. HH211 is a newly born Class 0 protostar (e.g. Enoch et al. 2006; Lee et al. 2018) at a distance of 321 pc (Ortiz-León et al. 2018), which presents high degrees of deuterium fractionation (Emprechtinger et al. 2009). It hosts a jet-driven molecular outflow (McCaughrean et al. 1994; Gueth & Guil-loteau 1999), and has a Keplerian disk (Segura-Cox et al. 2016). The surrounding envelope is elongated in a direction roughly perpendicular to the jet and outflow axis (Gueth & Guilloteau 1999). In Chantzos et al. (2018), it is shown that the deuteration of cyclopropenylidene, c-C₃H₂, is more efficient towards HH211 (N(c-C₃HD)/N(c-C₃H₂)=20%) than in L1544 (N(c-C₃HD)/N(c-C₃H₂)=10%). Gas phase deuteration processes are sufficient to reproduce the N(c-C₃HD)/N(c-C₃H₂) ratio observed in L1544, so the enhancement of deuteration in HH211 is assumed to be due to efficient deuteration happening on the surface of dust grains in the pre-stellar phase, followed by the release from the grains in the protostellar stage.

We present a survey of ground state-rotational lines of deuterated molecules towards L1544 and HH211, observed with the single-dish 20 m radio telescope at the Onsala Space Observatory. Our observations cover CCH, HCN, HNC, and HCO⁺, and their deuterated isotopologues. To avoid optical depth limitations with the main isotopologues, we also observed ¹³CCH, C¹³CH, H¹³CN, HN¹³C, H¹³CO⁺, and HC¹⁸O⁺ This data set allows for a unique comparison between the deuteration in a very dynamically evolved pre-stellar core on the verge of forming a low-mass star, and a young Class 0 protostellar core, and thus, to study the influence of the evolutionary stage on the deuterium fractionation. Such a comparison is important to understand how the deuterated molecules are inherited in the earlier stages of formation of a low-mass star, and eventually in its planetary system.

In Sect. 2, we describe the observations, followed by the results in Sect. 3. The analysis in Sect. 4 covers the derivation of the column densities and deuterium fraction assuming local thermal equilibrium (LTE) conditions. In Sect. 5, we use the non-LTE radiative transfer code LOC to model our observations. We discuss the results of the LTE and non-LTE analysis in Sect. 6 and present our conclusions in Sect. 7.

2 Observations

The data presented in this paper were obtained with the 20m radio telescope at the Onsala Space Observatory. The observed spectra are centred on the dust emission peaks of the sources (L1544: α₂₀₀₀ = 05^h04^m17^s.21, δ₂₀₀₀ = +25°10′42″.8; HH211: α₂₀₀₀ = 03^h43^m56^s.80, δ₂₀₀₀ = +32°00′50″.0). The first part of the observations was done in May 2020 (project ID: O2019b-04). A follow-up project observed the second part from January to May 2021 (project ID: O2020b-05). For the observations, we used the 3 mm (Belitsky et al. 2015) and the 4mm receiver (Walker et al. 2016) in combination with the OSA spectrometer at 19 kHz resolution, with 625 MHz bandwidth and dual polarisation. We applied the dual beam-switched mode, with the off position shifted by 11.8 arcmin in the azimuth direction. The intensity calibration was done using standard Dicke-switched type calibration, switching between a hot load and a sky direction close to the source. The 1σ pointing accuracy is estimated to be 3 arcsec. The frequency resolution of the data corresponds to a velocity resolution of about 0.07 km s⁻¹ at 84 GHz. The weather conditions during the observations were variable, with the system temperature ranging between 120 and 300 K.

The observed transitions are summarised in Table A.1. The corresponding spectra towards L1544 and HH211 are presented in Figs. A.1 and A.2, respectively. A selection of the observed transitions is given in Table 1 and Fig. 1. The average beam size of the observations is about 50 arcsec, which corresponds to a physical size of 8500 au for L1544 and 16 000 au for HH211.

The data processing was done using the GILDAS software (Pety 2005) and the Python package PYSPECKIT (Ginsburg & Mirocha 2011; Ginsburg et al. 2022). The antenna temperature $T_{\text{A}}^{*}$ was converted to the main beam temperature T_mb using the relation $T_{\text{mb}}=T_{\text{A}}^{*}/B_{\text{eff}}$. The corresponding values for the main beam efficiencies (B_eff) of the 20 m telescope are given in Table 1.

3 Results

Many species observed in this survey show either resolved (CCH, ¹³CCH, C¹³CH, CCD, HCN, H¹³CN, DCN), or unresolved (HNC, HN¹³C, DNC, H¹³CO⁺, DCO⁺ hyperfine structure. Due to the amount of hyperfine lines in the resolved case, we decided to select one hyperfine component per species for the LTE analysis, based on optical depth and signal-to-noise ratio (S/N). The properties of the selected lines are given in Table 1, determined by fitting a Gaussian profile to the lines (see Figs. A.1 and A.2). The corresponding spectra are shown in Fig. 1, where lines observed towards L1544 and HH211 are plotted in black and red, respectively. In the following, we describe the observed spectra and discuss the line selection.

3.1 L1544

The radical CCH and its isotopologues present lines with multiple hyperfine (hf) components. In L1544, we detected all six hf components of CCH, where four out of six likely show a dip due to self-absorption. This is caused by absorption in the outer layers of the core and dilutes the signal from the high-density central regions. We selected the line with the lowest optical depth (N=1–0, J=1/2–1/2, F=1–0), to minimise any issues with optical depth effects. In L1544, we successfully observed the components of the ¹³C isotopologues of CCH: For ¹³CCH we detected three of the eight hyperfine components with a signal-to-noise ratio > 3; C¹³CH shows five of the seven components with S/N > 4. For both, we selected the transition with the highest S/N, (N=1–0, J=3/2–1/2, F₁=2–1, F=5/2–3/2). The deuterated isotopologue, CCD is much brighter, and seven of the nine hyperfine components are detected (S/N > 8). Some of them show self-absorption, thus, we selected the line with the lowest S/N (≈8), N=1–0, J=3/2–1/2, F=5/2–3/2, to avoid optical depth issues. The Gaussian fit parameters of the ¹³CCHnd C¹³CHf components show that the lines are shifted by approximately 0.1 km s⁻¹ (V_LSR =7.3 km s⁻¹) compared to CCH and CCD7.2– 7.25 km s⁻¹), which are located at the typical system velocity observed for L1544 (7.2 kms⁻¹ ). However, this shift is close to the spectral resolution of the data (0.07 kms⁻¹) and might therefore not be significant.

HCN is optically thick, showing a double-peaked line profile in each hyperfine component with a dip reaching down to almost zero level, caused by self-absorption in the outer layers of the core. Due to its high optical depth, the line profile provides information on the kinematics in the outer layers of the core. A blue asymmetry in a self-absorbed peak is a typical signature of infall motions - extended inward velocities in the outer regions of a core cause enhanced self-absorption at redshifted velocities (e.g. Evans 1999). The intensity ratio between the two peaks can give information on the infall rates itself (e.g. Keto et al. 2015), because with increasing infall motions, the redder peak decreases in intensity and evolves to a shoulder of the blue peak (like it is seen with HCO⁺n Fig. 1 and for example in Redaelli et al. 2022). However, the peak of the central HCN hf component shows a red asymmetry. Following the argumentation above, this is caused by an expansion motion of the core in its outer regions, more specifically, in the layers traced by the HCN transition.

It has to be noted that the red asymmetry of HCN has a different origin than the spectra of HC₃N (1–0) towards L1544 observed by Bianchi et al. (2023). In the case of the optically thin HC₃N, the red asymmetry is caused by the molecule tracing the southern region of the core, which is redshifted with respect to the rest velocity of L1544 (see the V_LSR map in Fig. B.1 in Spezzano et al. 2016). By overplotting the central hyperfines of HCN and HC₃N, Fig. 2 displays the large line width of HCN. This additionally shows the optical thickness of the transition, as in L1544 the line width observed for optically thin species is approximately 0.4 km s⁻¹ (e.g. Caselli et al. 2002a). The red asymmetry is not observed for the other two hyperfine components of HCN. They show a rather symmetric double-peak profile, possibly tracing a more static or non–moving layer, which was also observed in CS (2–1) by Tafalla et al. (1998). This shows that the hyperfine components clearly trace different layers, depending on their optical depth. In Sect. 5, we analyse this behaviour by applying non–LTE radiative transfer modelling. Due to its optical thickness, the main species cannot be used to derive the column density and deuterium fraction of HCN, so we focus instead on the ¹³C isotopologue.

The hyperfine structures of both the ¹³C- and the deuterated isotopologue of HCN show some self-absorption. A hyperfine structure fitting, using the HFS method provided by the GILDAS package CLASS, reveals that in both cases the weakest component (F = 0–1) is moderately optically thin, with an optical depth of τ = 0.48(9) and τ = 0.9(2) for H¹³CNnd DCN, respectively (numbers in brackets give the uncertainty on the last digit). Therefore, we selected these lines for the further analysis.

Our observations cannot resolve the hyperfine structure of HNC, but only show blended components. In addition, the lines are expected to be heavily self-absorbed like in the case of HCN, which makes it impossible to reconstruct the components.

Therefore, we included this molecule in the non-LTE analysis in Sect. 5 but do not consider it in the derivation of column densities. The lines of the isotopologues, HN¹³Cnd DNC, are split into multiple hyperfine components, mainly due to the non-zero nuclear spins of the D and N nuclei. For HN¹³C van der Tak et al. (2009) identified four effective hyperfine components, with relative intensities of 1.00, 4.04, 6.63, and 3.66, from low to high frequency. In Fig. 3, we fitted the four effective components of HN¹³C to our observed spectrum by applying a multi-component Gaussian. The amplitudes of the fits result in relative intensities of 1.00, 2.18, 2.41, 4.95 (from high to low velocity components), corresponding to optical depths of 0.3, 0.7, 0.8, 4.5, respectively. Therefore, we conclude that all components of HN¹³C except for the weakest component at the reddest velocity, are optically thick and most probably self-absorbed. Hence, we only used the hf component at 87090.675 MHz to derive the column density of HN¹³C For DNC, van der Tak et al. (2009) identified six effective hyperfine components, reporting relative intensities of 0.33, 2.67, 1.67, 2.33, 0.99, 1.00, from low to high frequency. A CLASS HFS fitting gives out optical depths of 0.24, 1.94, 1.21, 1.69, 0.72, 0.72, from low to high frequency. Therefore, we selected the lowest frequency component of the DNC hyperfine structure for our further analysis.

The J=1–0 transition of HCO⁺ shows strong self-absorption in L1544, with a blue asymmetry caused by infall motions in the outer layers of the core. Therefore, we excluded it in the LTE analysis. A detailed discussion and non-LTE modelling of this line can be found in Redaelli et al. (2022). The ¹³C isotopologue of HCO⁺ also shows asymmetric self-absorption. However, we used the line to derive a lower limit for the integrated intensity and column density of HCO⁺ The spectrum of the ¹⁸O iso-topologue shows a flattened top. With higher spectral resolution (20 kHz) and a smaller telescope beam (31″), a small dip is visible (e.g. Caselli et al. 2002a; Redaelli et al. 2019) that is caused by kinematics and CO depletion in the central regions of the core. Due to the lower spectral resolution and larger telescope beam of our data, however, the emission is diluted. The line does not show any asymmetry, as the abundance of the molecule in the outer envelope is too low to be affected by self-absorption. The deuterated isotopologue consists of three blended hyperfine components that we cannot resolve with our spectral resolution. Therefore, we treated the line as one component.

Fig. 1 Selection of the spectra observed towards the pre-stellar core L1544 (black) and the protostellar core HH211 (red). The dotted vertical lines indicate the respective rest velocities of the cores. For HCN, the F = 2–1 transition is shown. For CCH, C13CH CCD H13CN and DCN, the transitions listed in Table 1 are shown.

Fig. 2 Central hyperfine (F=2–1) of the (1–0) transitions of HCN (black) and HC3N (orange), observed towards L1544. The data for HC3N are taken from Bianchi et al. (2023).

Fig. 3 Multi-component Gaussian fit of the four effective hyperfine components found by van der Tak et al. (2009) to the spectrum of HN13Cbserved towards L1544.

3.2 HH211

The spectra observed towards HH211 show an additional velocity component centred at around 8.22 km s⁻¹, which is significantly blue-shifted with respect to the main emission peak (~9.2km s⁻¹). This emission likely originates from the outskirts of a lower velocity plateau which might be part of or influenced by the outflow. It is located to the south-east of the core (Pineda & Friesen, in prep.), which is partly covered by the telescope beam (40–60″, depending on the frequency of the transition) of the Onsala 20 m telescope. To account for the extra emission of the additional velocity component, we apply a two-component Gaussian fit to the spectral lines observed towards HH211, wherever possible. In the analysis, we used the fit parameters of the main component emission to derive the column density and deuterium fraction of the molecular tracers. The only exception are the two ¹³C isotopologues of CCH, where the signal-to-noise ratio is too low to resolve the second component. In this case, a one-component Gaussian fit is sufficient to obtain the line properties. In Table B.1, we list the best-fit parameters and the derived column densities for the additional velocity component. There is no clear trend visible in the FWHM of the additional component compared to the main component. For CCH, CCDnd HCO⁺ the second component appears to be broader than the main component, for the rest it is narrower. However, the line width has a rather large error and in some cases, the two line widths are, within the error bars, the same.

The main emission components of the observed lines all show a Gaussian line shape and do not exhibit any dips or asymmetries caused by self-absorption. This is possibly caused by a combination of infall motions and turbulent motions induced by the protostellar outflow in the less dense regions of the core, which is enlarging the velocity range in the core and thus, decreasing the possibility for self-absorption. In addition, the lines are also tracing outflow material, for example HCO⁺r HNC (e.g. Arce & Sargent 2004). In general, the increased level of turbulence in the core causes a broadening of the line widths (on average 0.2 km s⁻¹ larger than in L1544).

For CCH, we detected all six hf components, and selected the transition with the lowest optical depth to calculate the column density (N=1–0, J = 3/2–1/2, F = 1/1). Towards the proto-stellar core, the ¹³C isotopologues of CCH are less bright than towards L1544, and we detected only two hyperfine components with S/N > 3 for each of them. We chose the components with the higher S/N, N=1–0, J=3/2-1/2, F₁=1–0, F=3/2–1/2 for ¹³CCH and N=1–0, J=3/2–1/2, F₁=2–1, F=5/2–3/2 for C¹³CH For the deuterated isotopologue of CCH, we detected seven out of nine hyperfine transitions. Among them, the weaker lines do not show the additional velocity component. We selected the transition with the highest signal-to-noise ratio (S /N > 13), which is N=1–0, J=3/2–1/2, F=5/2–3/2. The lines of the main isotopologue, CCH, seem to be shifted to lower velocity (~9.0km s⁻¹), while the lines of the rarer isotopologues are located around the rest velocity of the system (~9.1 km s⁻¹).

For HCN, H¹³CN and DCN, we applied the hyperfine-fitting method of CLASS to check the optical depth. As in L1544, all components of HCN are optically thick. We therefore excluded the molecule from the LTE analysis, and only considered it in the radiative transfer modelling in Sect. 5. The optical depths of DCN are similar to those towards L1544, only the weakest hyperfine is optically thin, with τ = 0.35. For H¹³CN however, all hyperfine components turn out to be optically thin, with τ = 0.08–0.4. Therefore, we selected the brightest hyperfine (F=2–1) for the further analysis.

In the spectra of HNC and isotopologues, the hyperfine components are completely blended. There is no clear evidence of self-absorption as it was seen in L1544. However, the main species is expected to be optically thick, thus, we focus on the ¹³C isotopologue. Towards HH211, the hf components of HN¹³Cnd DNC are less distinct than towards L1544, and the reddest component is less pronounced. Due to this and the additional velocity component at lower velocity (~8.2 km s⁻¹ ), we were not able to fit the four effective components of HN¹³C and the six effective components of DNC reported by van der Tak et al. (2009). Therefore, we applied a three-component Gaussian fit to the observed spectrum to account for the additional velocity component and the high velocity hf component separately from the remaining, optically thick emission.

The observed emission of HCO⁺ shows a very pronounced additional velocity component, with almost half of the intensity of the main emission. To avoid any optical depth issues, we did not use HCO⁺ for the LTE analysis. The main velocity components of the ¹³C- and D-bearing isotopologues have an optical depth of 0.9, and therefore have to be interpreted with caution.

4 Analysis

4.1 Excitation temperature and column density

In this work, all column densities were calculated under the assumption of optically thin emission, using the formula presented in Mangum & Shirley (2015). Furthermore, we applied the approximation of a constant excitation temperature throughout the core (CTex), following Caselli et al. (2002b) and Redaelli et al. (2019):

$$N_{\text{tot}}^{\text{thin}}=\frac{8\pi v^2}{c^3}\frac{Q_{\text{rot}}\left(T_{\text{ex}}\right)}{g_u{A}_{\text{ul}}}{\left[J_v\left(T_{\text{ex}}\right)-J_v\left(T_{\text{bg}}\right)\right]}^{-1}\frac{{\text{e}}^{\frac{E_u}{k{T}_{\text{ex}}}}}{{\text{e}}^{\frac{hv}{k{T}_{\text{ex}}}}-1}{\displaystyle \int T_{\text{mb}}\text{d}\upsilon },$$(2)

where Q_rot(T_ex) is the partition function of the molecule at an excitation temperature T_ex, ɡ_u and E_u are the degeneracy and energy of the upper level of the transition, respectively, A_ul the Einstein coefficient for spontaneous emission, T_bg = 2.73 K the temperature of the cosmic microwave background, J(T) the Rayleigh-Jeans equivalent temperature and T_mb the main beam temperature. Assuming a Gaussian line profile, the integrated main beam temperature is calculated by $\frac{1}{2}\sqrt{\pi /\text{In\hspace{0.17em}}2}\cdot \text{Δ}V\cdot T_{\text{mb},\text{peak}}$. The corresponding parameters for each transition are listed in Table 1 and Table 2. We derived the partition function from the energies and degeneracies of the rotational levels (see Mangum & Shirley 2015), taken from the Cologne Database for Molecular Spectroscopy (CDMS Müller et al. 2001).

However, as some of our transitions are only moderately optically thin, we applied the optical depth correction factor (Goldsmith & Langer 1999) to all of our lines to derive the total column density:

$$N_{\text{tot}}=N_{\text{tot}}^{\text{thin}}\frac{\tau }{1-{\text{e}}^{-\tau }}.$$(3)

The optical depth is derived using:

$${\tau }_v=-\text{In\hspace{0.17em}}\left[1-\frac{T_{\text{mb}}}{f\left[J\left(T_{\text{ex}}\right)-J\left(T_{\text{bg}}\right)\right]}\right],$$(4)

where we assumed a filling factor f = 1.

For CCH and isotopologues, we assumed an excitation temperature of 6 K in both sources. This temperature has been applied to similar carbon chain molecules such as c-C₃H₂ in Giers et al. (2022) and Chantzos et al. (2018), following Crapsi et al. (2005) and Emprechtinger et al. (2009). It is also consistent with Taniguchi et al. (2019) who derived T_ex ≈ 6 K for CCH in the starless cores L1521B and L134N. To estimate the excitation temperature of H¹³CN DCN, HN¹³Cnd DNC, we used the results of the hyperfine-fitting done earlier with the HFS method in CLASS. The derived values are listed in Table 2 and are consistent with excitation temperatures used by Padovani et al. (2011), Hily-Blant et al. (2013) and Quénard et al. (2017). For HC¹⁸O⁺ and DCO⁺ we used the excitation temperatures derived by Redaelli et al. (2019) using radiative transfer. For H¹³CO⁺ we assumed an excitation temperature of 6 K, which is an intermediate value and similar to what we used for CCH and isotopologues. In general, we assumed the same or similar excitation temperatures for both the pre-stellar and the protostellar core, as our observations mainly trace the cold outer envelopes of the two cores, where similar conditions hold. Indeed, the effect of the excitation temperature on the derived column density was found to be moderately small. Changing T_ex by one Kelvin introduces an error of about 20%.

The derived column densities of all species observed towards the two sources are presented in Table 3. To analyse the nature of the additional velocity component in HH211, we derived the column density of it separately. The results are presented in the last column of Table B.1. The uncertainties for the column densities were derived by propagating the 1σ errors on the integrated intensity and the optical depth, and adding an additional error of 15% to account for the uncertainties in the flux calibration.

For CCH, the derived column densities are consistent with previous measurements within a factor of two (Sakai et al. 2008; Zhang et al. 2021). The ¹³C isotopologues of CCH are very weak and difficult to observe. Therefore, they have only been detected towards the sources L1521B, TMC-1, L134N/L183, L483 and L1527 so far (Taniguchi et al. 2019; Sakai et al. 2010; Agúndez et al. 2019; Yoshida et al. 2019).

HCN, HNC and isotopologues have been observed multiple times towards L1544 (e.g. Hirota et al. 2003; Hily-Blant et al. 2010; Quénard et al. 2017; Spezzano et al. 2022) and HH211 (Imai et al. 2018) in previous studies, reporting results which are consistent with our measurements. Previous measurements towards HH211 by Roberts et al. (2002) reported column densities for H¹³CN and DCN that are smaller by a factor of two compared to our results. These differences are likely due to different approximations used in the derivation of the column densities, as well as a different telescope beam.

By assuming the carbon isotopic ratio for the local interstellar medium, ¹²C/¹³C = 68 (Milam et al. 2005), we derived a HNC/HCN ratio of 0.3(1) and 0.4(1) for L1544 and HH211, respectively. Previous measurements of this ratio towards L1544 and other cores have revealed a value revolving around 1 (L1544: 1.0(5), Quénard et al. 2017; L1498: 0.9(1), L1521E: 1.1(1), TMC-2: 1.1(1), Padovani et al. 2011), suggesting that HNC and HCN have similar abundances. However, as this is derived using the ¹³C isotopologues, it implies that the fractionation of the two molecules is the same. In Colzi et al. (2020) it is shown that the carbon fractionation of HCN and the carbon fractionation of HNC do behave similarly with a ratio close to 68. However, most of the time, the isotopic ratio for HNC is slightly larger, which can change the resulting HNC/HCN ratio by up to 40%. Therefore, using the same ¹²C/¹³C ratio for the two molecules might not be appropriate.

In L1544, HCO⁺ and its deuterated and non-deuterated isotopologues have been the subject of multiple studies (e.g. Caselli et al. 2002b; Jørgensen et al. 2004; Vastel et al. 2006; Redaelli et al. 2019). Recently, Redaelli et al. (2019) analysed high-sensitivity maps of HC¹⁸O⁺nd DCO⁺ where the peak column densities are consistent with our measurements, taking into account that our observations are taken towards the dust peak and are averaged over the whole beam. To our knowledge, there is no literature data of HCO+nd isotopologues towards HH211 so far.

4.2 Deuterium fraction

The deuterium fractions are derived by dividing the column density of the deuterated isotopologues by the column density of the main species. As the main species of all observed molecules, except CCH, are optically thick, we used the ¹³C and¹⁸ isotopo-logues to estimate their column densities. For this, we multiplied the column densities of the ¹³C isotopologues by ¹²C/¹³C = 68 (Milam et al. 2005). Similarly, the column density of HC¹⁸O⁺ was multiplied by 557 (Wilson 1999) to obtain an estimate for N(HCO⁺. The resulting deuterium fractions are summarised in Table 4. Furthermore, we derived the deuteration level in the additional velocity component of HH211. The results are presented in Table B.2. Despite the high uncertainties in some cases, the additional component and the main emission generally show quite similar deuteration levels. This likely suggests that they both trace the same gas, with similar properties along the line-of-sight.

Colzi et al. (2020) show that in the local interstellar medium the carbon isotopic ratio is molecule-dependent and varies with time, volume density, and temperature, and can deviate significantly from 68. For example, for a fixed density and temperature, it may change by a factor of two between 10⁵ and 10⁶ yr of core evolution, depending on the molecule. In this work, we applied the canonical value, taking note of the uncertainties it might introduce.

When interpreting the deuterium fraction, one has to keep in mind that, first, with single-dish telescopes the observations of deuterated and non-deuterated species correspond to a different beam size and therefore can cover different regions of the core. Inourcase, the beamsizes forthe deuterated isotopologues on average are generally 10″ larger than the telescope beams of the normal and ¹³C isotopologues (see Table A.1). This might dilute the real deuterium fraction of a molecule and cause systematic errors in the ratios. Moreover, the species themselves are abundant in different regions or shells of the cores. In pre-stellar cores, deuterated species are expected to trace higher-density regions close to the centre of the core, where the deuteration process is most efficient (e.g. Giers et al. 2022). Non-deuterated isotopologues, on the other hand, are additionally also abundant in the outer parts of a core. In L1544, this is observed for example in HCO⁺ CH₃OH, H₂CO, and c-C₃H₂ (Redaelli et al. 2019; Chacón-Tanarro et al. 2019; Giers et al. 2022). Towards protostellar cores, various single dish observations have shown high levels of deuteration (van Dishoeck et al. 1995; Parise et al. 2004, e.g.). However, Persson et al. (2018) show that with interferometric observations the D/H ratio is actually lower in the inner regions, as they are less affected by optical depth effects and beam dilution. These uncertainties are important to keep in mind when interpreting single dish observations of the inner warm regions of protostellar cores. However, due to the size of the telescope beam, these regions are very much diluted, and with our data, we are mainly sensitive to the outer (and cold) protostellar envelope.

5 Radiative transfer modelling

In this section, we use non-LTE radiative transfer simulations to model the observed molecular lines. We applied the software LOC (Line Transfer with OpenCl, Juvela 2020), which is based on ray tracing. To reduce the complexity of the simulations, we made the simplifying assumption of a one dimensional, spherically symmetric grid for both of our cores.

For CCH, HCN, H¹³CN HCO⁺nd DCO⁺ we used the latest collision rate coefficients available on the Leiden Atomic and Molecular Database (LAMDA, Schöier et al. 2005; van der Tak et al. 2020; Dagdigian 2018; Hernández Vera et al. 2017; Denis-Alpizar et al. 2020; Pagani et al. 2012). In the case of HCN, H¹³CN and DCO⁺ we applied the hyperfine-resolved rate coefficients. For the collisional coefficients of H¹³CO⁺nd HC¹⁸O⁺ we used the rates of HCO⁺ from LAMDA and scaled them with the corresponding reduced mass of the isotopologues. For the ¹³C-and the deuterated species of CCH, a scaling of the CCH rates is not possible, as the hyperfine structures are different due to the different nuclear spins. Therefore, we excluded these molecules from the non-LTE modelling.

To the best of our knowledge, hyperfine resolved rate coefficients for DCN, HNC, HN¹³C and DNC do not exist¹. A simple approach to simulating lines of these species using radiative transfer would be to use the rate coefficients of the dominant isotopologues (HCN and HNC), neglecting the effect of isotopic substitution. However, it has been shown recently (Navarro-Almaida et al. 2023) that isotope effects can be important in the case of hydrogen cyanides and isocyanides. Indeed, the comparison of the different sets of collisional data revealed that isotopic substitution can lead to substantial changes (up to ~50%) at 10 K (the typical temperature of cold cores).

Accurate hyperfine-resolved rate coefficients for the excitation of DCN, HNC, HN¹³C and DNC induced by collisions with H2 are computed in this work. The scattering calculations are based on the HCN-H₂ (Denis-Alpizar et al. 2013) and HNC-H₂ (Dumouchel et al. 2011) interaction potentials corrected to consider the effect of isotopic substitution. The quantum mechanical close-coupling approach (Green 1975) and the almost exact recoupling method (Lanza & Lique 2014) was used to compute the collisional data for temperature ranging between 5 and 30 K. In the calculations, only the excitation by para-H2, the strongly dominant form of H₂ in cold ISM, has been considered. Details on the scattering calculations are provided in Appendix C.

The data underlying this article will be made available through the EMAA², LAMDA (Schöier et al. 2005; van der Tak et al. 2020), and BASECOL (Dubernet et al. 2013) data bases. They are also available on request.

Fig. 4 Profiles of the gas temperature (red), H2 number density (orange, in logarithmic scale) and infall velocity (blue, in units of 0.1 kms−1) for the Keto-Caselli model (dashed, KC Keto et al. 2015), the HDCRT model (dotdashed, HD Sipilä et al. 2022) of L1544, and the physical model of HH211 (dotted) derived in this work.

5.1 Physical structure of L1544

The physical structure of L1544 is well-studied (e.g. Keto et al. 2015, and references therein). To model the physical structure of the core, we applied the physical model presented by Keto et al. (2015, hereafter Keto-Caselli model). It describes an unstable quasi-equilibrium Bonnor-Ebert sphere with a peak central H₂ volume density of n₀ ≈ 10⁷ cm⁻³ and a central gas temperature of 6 K. The model provides the infall velocity, density and gas temperature structure as a function of distance from the centre of the core (see Fig. 4).

Recent studies prove the Keto-Caselli model to be unable to explain double-peaked profiles of optically thin lines (e.g. Redaelli et al. 2019; Ferrer Asensio et al. 2022). On top of that, Redaelli et al. (2022) show that to reproduce the blue asymmetry and high level of self-absorption of HCO⁺ (1–0), it is necessary to add a low-density and contracting envelope around the Keto-Caselli model. To test the effect of velocity variations in the outer core, we additionally used the HydroDynamics with Chemistry and Radiative Transfer model introduced by Sipilä et al. (2022, hereafter HDCRT model). The corresponding infall velocity, density and gas temperature structure are shown alongside the Keto-Caselli model in Fig. 4. The HDCRT model was designed to recreate the conditions of the Keto-Caselli model in the inner regions of the core where the gas-dust thermal coupling is strong. Unlike the Keto-Caselli model, the HDCRT model allows for expansion motions of the gas in the outer regions, caused by photoelectric heating.

5.2 Physical structure of HH211

There is no model of the physical structure of HH211 available in the literature. We computed it using data from the Herschel Science Archive (HSA). Following Harju et al. (2017), we used the archive Herschel/SPIRE dust continuum emission maps to derive the temperature and column density maps. From those, we extracted radial temperature and volume density profiles for the outer regions of the core. For the inner regions, we used power-laws following Crimier et al. (2010) and Motte & André (2001). For the velocity, we assumed free fall with the literature mass of the protostar (Lee et al. 2018). A detailed derivation of the physical structure of HH211 can be found in Appendix D.

5.3 Molecular abundances

We used the state-of-the-art gas-grain chemical model of Sipilä et al. (2019) to obtain an estimate of the fractional molecular abundances in the two cores. We applied the two-phase model of the chemical code, where the gas phase chemistry and the entire ice on the grains are active.

To account for the fact that L1544 is embedded in a molecular cloud, we assumed an external visual extinction of A_V = 2 mag. We employed the same initial abundances and chemical networks as in Sipilä et al. (2019) and the same standard values for various model parameters as in Giers et al. (2022), not recounted here for brevity. The abundance profiles are obtained by using the Keto-Caselli model and the HDCRT model as a static physical model of the core. We ran the chemical simulation, extracting the abundance profiles for various evolutionary times. In addition, we tested the abundances at the time of the best match between the Keto-Caselli model and the fiducial hydrodynamical simulation, as presented in Sipilä et al. (2022).

To derive the molecular abundances for HH211, we first adopted the physical conditions that roughly correspond to the edge of the HH211 model (n_H = 5 × 10³ cm⁻³, T = 15 K, A_v = 6 mag) and let the chemistry evolve over 10⁶ yr in those conditions, extracting the abundances at that time to use as the initial abundances for the protostellar core model.

The deuterium chemical model does not consider the isotope chemistry of carbon, and therefore we scaled the abundance profiles by the carbon isotopic ratio of the local ISM (68, Milam et al. 2005) to obtain estimates for the ¹³C-species. The same applies for the ¹⁸O isotopologue of HCO⁺ (557, Wilson 1999).

In addition to the molecular abundance profiles derived by chemical modelling, we applied a simple approach by assuming a constant abundance of the molecules throughout the core. This simplification allows to estimate the abundance levels of the molecules, especially in the case of the HDCRT model, and can be used for parameter-space explorations in future studies. However, in the case of L1544, Caselli et al. (2022) find evidence of almost complete freeze out within the central 1800 au. We took this into account by setting the abundance to zero for radii within 1500 au.

5.4 Results

For the comparison with the observed spectra, the synthetic ones were convolved with the respective corresponding telescope beam size (see Table A.1). To find the best fits to the observed spectra, we tested the different timesteps of the abundance profiles (10 logarithmically spaced steps between 10⁵ and 10⁶ yr of chemical evolution) and run through a range of constant abundances in log scale. The best fits were determined by optimising the line shape, line width and peak intensity of the synthetic spectra in accordance with the observations. A selection of the best-fit synthetic spectra derived with radiative transfer modelling using LOC are presented in Figs. 5–9. The respective labelling of the applied models is explained in Table 5.

For L1544, the results are split between the Keto-Caselli model and the HDCRT model. The labels KC (red solid line) and HD (blue solid line) denote the results of a constant abundance applied to the respective model, while STKC (cyan dashed line) and STHD (orange dashed line) refer to the abundance profiles derived with chemical modelling. The results produced by the abundance profile from the best match between Keto-Caselli and fiducial hydrodynamical model are described by the label HD0 (green dot-dashed line). For HH211, we plot the best-fit results of applying a constant abundance (blue solid line) and an abundance profile derived with the two-phase chemical model (red dotdashed line). The observed spectra are overplotted in black. The corresponding best-fitting abundances are plotted in Figs. E.1 and E.2 for L1544 and HH211, respectively.

To account for the level of turbulence in the cores, we used the average of the non-thermal components of the observed velocity dispersions. The values that best reproduce the line widths of the optically thin species are σ_turb = 0.1 km s⁻¹ (L1544, Keto-Caselli model), σ_turb = 0.17km s⁻¹ (L1544, HDCRT model) and σ_turb = 0.23km s⁻¹ (HH211). Separate values for the two different physical models of L1544 are necessary, because the outer radius of the HDCRT model is a factor of two larger compared to the Keto-Caselli model. Therefore, the model covers more turbulent gas.

To account for the additional velocity component at 8.22 km s⁻¹ observed towards HH211, we approximated it by a static Gaussian that is added to the modelling results afterwards. For this, we used the results of the two-component Gaussian fits listed in Table B.1. In the case of HNC, where a two-component Gaussian fit was not successful, we estimated the static Gaussian by eye to improve the readability of the modelled spectra.

5.4.1 L1544

The synthetic spectra of CCH are shown in Fig. 5, with each hyperfine component presented in a separate plot (numbered from low to high frequency). Towards L1544, the constant abundances applied to the two different physical models are able to reproduce the intensities of the more optically thin hyperfines. The line profiles of the more optically thick components are not reproduced by the Keto-Caselli model. The observed lines show rather symmetric double peaks, whereas the synthetic lines display a blue asymmetry. This behaviour was also observed for CS (2–1) by Tafalla et al. (1998) and might be an indication that these transitions trace a static envelope. Using the HDCRT model and a constant abundance, we are indeed able to reproduce the symmetric double-peak profiles, except for hyperfine component 5. The abundance profiles derived with chemical modelling largely underestimate the amount of CCH, and therefore fail to reproduce the lines in both the STKC and the STHD case. To reproduce the same intensities as in the KC case, the STKC abundance profile has to be multiplied by a factor of 20. The discrepancy between the simulated and observed CCH abundances is most likely related to the uncertainties in the chemistry of CCH, though the use of a simple 1D physical model may contribute as well. Future studies of the formation and destruction pathways of CCH are required to improve the chemical simulations, as well as the adoption of 3D models (Jensen et al. 2023).

The synthetic spectra of HCN, H¹³CN and DCN are shown in Fig. 6 and 7. The Keto-Caselli model fails at reproducing HCN with both a constant abundance and an abundance profile. The synthetic spectra show blue asymmetries in all hyperfines, indicating infall motion, and do not recreate the red asymmetry observed in the central hyperfine. This shows that the applied abundance profiles overestimate the amount of HCN in regions of infall and indicates that HCN might have a higher abundance in the low-density regions of the core. In addition, the modelled dip caused by the self-absorption does not reach down to zero level. This might be solved by adding an extended low-density envelope to the Keto-Caselli model, such as in the case of HCO⁺ (Redaelli et al. 2022). The HDCRT model, on the other hand, is able to produce the red asymmetry profiles for HCN by applying the STHD abundance profile. However, it results in red asymmetries for all hyperfines, reproducing the intensity of two hyperfines, but underestimating the weakest component by a factor of two. Applying a constant abundance results in a rather symmetric line profile for all components. The self-absorption dip of the synthetic spectrum goes down to almost zero level, as observed. A comparison of the STHD and HD0 abundance profiles (see Fig. E.1) shows that, in the case of HD0, HCN is more abundant at larger radii than in the case of STHD. This results in a deeper self-absorption dip. The distance between the two peaks is not reproduced (see Fig. 6), indicating that for HCN the turbulence level of 0.17 km s⁻¹ is overestimated in the hydrodynamical case.

In the case of H¹³CNsee Fig. 6), the Keto-Caselli model works very well and is able to reproduce the spectrum with both the KC and STKC abundance profile. The HDCRT model, however, struggles to reproduce the correct shapes of the self-absorbed hyperfines. Both the HD and STHD abundance profiles result in symmetric rather than blue asymmetric lines, which might hint at an overestimation of H¹³CN in the outer regions and/or a radially changing carbon isotopic ratio. In the case of HD0, the spectrum is severely underestimated. This might be a hint that at this specific timestep chosen by Sipilä et al. (2022), the carbon isotopic ratio is actually below 68. For DCN (see Fig. 7), the Keto-Caselli model can reproduce the intensities of the hyperfines within a factor of two. However, it fails to reproduce the shapes of the lines. Both the KC and the STKC abundance profiles predict a strong red asymmetry that is not observed. In addition, the synthetic spectra show higher self-absorption than the observed spectrum. This shows that the abundance profiles are actually overestimating the amount of DCN and that in reality there is less material available to trace the infall motions. The HDCRT model cannot reproduce the spectrum, neither the line shapes nor the intensities match the observations.

The synthetic spectra of HNC, HN¹³C and DNC are shown in Fig. 8. The collision rate coefficients for HN¹³C and DNC only consider the spin of the¹⁴N for the hyperfine splitting. Therefore, the modelling fails to reproduce the lowest frequency component of these two molecules. However, work is underway to address this issue, where the collisional rate coefficients are computed by also including the effect of ¹³C and D, respectively, which will allow to reproduce the full spectrum. Towards L1544, HNC can be reproduced quite well, though the Keto-Caselli model seems to perform better than the HDCRT model at reproducing the observed intensity peaks. However, due to the greater physical size of the structure in the HDCRT model, the dip at 7.2 km s⁻¹ is better reproduced than in the case of the Keto-Caselli model. At high velocity there seems to be an issue with the hyperfine component that has the lowest statistical weight. This might be caused by an underestimation of the general abundance of the molecule and subsequently of the depth of self-absorption in the lower velocity hyperfine components.

In the case of HN¹³C the Keto-Caselli model is able to reproduce the spectrum with both constant abundance and abundance profile, except for the high velocity component. The HDCRT model, however, results in a stronger self-absorption and therefore underestimates the observed intensity. In the case of DNC, the STKC abundance profile, applied to the Keto-Caselli model, results in a good match with the observed line shape and intensity (excluding the high velocity component). The synthetic spectrum produced by the constant abundance largely overestimates the amount of DNC in the outer regions of the model, resulting in a deep dip caused by self-absorption of the strongest hyperfines that is not observed. The HDCRT model is not able to reproduce DNC, it underestimates both intensity and line width. A comparison of the abundance profiles (see Fig. E.1) shows that the peaks of the STHD and the HD0 profiles are at least one magnitude lower than for the STKC profile, which leads to an underestimation of the molecular abundance.

The synthetic spectra of HCO⁺ are shown in Fig. 9. Towards L1544, our models fail to reproduce the strong self-absorption and asymmetry of HCO⁺ However, the case of this molecule is extensively studied in Redaelli et al. (2022). By adding a low-density (27 cm⁻³), contracting envelope extending out to 1 pc, they manage to reproduce the line shape of the transition. In the case of H¹³CO⁺ a constant abundance profile applied to the Keto-Caselli model can reproduce the observed spectrum. The modelled line profile shows a stronger blueshift and a higher intensity than observed. However, a constant abundance of H¹³CO⁺ with complete freeze out in the central 1500 au is clearly too simplistic for this molecule, present not only in the core, but also in the surrounding cloud. In the case of the HDCRT model, a constant abundance of H¹³CO⁺ results in an almost symmetric double-peak profile. The abundance profiles derived by chemical modelling severely underestimate the observed intensity and cannot reproduce the spectrum. This might indicate that in this case the carbon isotopic ratio is less than 68.

For HC¹⁸O⁺ the line cannot be reproduced by either of the models. The HD abundance profile roughly reproduces the line width, but fails at the flattened top caused by CO depletion in the central regions of the core. Like in the case of H¹³CO⁺ the synthetic spectra derived with the abundance profiles underestimate the intensity of the line and indicate that the assumption of the isotopic ratio to be 557 might be not correct in this case. For DCO⁺ the observed line shape cannot be reproduced. The simulation is considering the overlapping hyperfine structure of the molecule and the modelled spectra show a rather flat peak. This is likely caused by individual, self-absorbed hyperfine components.

Fig. 5 Synthetic spectra of CCH obtained with LOC (top rows: L1544, bottom row: HH211). Each hyperfine line is plotted separately, numbered from low to high frequency. For L1544, the results are splitted between using the Keto-Caselli model (top) and the HDCRT model (centre) as approximation for the physical structure of the core. Towards L1544, the best-fitting abundances and timesteps are 2 × 10−8 (KC), 1e5 yr (STKC), 3 × 10−8 (HD), 1e5 yr (STHD). Towards HH211, the best-fitting abundances and timesteps are 2 × 10−9 (const), 1e6 yr (2-phase).

Fig. 6 Synthetic spectra of HCN (left) and H13CN (right) obtained with LOC (top rows: L1544, bottom row: HH211). For L1544, the results are splitted between using the Keto-Caselli model (top) and the HDCRT model (centre) as approximation for the physical structure of the core. Towards L1544, the best-fitting abundances and timesteps are 2 × 10−8 (KC), 7.9e5 yr (STKC), 1 × 10−8 (HD), 1.3e5 yr (STHD) for HCN; 2 × 10−10 (KC), 6e5 yr (STKC), 3 × 10−10 (HD), 1e5 yr (STHD) for H13CN. Towards HH211, the best-fitting abundances and timesteps are 2 × 10−9 (const), 1e5 yr (2-phase) for HCN; 3 × 10−11 (const), 1.3e5 yr (2-phase) for H13CN.

Fig. 7 Synthetic spectra of DCN obtained with LOC (top rows: L1544, bottom row: HH211). Towards L1544, the best-fitting abundances and timesteps are 1 × 10−9 (KC), 1.3e5 yr (STKC), 1 × 10−9 (HD), 9.8e5 yr (STHD) for DCN. Towards HH211, the best-fitting abundances and timesteps are 4 × 10−10 (const), 1e5 yr (2-phase) for DCN.

5.4.2 HH211

The synthetic spectra of CCH are shown in Fig. 5, with each hyperfine component presented in a separate plot (numbered from low to high frequency). The observed spectrum can be reproduced assuming a constant molecular abundance of 2 × 10⁻⁹ with respect to H₂ throughout the core. The intensity ratios between the hyperfine components are reproduced within a factor of two in the abundance. The synthetic spectra resulting from abundance profiles derived with the 2-phase chemical model, where gas phase and grain chemistry are considered, underestimate the observed lines by a factor of two to five. This supports the need to constrain the reactions involved in the chemistry of CCH to improve the chemical simulations.

The synthetic spectra of HCN, H¹³CN and DCN are shown in Figs. 6 and 7. Towards HH211, we are not able to reproduce the observed spectrum of HCN. Applying the derived physical model of the protostellar core produces a self-absorbed line profile with blue asymmetry for HCN that is not observed. The¹³ C isotopologue, however, can be reproduced with the correct intensities, line shapes and line widths with both the constant abundance and the abundance profile. In the case of DCN, only the weakest hyperfine component is reproduced. The more optically thick hyperfines, however, are overestimated by the model by a factor of two to three. This might indicate that DCN is more extended than what the model covers currently, resulting in more self-absorption than produced by the model. This should be addressed in future studies.

The synthetic spectra of HNC, HN¹³C and DNC are shown in Fig. 8. The modelled HNC shows an asymmetric self-absorption that is not observed. Also here, the linewidth towards high velocity is underestimated, similar to HNC in L1544. Taking into account the missing hyperfine component at high velocity (~ 10 km s⁻¹), HN¹³Cnd DNC are well reproduced by a constant abundance of 1 × 10⁻¹¹ and 7 × 10⁻¹¹, respectively. For HN¹³C the synthetic spectrum derived with an abundance profile shows that also in this case, the assumption of ¹²C/¹³C = 68 is underestimating the amount of HN¹³C present. For DNC, the use of an abundance profile results in a self-absorbed line profile, most likely due to an overestimation of the molecular abundance in the outer regions.

The synthetic spectra of HCO⁺ and isotopologues are shown in Fig. 9. Towards HH211, all isotopologues of HCO⁺ can be reproduced in shape and intensity, using a constant abundance profile. The line shape of the main species can be modelled, however, the constant abundance used to achieve this is of the same magnitude as for DCO⁺ This indicates that the spectrum of HCO⁺ is heavily self-absorbed. When applying an abundance profile, this self-absorption is shown in the synthetic spectrum. The same happens in the case of DCO⁺ indicating that this molecule might be slightly optically thick as well. The ¹³C and ¹⁸O isotopologues cannot be reproduced with the HCO⁺ abundance profile scaled down by the respective isotopic ratio. This highlights the necessity of new constraints on the isotopic abundance ratios in varying physical conditions.

6 Discussion

6.1 Dilution of¹³ in CCH

We find that C¹³CH is more abundant than ¹³CCH in L1544 and in HH211 (see Table 4). This characteristic appears to be common in cold dark clouds, and was observed already in the starless cores TMC-1, L1521B and L134N (Turner 2001; Sakai et al. 2010; Taniguchi et al. 2019), and in the protostellar cores L1527 and L483 (Sakai et al. 2010; Agúndez et al. 2019). Taniguchi et al. (2019) conclude that this difference in abundances is caused during the formation pathway of the molecule and the isotopomer-exchange reaction after it is formed. This reaction is exothermic:

$${}^{13}\text{CCH}+\text{H}\rightleftharpoons {\text{C}}^{13}\text{CH}+\text{H}+8\text{\hspace{0.17em}K},$$(5)

thus, at low temperatures the forward reaction is more efficient, increasing the abundance of the C¹³CH isotopomer. Chemical models of Furuya et al. (2011) predicted an increase of the C¹³CH/¹³CCH ratio with pre-stellar evolution. Literature values of this ratio are collected in Fig. 10, which shows the ratios observed towards the starless cores L1521B, TMC-1, L134N, the very evolved pre-stellar core L1544, and the protostellar cores HH211 and L1527. However, due to the large uncertainties, no conclusive statement can be made. More observations are needed to confirm any trends.

The carbon isotopic ratios observed for CCH are 200 ± 50 and 90 ± 20 (L1544), and 160 ± 60 and 100 ± 30 (HH211) for ¹³CCH and C¹³CH respectively. This is significantly higher than the canonical value for the ISM, 68, and is caused by the dilution of ¹³C in carbon-chain molecules in the local ISM. In dark clouds, ¹³C is mainly locked in ¹³CO due to the reaction

$$\text{CO}{+}^{13}{\text{C}}^{+}{\to }^{13}\text{CO}+{\text{C}}^{+}+35\text{\hspace{0.17em}K}.$$(6)

As this reaction is exothermic, the forward reaction is more efficient at low temperatures. Therefore, molecules such as CCH or c-C₃H₂ that are produced from C⁺ show a significantly higher ¹²C/¹³C ratio in molecular clouds. This dilution of ¹³C species was predicted by Langer et al. (1984), and confirmed by previous studies (Sakai et al. 2010; Agúndez et al. 2019; Taniguchi et al. 2019; Yoshida et al. 2019). In Colzi et al. (2020), this behaviour is also predicted for CH.

Fig. 8 Synthetic spectra of HNC and isotopologues obtained with LOC (top rows: L1544, bottom row: HH211). For L1544, the results are splitted between using the Keto-Caselli model (top) and the HDCRT model (centre) as approximation for the physical structure of the core. Towards L1544, the best-fitting abundances and timesteps are 2 × 10−8 (KC), 1e5 yr (STKC), 5 × 10−8 (HD), 1e5 yr (STHD) for HNC; 2 × 10−10 (KC), 1.3e5 yr (STKC), 2 × 10−10 (HD), 9.8e5 yr (STHD) for HN13C; 3 × 10−9 (KC), 1e5 yr (STKC), 8 × 10−10 (HD), 9.8e5 yr (STHD) for DNC. Towards HH211, the best-fitting abundances and timesteps are 5 × 10−10 (const), 1e5 yr (2-phase) for HNC 1 × 10−11 (const), 1e5 yr (2-phase) for HN13C; 7 × 10−11 (const), 1e6 yr (2-phase) for DNC.

Fig. 9 Synthetic spectra of HCO+ and isotopologues obtained with LOC (top rows: L1544, bottom row: HH211). For L1544, the results are splitted between using the Keto-Caselli model (top) and the HDCRT model (centre) as approximation for the physical structure of the core. Towards L1544, the best-fitting abundances and timesteps are 7 × 10−10 (KC), 9.8e5 yr (STKC), 3 × 10−9 (HD), 9.8e5 yr (STHD) for HCO+; 2 × 10−10 (KC), 1e5 yr (STKC), 4 × 10−10 (HD), 1e5 yr (STHD) for H13CO+; 6 × 10−12 (KC), 1e5 yr (STKC), 1 × 10−11 (HD), 1e5 yr (STHD) for HC18O+ ; 2 × 10−10 (KC), 1.7e5 yr (STKC), 4 × 10−10 (HD), 1.7e5 yr (STHD) for DCO+. Towards HH211, the best-fitting abundances and timesteps are 1 × 10−11 (const), 1e6 yr (2-phase) for HCO+; 7× 10−12 (const), 1e5 yr (2-phase) for H13CO+; 4× 10−13 (const), 1e5 yr (2-phase) for HC18O+; 2× 10−11 (const), 1e6 yr (2-phase) for DCO+.

6.2 Deuterium fractionation

Due to the dilution of the ¹³C isotopologues of CCH, the derived deuterium fractions have to be interpreted as upper limits, as the applied ¹²C/¹³C ratio of 68 is underestimating the real value. However, as it was possible to derive the column density of CCH using emission of the main species directly, we use these results for a comparison of the deuteration in the two cores. CCH shows a high and similar level of deuteration in both cores (≈ 10%) and is therefore consistent with other carbon chains such as c-C₃H₂ (e.g. Chantzos et al. 2018). Measurements towards other cores report lower deuteration levels (L183: 0.06(4), Turner 2001; TMC-1: 0.05(2), Turner 2001; L1527: 0.04(1), Yoshida et al. 2019). Evolutionary and environmental effects might play a role in this difference, as well as differences in the excitation temperatures used to derive the column densities.

The deuterium fraction of HCN appears to be slightly larger in HH211 with respect to L1544. However, due to the uncertainties in the carbon isotopic ratio of HCN, these numbers have to be interpreted with caution. Therefore, no conclusive statement can be made based on our errorbars, and further studies are necessary to eventually confirm this trend. The D/H ratio for HCN of the cold dark cloud L183 (0.05(2), Turner 2001) is consistent with our measurements towards L1544. Roberts et al. (2002) derived slightly lower values in HH211, N(DCN)/N(HCN) = 0.038(8), which can be explained by different approximations used in the calculations. The measurements towards the protostellar core L1527 are consistent with our observations (0.05(1), Yoshida et al. 2019).

For HNC, we derive a deuterium fraction of 0.08(2) for both L1544 and HH211. Previous measurements of the deuteration in the two cores show moderate deuterium fractions (L1544: 0.03, Hirota et al. 2003; HH211: 0.07, Imai et al 2018). However, the authors considered the emission of HN¹³C and DNC to be optically thin and apply simple one-component Gaussian fits to obtain the line parameters. Especially in the case of DNC this can lead to an underestimation of the column density, and subsequently of the deuterium fraction. Observations towards other cores show moderate deuterium fractions (L183: 0.05(2), Turner 2001; L1527: 0.045, Yoshida et al. 2019).

The deuterium fractions of HCO⁺ derived towards L1544 and HH211 are consistent with previous measurements in similar objects (see e.g. Turner 2001; Jørgensen et al. 2004; Koumpia et al. 2017; Yoshida et al. 2019). For L1544, the D/H ratios we derived using H¹³CO⁺ and HC¹⁸O⁺ are in agreement with previous observations by Caselli et al. (2002a) and Redaelli et al. (2019), who report a ratio of around 4%.

Towards HH211, the result derived using the ¹⁸O isotopo-logue is increased by a factor of two compared to the D/H ratio derived from H¹³CO⁺ This is most likely caused by an overesti-mation of the elemental carbon isotopic ratio. Chemical models predict that the HCO⁺/H¹³CO⁺ ratio is actually lower than the canonical ¹²C/¹³C ratio (Colzi et al. 2020), because of isotopic exchange reactions that are important at low temperatures. The canonical ratio of ¹³C/¹⁸O is 8.2, whereas our observations give 10(1) and 20(4) for L1544 and HH211, respectively. This shows that especially towards HH211, the carbon fractionation of HCO⁺ is higher than in the local ISM, resulting in HCO⁺/H¹³CO⁺ < 68. Hence, by using 68, the deuterium fraction derived from H¹³CO⁺ is underestimated.

The comparison of the D/H ratio derived from HC¹⁸O⁺ in the two cores shows that the ratio towards the protostellar core is higher by a factor of almost two. This difference is caused by some effect that is either increasing or decreasing the observed column densities of the isotopologues. The latter one could be an effect of isotope-selective photodissociation which is impacting the ¹⁸O more than the D isotopologue. The higher abundance of CO leads to a stronger self-shielding of the molecule compared to C¹⁸O. This results in a local enhancement of CO/C¹⁸O, that is likely larger in HH211 than in L1544. This is reflected in an increased D/H ratio towards HH211, as DCO⁺ is created directly from CO, and therefore it is more efficiently formed than HC¹⁸O⁺

In conclusion, we do not see a general trend in the deuteration efficiency in simple molecules going from a pre-stellar to a protostellar core, as it was indicated previously by observations of c-C₃H₂ (Chantzos et al. 2018). This might be due to a more efficient deuteration of c-C₃H₂ happening on the surface of the dust grains in the pre-stellar phase that is not applicable for the simple molecules studied in this work. Additionally, the fact that the deuterium fraction is similar in two cores in different environments seems to suggest that the deuterium fractiona-tion is not sensitive to the initial conditions of its surroundings. Similar results have been found for N₂H⁺ (Crapsi et al. 2005; Emprechtinger et al. 2009) and c-C3H2 (Chantzos et al. 2018). Nevertheless, more comprehensive surveys including multiple pre-stellar and protostellar cores and covering more molecules are necessary to provide larger statistics and confirm any trends.

Fig. 10 C13CH/13CCH ratio of starless/pre-stellar (black) and protostellar (red) cores. L1521B: Taniguchi et al. (2019); L134N: Taniguchi et al. (2019); TMC-1: Sakai et al. (2010); L1527: Yoshida et al. (2019).

6.3 Radiative transfer modelling

6.3.1 L1544

Our radiative transfer simulations show that the Keto-Caselli model works fine to explain moderately optically thin lines such as H¹³CN HN¹³C DNC, and H¹³CO⁺ hat trace regions with moderate to high density (few 10⁵ cm⁻³). The model shows problems when it comes to optically thin lines like HC¹⁸O⁺ where the line shape is affected by the freeze-out of CO in addition to contraction motions. In the case of CCH, the model fails to reproduce the line shapes of the optically thin hf components, which show a red asymmetry. As the lines are optically thin, they do not trace any expansion motion like it was seen for HCN. Instead, this is likely an effect of the asymmetric distribution of the molecule across the core. In L1544, carbon chains are known to peak towards the south-east of the core (Spezzano et al. 2017; Giers et al. 2022). Due to the velocity gradient across the core (Spezzano et al. 2016), the spectra observed towards the south are redshifted with respect to the rest velocity of the system. As CCH is most probably peaking in the south, where the gas motion is redshifted, the contribution of the redshifted lines is more prominent, resulting in an asymmetric line shape. This behaviour is also observed for cyanopolyynes by Bianchi et al. (2023). As we assume a spherically symmetric core in the simulations, and therefore a symmetric distribution of the molecule, this line shape is not reproducible.

The HDCRT model, on the other hand, struggles to explain the blue asymmetry observed in the moderately optically thin lines. This is most likely caused by an overestimation of the molecular abundances in the outer regions, resulting in symmetric double-peak profiles that trace more static rather than infalling layers.

In most cases, the simulations fail to reproduce the observed intensities of the rarer isotopologues when applying the abundance profiles of the main isotopologues scaled by the respective isotope ratio. This highlights the importance of considering the carbon and oxygen isotope chemistry in the chemical modelling, which will be addressed in future work.

In the case of the optically thick and more widespread main isotopologues, the results are better reproduced when applying the HDCRT model, while the Keto-Caselli model struggles to explain the lines. The only exception here is HNC, where the Keto-Caselli model is sufficient to reproduce the observed intensities. This might be an indication that HNC is actually tracing slightly different layers compared to HCN. In the case of CCH, the HDCRT model works fine for the optically thick components that trace a static inbetween layer.

In the case of HCN, the radiative transfer modelling in this work showed that the Keto-Caselli model is not sufficient to explain the line shapes and intensities of the molecule’s hyperfine components observed towards L1544. Instead, it is necessary to use a physical model that shows expansion motions in the outer layers to reproduce the red asymmetry in the central hyperfine component of HCN. Contrary to this, Redaelli et al. (2022) show that for the blue asymmetry and strong self-absorption observed in HCO⁺ it is necessary to add an extended low-density and contracting envelope to the Keto-Caselli model. One possibility to solve this discrepancy could be that HCN is tracing an even outer layer than HCO⁺ As HCN can be efficiently formed also without CO, the layers where this process occurs could be warm enough to be associated with expansion motions.

6.3.2 HH211

The physical model derived in this work struggles to explain the very optically thick main isotopologues (HCN, HNC, HCO⁺). Both a constant abundance across the core and the abundance profiles predict high molecular abundances in the outer regions of the core. This leads to a strong self-absorption in the synthetic spectra that is not observed. To get a more realistic estimate of the molecular abundances across the core, it would be necessary to run more tests of the physical conditions and the parameter space used in the chemical modelling. However, this is beyond the scope of this paper.

Furthermore, very abundant and widespread molecules like HCO⁺ HCN, and possibly also HNC also trace the outflow of the protostar and the extended structures of the core’s environment. These kind of large-scale structures are not considered in the simplistic and spherically symmetric physical model derived in this work and will likely be different for the sources studied here. L1544 is located in an isolated region in Taurus. In contrast, HH211 is forming in an active region within Perseus, with low density material traced by ¹³CO spreading over several km s⁻¹ around HH211 (Sun et al. 2006). The complex and turbulent environment towards HH211 is reflected in the line profiles for the more abundant species. For instance, our observations detect a second velocity component blueshifted by 1 km s⁻¹, which is associated with large-scale emission extending to the north-east of HH211 (Sun et al. 2006). To correctly model the observed line shapes of molecular tracers like HCO⁺ HCN, and HNC, these structures have to be taken into account in the density, velocity, and temperature profiles of the core.

On the other hand, we are able to reproduce all optically thin and moderately optically thick molecular lines observed towards the protostellar core. For these spectra, the simplifying assumptions of spherical symmetry and constant abundance throughout the core seem to be sufficient to reproduce the molecular emission. This is most likely due to the fact that these molecules are less abundant and trace higher-density regions, and are therefore less affected by the kinematics of the large-scale structures surrounding the core as well as the extended outflow emission.

7 Conclusions

We presented a survey of ground-state rotational lines of simple molecules, which allow us to compare the levels of deuteration in the very evolved pre-stellar core L1544 and the very young protostellar core HH211. In a non-LTE approach, we used radiative transfer simulations and molecular abundance profiles derived from chemical modelling to reproduce the observed molecular spectra. By applying new hyperfine-resolved collisional rate coefficients, we take into account the specific spectroscopy of the respective isotopologues.

Our main results can be summarised as follows:

• The similar levels of deuteration show that the deuterium fractionation seems to be equally efficient towards both cores. We do not see a general trend in the level of deuteration when going from a pre-stellar core to a protostellar core, suggesting that the protostellar envelope still retains the chemical composition of the original pre-stellar core. In addition, the deuterium fraction seems to be independent of the initial conditions present in the molecular clouds where the two cores are embedded;

• Towards HH211, the D/H ratio of HCO⁺ derived using the¹⁸ O isotopologue is higher by a factor of two compared to L1544. This is likely an effect of isotope-selective photodissociation, creating a local enhancement of CO/C¹⁸O in the protostellar core;

• The ¹³C dilution of carbon-chain species in dark clouds leads to an increased ¹²C/¹³C ratio for CCH. On the other hand, the increased ¹³C/¹⁸O ratio of HCO⁺ indicates a decreased ¹²C/¹³C ratio for the molecule, caused by the isotopic exchange reaction being more effective at low temperatures. This highlights the uncertainties when dealing with ¹³C isotopologues and the influence of the applied carbon isotopic ratio;

• The central hyperfine component of HCN (1–0) observed towards L1544 shows a red asymmetry, tracing expansion motions caused by external heating in the outer layers of the core. This stands in contrast to HCO⁺ that has been shown to trace infall motion in the outer layers and low-density envelope of L1544. This discrepancy could be explained by HCN tracing an even outer layer than HCO⁺

• The radiative transfer modelling showed that the new colli-sional rate coefficients for HN¹³C and DNC are incomplete because they only consider the splitting by ¹⁴N and therefore miss a part of the hyperfine structure. This issue is addressed by work underway that includes the effect of ¹³C and D;

• The radiative transfer modelling of HH211 using the physical structure derived in this work is successful for optically thin emission lines, but shows problems with optically thick lines that also trace more complex structures and are influenced by them;

• The radiative transfer modelling of L1544 applying the Keto-Caselli model works well for moderately optically thin lines that trace inner layers of the core, but shows problems with optically thin lines influenced by CO depletion (HC¹⁸O⁺ ) and asymmetric distribution across the core (CCH);

• The radiative transfer modelling of L1544 applying the HDCRT model is successful at reproducing the optically thick hyperfine components of CCH, which show symmetric double-peak profiles, and the red asymmetry observed in the central hyperfine component of HCN (1–0);

• The modelling results of both cores show that to correctly model emission lines, it is crucial to include the outer layers of the cores to consider the effects of extended structures.

Future projects with more comprehensive surveys, including multiple pre-stellar and protostellar cores and covering more molecules, will be able to provide larger statistics and confirm the trends reported in this work. New detailed chemical models including the time-dependent variations in ¹²C/¹³C ratios (Colzi et al. 2020; Sipilä et al. 2023) will aid in further constraining the lines of C-containing molecules, including isotopologues. Additionally, modelling of molecules such as CS or H₂CO will help to further constrain the diffuse envelope surrounding the pre-stellar core L1544.

Acknowledgements

We wish to thank the anonymous referee for their constructive comments. K.G. thanks Maria José Maureira for useful discussions. S.S. and K.G. wish to thank the Max Planck Society for the Max Planck Research Group funding. All others authors affiliated to the MPE wish to thank the Max Planck Society for financial support. The Onsala Space Observatory national research infrastructure is funded through Swedish Research Council grant no. 2017-00648. C.T.B. and F.L. acknowledge financial support from the European Research Council (Consolidator Grant COLLEXISM, Grant No. 811363).

Appendix A Observed lines

Fig. A.1 Ground-state rotational lines observed towards the pre-stellar core L1544. The respective hyperfine transition is given in the upper left corner of each plot. The dotted line indicates the rest velocity of the system (7.2 km s−1 ). The red curves show the fits used in the LTE analysis.

Fig. A.2 Ground-state rotational lines observed towards the protostellar core HH211. The respective hyperfine transition is given in the upper left corner of each plot. The dotted line indicates the rest velocity of the system (9.2 km s−1). The red curves show the fits used in the LTE analysis.

Table A.1 lists all emission lines observed towards L1544 and HH211, along with the respective transition, frequency and telescope beam size. Fig. A.1 and Fig. A.2 show the corresponding spectra for L1544 and HH211, respectively, with the hyperfine transition given in the top left corner.

Appendix B Parameters of the additional velocity component towards HH211

Table B.1 presents the best–fit parameters of the additional velocity component observed towards HH211, along with the corresponding column density. The parameters were obtained by applying a two-component Gaussian fit to the spectral lines. In the case of HN¹³C and DNC, we applied a three–component Gaussian to account for their more separate hyperfine component at high velocity. The additional velocity component is located at approximately 8.22km s⁻¹.

Table B.2 compares the level of deuteration between the additional component (add.) and the main component. Within errorbars, the deuterium fractionation is similar for CCH, HCN, and HNC, indicating that main and additional component trace the same gas along the line–of–sight. For HCO⁺ the D/H ratios differ by a factor of two between the two components. However, with this little amount of data on the additional velocity component, we cannot make any conclusive statement about this. Most likely, the additional velocity components of HCO⁺ and isotopo– logues are tracing the same gas as the other molecules in this study, plus some additional gas or structure.

Appendix C Collisional rate coefficients

Within the Born-Oppenheimer approximation, the interaction-potential is the same for HCN-H₂ and DCN-H₂ (HNC-H₂, DNC-H₂, HN¹³C-H₂) and only depends on the mutual distances of the atoms. Hence, the only difference between the interaction potentials of two different isotopologues with H₂ is the position of the centre of mass taken for the origin of the Jacobi coordinates used to describe the geometries of the system. Then, we employed the HCN-H₂ (Denis-Alpizar et al. 2013) and HNC-H₂ (Dumouchel et al. 2011) potential energy surfaces (PESs) corrected to consider the effect of isotopic substitution.

In the following, the rotational level of the targets (DCN, HNC, HN¹³C, and DNC) will be denoted j. Only the coupling of the molecular rotation with the nuclear spin ((I = 1) of the nitrogen atom will be considered. The coupling results in a weak splitting of each rotational level j, into three hyperfine levels (except for the j = 0 level which is split into a single level). Each hyperfine level is designated by a quantum number F (F = I + j) varying between |I - j| and I + j. Inly collisional excitation by para-H₂ in its ground rotational state has been considered.

Rotational (i.e. nuclear-spin free) scattering matrix, cross sections and rate coefficients for the excitation of DCN, HNC, HN¹³C and DNC (the isotopologues of HCN and HNC) induced by collisions with H₂ were computed by Navarro-Almaida et al. (2023) for temperatures up to 30 K. These data were calculated using the quantum mechanical close-coupling approach (Green 1975). For more computational details, we refer the readers to the work of Navarro-Almaida et al. (2023).

To account for the hyperfine structure due to the nitrogen nuclear spin of the isotopologues mentioned above, we used the scattering matrix computed by Navarro-Almaida et al. (2023) and applied a recoupling method (Alexander & Dagdigian 1985; Lanza & Lique 2014). Therefore, we derived hyperfine resolved rate coefficients for the 25 low-lying energy levels, that is (j, F) ≤ (8, 9), for temperatures up to 30 K.

Fig. C.1 displays the temperature dependence of the hyper-fine resolved excitation rate coefficients of HCN, DCN, HNC, DNC and HN¹³C induced by collision with H₂ for the (j, F) = (2, 1) → (1, F′) transitions, where F′ = 0, 1, 2. The upper panels show that the rate coefficients for DCN-H₂ collisions exhibit systematically a higher magnitude than those for HCN-H₂ collisions for all transitions, the differences being of the order of roughly 25% – 35%. At the opposite, the collisional data for DNC-H₂ collisions can either overestimate or underestimate the ones for HNC-H₂ collisions depending on the transitions (see lower panels). It is also interesting to note that the HN¹³C-H₂ collisional data increase with increasing temperature while those for HNC-H₂ (DNC-H₂) collisions decrease. All this comparison clearly demonstrates the need of using isotopologues specific collisional data for modelling HCN/HNC isotopologues’ observations.

Fig. C.1 Temperature dependence of hyperfine resolved rate coefficients of HCN and DCN (upper panels) and HNC, DNC and HN13C (lower panels) for the ( j, F) = (2, 1) → (1, F′) transitions. The data of HCN are from Magalhães et al. (2018).

Appendix D Derivation of the physical structure of HH211

To derive the physical structure of HH211, we used the Herschel SPIRE dust continuum emission maps of the core from the Her-schelScience Archive (HSA). Using the three SPIRE bands at 250, 350 and 500 µm, we derived the H₂ column density and the dust temperature map, following Harju et al. (2017). Therefore, we smoothed the 250 µm and 350 µm images to the resolution of the 500 µm image (35″), and then fitted a modified black-body to each pixel. We applied a dust emissivity spectral index of β = 2, and a dust emissivity coefficient per unit mass of gas of ĸ_250µm= 0.1 cm² g⁻¹ (Hildebrand 1983). The two resulting maps are shown in Fig. D.1, with the direction and size of the jet indicated by red and blue arrows.

From the H2 column density map we derived a radial column density profile. For this, we masked out the structures north and south-east of the core, and the filamentary structure trailing the core in the south-west (see Fig. D.1). To characterise the radial profile, we fitted a Plummer-like profile (Plummer 1911; Whit-worth & Ward-Thompson 2001), that is modified by a constant term to account for the background column density (Launhardt et al. 2013):

$$N_{H_2}\left(r\right)=\frac{N_{H_2,0}}{{\left[1+{\left(r/R_{\text{flat}}\right)}^2\right]}^{\frac{p-1}{2}}}+N_{\text{out}},$$(D.1)

where $N_{H_2}$,0 is the central column density, R_flat the characteristic radius of the flat inner portion of the density profile, and power-law index p. The three fit parameters were then used to derive the radial volume density profile (André et al. 2016):

$$n_{\text{outer}}\left(r\right)=\frac{n_{H_2,0}}{{\left[1+{\left(r/R_{\text{flat}}\right)}^2\right]}^{p/2}}.$$(D.2)

The central volume density $n_{H_2}$,0 is related to the central column density by $n_{H_2}$,0 = $N_{H_2}$,0/ (A_PR_flat), where A_P is a constant factor related to the inclination angle of the core to the plane of sky, $A_P=\frac{1}{\cos i}\times B\left(\frac{1}{2},\frac{p-1}{2}\right)$, and B is the Euler Beta function. Assuming i = 0°, we derived N₀ = 9.7(6) × 10²² cm⁻², flat = 8(1) × 10⁻² pc, p = 3.0(4), N_out = 6(1) × 10²¹ cm⁻², and n₀ = 1.9(3) × 10⁵ cm⁻³.

For the inner regions of the core, we adapted the frequently used physical model for IRAS 16293-2422 A/B (Crimier et al. 2010), where the radial density distribution is described by a power law, n(H₂) = r^−1.8. After normalising the profile with the central pixel of the Herschel column density map, we arrived at a volume density profile of

$$n_{\text{inner}}\left(r\right)=n_1\cdot {\left(\frac{r}{r_1}\right)}^{-1.5},$$(D.3)

where n₁ = 1.89 × 10⁵ cm⁻³ and r₁ = 0.0128 pc are the starting values of n_outer. This Shu-like density distribution is for example also used in Quénard et al. (2018).

Fig. D.1 Analysis of Herschel SPIRE data towards HH211. Top: H2 column density map of HH211, derived using the Herschel SPIRE images at 250, 350 and 500 µm. The red and blue arrows indicate the direction and size of the molecular outflow/jet from the protostellar core (marked by a black triangle). The Herschel beam size is shown as a white circle, the average beam size of our observations is shown as red circle. Bottom: Dust temperature map of HH211, derived using the Herschel SPIRE images at 250, 350 and 500 µm.

To parameterise the radial temperature profile, we applied an empirical fit following Arzoumanian et al. (2011):

$$T\left(r\right)=T_{\text{out}}-\frac{\text{Δ}T}{{\left(1+{\left(\frac{r}{r_{\text{flat}}}\right)}^2\right)}^{q/2}}.$$(D.4)

We fixed the values for T_out and ΔT to 21 K and 8 K, respectively, resulting in best-fit parameters of r_flat = 0.17(1) pc and q = 0.38(3). The profile is decreasing towards the centre, with a central temperature of around 13 K, and does not show an internal heating source. This is caused by the low spectral resolution of the Herschel temperature map, where one pixel covers around 6000 au. For example, Launhardt et al. (2013) show that embedded protostars only influence the temperature in the inner 5000 au of a globule.

To account for the central heating, we adopted the radial dust temperature profile for protostellar envelopes introduced by Motte & André (2001):

$$T_{\text{dust}}\left(r,L_{*}\right)\approx 38\text{\hspace{0.17em}K}\times {\left(\frac{r}{100\text{\hspace{0.17em}AU}}\right)}^{-0.4}{\left(\frac{L_{*}}{1L_{\odot }}\right)}^{0.2}.$$(D.5)

Fig. D.2 Analysis of the H2 column density map and the temperature map derived from Herschel data. Top: Radial H2 column density profile (red), overlaid with the fitted Plummer-like profile (blue). Given in green are the data points averaged over shells of 11 arcsec, the size corresponding to 1/3 of the Herschel beam size. Bottom: Radial temperature profile (red), overlaid with the best fit (blue) and the averaged shells (green).

Assuming that internal heating by the accreting protostar dominates the thermal balance of the envelope, and the envelope is optically thin to the bulk of the infrared radiation. To derive the internal luminosity of the protostar, we followed Dunham et al. (2008), using the 70 µm flux as an estimate:

$$L_{\mathrm{int}}=3.3\times {10}^8{F}_{70}^{0.94}{L}_{\odot }.$$(D.6)

Using the 70 µm flux of HH211 published by Dunham et al. (2015), we derived a protostellar luminosity of L_* = 1.2(2) L_⊙. For the purpose of our radiative transfer simulations, we assumed the gas temperature in the core to be equal to the dust temperature. This is reasonable, as gas and dust are closely coupled at densities above 10⁴ cm⁻³ (Goldsmith 2001).

For the infall velocity profile, we followed Crimier et al. (2010) and assumed a free fall velocity gradient,

$$\upsilon \left(r\right)=\sqrt{\frac{2gM}{r}},$$(D.7)

where the mass M is given by the protostellar mass, 0.05 M_⊙ (Lee et al. 2018). Pineda et al. (2019) show that HH211 is in fact not in free fall, but is rotating. However, as we are working with a 1D model of the source, the freefall velocity profile is somewhat mimicing the rotation by the steep increase of velocity towards the centre (see Pineda et al. in prep).

Appendix E Molecular abundance profiles of the best-fit modelling results

Fig. E.1 and E.2 show the corresponding molecular abundance profiles of the best-fit modelling results derived with LOC and presented in this work, for L1544 and HH211, respectively.

Fig. E.1 Fractional molecular abundance (with respect to H2) profiles of the best-fit results produced with LOC for the spectra observed towards L1544. In the case of the 13C and 18O isotopologues, the abundance profiles correspond to the profiles of the main species, scaled down by the isotopic ratio (68 and 557, respectively).

Fig. E.2 Fractional molecular abundance (with respect to H2) profiles of the best-fit results produced with LOC for the spectra observed towards HH211. In the case of the 13C and 18O isotopologues, the abundance profiles correspond to the profiles of the main species, scaled down by the isotopic ratio (68 and 557, respectively).

References

All Tables

All Figures

	Fig. 1 Selection of the spectra observed towards the pre-stellar core L1544 (black) and the protostellar core HH211 (red). The dotted vertical lines indicate the respective rest velocities of the cores. For HCN, the F = 2–1 transition is shown. For CCH, C13CH CCD H13CN and DCN, the transitions listed in Table 1 are shown.
In the text

	Fig. 2 Central hyperfine (F=2–1) of the (1–0) transitions of HCN (black) and HC3N (orange), observed towards L1544. The data for HC3N are taken from Bianchi et al. (2023).
In the text

	Fig. 3 Multi-component Gaussian fit of the four effective hyperfine components found by van der Tak et al. (2009) to the spectrum of HN13Cbserved towards L1544.
In the text

	Fig. 4 Profiles of the gas temperature (red), H2 number density (orange, in logarithmic scale) and infall velocity (blue, in units of 0.1 kms−1) for the Keto-Caselli model (dashed, KC Keto et al. 2015), the HDCRT model (dotdashed, HD Sipilä et al. 2022) of L1544, and the physical model of HH211 (dotted) derived in this work.
In the text

Fig. 5 Synthetic spectra of CCH obtained with LOC (top rows: L1544, bottom row: HH211). Each hyperfine line is plotted separately, numbered from low to high frequency. For L1544, the results are splitted between using the Keto-Caselli model (top) and the HDCRT model (centre) as approximation for the physical structure of the core. Towards L1544, the best-fitting abundances and timesteps are 2 × 10−8 (KC), 1e5 yr (STKC), 3 × 10−8 (HD), 1e5 yr (STHD). Towards HH211, the best-fitting abundances and timesteps are 2 × 10−9 (const), 1e6 yr (2-phase).

In the text

Fig. 6 Synthetic spectra of HCN (left) and H13CN (right) obtained with LOC (top rows: L1544, bottom row: HH211). For L1544, the results are splitted between using the Keto-Caselli model (top) and the HDCRT model (centre) as approximation for the physical structure of the core. Towards L1544, the best-fitting abundances and timesteps are 2 × 10−8 (KC), 7.9e5 yr (STKC), 1 × 10−8 (HD), 1.3e5 yr (STHD) for HCN; 2 × 10−10 (KC), 6e5 yr (STKC), 3 × 10−10 (HD), 1e5 yr (STHD) for H13CN. Towards HH211, the best-fitting abundances and timesteps are 2 × 10−9 (const), 1e5 yr (2-phase) for HCN; 3 × 10−11 (const), 1.3e5 yr (2-phase) for H13CN.

In the text

	Fig. 7 Synthetic spectra of DCN obtained with LOC (top rows: L1544, bottom row: HH211). Towards L1544, the best-fitting abundances and timesteps are 1 × 10−9 (KC), 1.3e5 yr (STKC), 1 × 10−9 (HD), 9.8e5 yr (STHD) for DCN. Towards HH211, the best-fitting abundances and timesteps are 4 × 10−10 (const), 1e5 yr (2-phase) for DCN.
In the text

Fig. 8 Synthetic spectra of HNC and isotopologues obtained with LOC (top rows: L1544, bottom row: HH211). For L1544, the results are splitted between using the Keto-Caselli model (top) and the HDCRT model (centre) as approximation for the physical structure of the core. Towards L1544, the best-fitting abundances and timesteps are 2 × 10−8 (KC), 1e5 yr (STKC), 5 × 10−8 (HD), 1e5 yr (STHD) for HNC; 2 × 10−10 (KC), 1.3e5 yr (STKC), 2 × 10−10 (HD), 9.8e5 yr (STHD) for HN13C; 3 × 10−9 (KC), 1e5 yr (STKC), 8 × 10−10 (HD), 9.8e5 yr (STHD) for DNC. Towards HH211, the best-fitting abundances and timesteps are 5 × 10−10 (const), 1e5 yr (2-phase) for HNC 1 × 10−11 (const), 1e5 yr (2-phase) for HN13C; 7 × 10−11 (const), 1e6 yr (2-phase) for DNC.

In the text

Fig. 9 Synthetic spectra of HCO+ and isotopologues obtained with LOC (top rows: L1544, bottom row: HH211). For L1544, the results are splitted between using the Keto-Caselli model (top) and the HDCRT model (centre) as approximation for the physical structure of the core. Towards L1544, the best-fitting abundances and timesteps are 7 × 10−10 (KC), 9.8e5 yr (STKC), 3 × 10−9 (HD), 9.8e5 yr (STHD) for HCO+; 2 × 10−10 (KC), 1e5 yr (STKC), 4 × 10−10 (HD), 1e5 yr (STHD) for H13CO+; 6 × 10−12 (KC), 1e5 yr (STKC), 1 × 10−11 (HD), 1e5 yr (STHD) for HC18O+ ; 2 × 10−10 (KC), 1.7e5 yr (STKC), 4 × 10−10 (HD), 1.7e5 yr (STHD) for DCO+. Towards HH211, the best-fitting abundances and timesteps are 1 × 10−11 (const), 1e6 yr (2-phase) for HCO+; 7× 10−12 (const), 1e5 yr (2-phase) for H13CO+; 4× 10−13 (const), 1e5 yr (2-phase) for HC18O+; 2× 10−11 (const), 1e6 yr (2-phase) for DCO+.

In the text

	Fig. 10 C13CH/13CCH ratio of starless/pre-stellar (black) and protostellar (red) cores. L1521B: Taniguchi et al. (2019); L134N: Taniguchi et al. (2019); TMC-1: Sakai et al. (2010); L1527: Yoshida et al. (2019).
In the text

	Fig. A.1 Ground-state rotational lines observed towards the pre-stellar core L1544. The respective hyperfine transition is given in the upper left corner of each plot. The dotted line indicates the rest velocity of the system (7.2 km s−1 ). The red curves show the fits used in the LTE analysis.
In the text

	Fig. A.2 Ground-state rotational lines observed towards the protostellar core HH211. The respective hyperfine transition is given in the upper left corner of each plot. The dotted line indicates the rest velocity of the system (9.2 km s−1). The red curves show the fits used in the LTE analysis.
In the text

	Fig. C.1 Temperature dependence of hyperfine resolved rate coefficients of HCN and DCN (upper panels) and HNC, DNC and HN13C (lower panels) for the ( j, F) = (2, 1) → (1, F′) transitions. The data of HCN are from Magalhães et al. (2018).
In the text

Fig. D.1 Analysis of Herschel SPIRE data towards HH211. Top: H2 column density map of HH211, derived using the Herschel SPIRE images at 250, 350 and 500 µm. The red and blue arrows indicate the direction and size of the molecular outflow/jet from the protostellar core (marked by a black triangle). The Herschel beam size is shown as a white circle, the average beam size of our observations is shown as red circle. Bottom: Dust temperature map of HH211, derived using the Herschel SPIRE images at 250, 350 and 500 µm.

In the text

Fig. D.2 Analysis of the H2 column density map and the temperature map derived from Herschel data. Top: Radial H2 column density profile (red), overlaid with the fitted Plummer-like profile (blue). Given in green are the data points averaged over shells of 11 arcsec, the size corresponding to 1/3 of the Herschel beam size. Bottom: Radial temperature profile (red), overlaid with the best fit (blue) and the averaged shells (green).

In the text

	Fig. E.1 Fractional molecular abundance (with respect to H2) profiles of the best-fit results produced with LOC for the spectra observed towards L1544. In the case of the 13C and 18O isotopologues, the abundance profiles correspond to the profiles of the main species, scaled down by the isotopic ratio (68 and 557, respectively).
In the text

	Fig. E.2 Fractional molecular abundance (with respect to H2) profiles of the best-fit results produced with LOC for the spectra observed towards HH211. In the case of the 13C and 18O isotopologues, the abundance profiles correspond to the profiles of the main species, scaled down by the isotopic ratio (68 and 557, respectively).
In the text