High sensitivity molecular line observations towards the L1544 pre-stellar core challenge current models
J. Ferrer Asensio, S. S. Jensen, S. Spezzano, P. Caselli, F. O. Alves, O. Sipil\"a, E. Redaelli
TL;DR
High-sensitivity molecular line observations of the L1544 pre-stellar core reveal discrepancies with existing models, highlighting the need for improved chemical networks and better understanding of core dynamics.
Contribution
This study provides detailed radiative transfer modeling of multiple molecular lines towards L1544, challenging current chemical models and emphasizing the importance of S-depletion processes.
Findings
Different molecular transitions trace distinct gas dynamics.
Current chemical models poorly reproduce sulphur chemistry.
Single transition modeling is insufficient for accurate abundance profiles.
Abstract
The increased sensitivity and spectral resolution of observed spectra towards the pre-stellar core L1544 are challenging the current physical and chemical models. With the aim of further constraining the structure of L1544 as well as assessing the completeness of chemical networks, we turn to radiative transfer modelling of observed molecular lines towards this source. We obtained high-sensitivity and high-spectral resolution observations of HCO+ (J = 1 - 0), CS (J = 2 - 1), C34S (J = 2 - 1), H2CO (J ,Ka,Kc = 2,1,2 - 1,1,1), c-C3H2 (J,Ka,Kc = 2,1,2 - 1,0,1), SO (N,J = 2,3 - 1,2) and 34SO (N,J = 2,3 - 1,2) with the IRAM 30m telescope towards the dust peak of L1544. A non-Local Thermodynamic Equilibrium radiative transfer code is coupled to the Markov Chain Monte Carlo method to model the observations. We find that with just one transition for each isotope, the modelling cannot find a…
| Molecule | Transition | Frequency | n @ 10 K | rms | HPBW | Main Beam efficiency | Velocity resolution | ||
|---|---|---|---|---|---|---|---|---|---|
| (MHz) | (K) | (x 10-5 s-1) | (cm-3) | (mK) | (”) | (%) | (km s-1) | ||
| \ceHCO+ | = 1 - 0 | 89188.525 | 4.28 | 4.19 | 1.7e5 | 64 | 27.6 | 80 | 0.033 |
| HC17O+ | = 1 - 0 | 87057.53 | 4.18 | 3.9 | 1.5e5 | 2.8 | 28.3 | 81 | 0.034 |
| \ceCS | = 2 - 1 | 97980.95 | 7.05 | 1.68 | 3.6e5 | 44 | 25.1 | 79 | 0.029 |
| C34S | = 2 - 1 | 96412.94 | 6.94 | 1.60 | 3.6e5 | 16 | 25.5 | 79 | 0.030 |
| \ceH2CO | = 21,2 - 11,1 | 140839.50 | 21.92 | 5.30 | 6.8e5 | 61 | 17.5 | 74 | 0.020 |
| c-C3H2 | = 21,2 - 10,1 | 85338.89 | 6.44 | 2.32 | 1.7e5 | 52 | 28.8 | 82 | 0.034 |
| SO | = 2,3 - 1,2 | 99299.87 | 9.22 | 1.12 | 1.1e5 | 16 | 24.7 | 79 | 0.029 |
| 34SO | = 2,3 - 1,2 | 97715.39b | 9.09 | 1.07 | 1.1e5 | 19 | 25.1 | 79 | 0.030 |
| Transition | ain | aout | fv | ||
|---|---|---|---|---|---|
| (au) | (km s-1) | ||||
| C34S (2 - 1) | 1.8410-12 | 1.8210-9 | 7697.42 | 0.93 | 0.10 |
| 34SO (2,3 - 1,2) | 1.1610-13 | 3.7210-9 | 14075.36 | 1.06 | 0.11 |
| HC17O+ (1 - 0) | 1.9010-13 | 4.8710-11 | 9172.02 | 1.45 | 0.07 |
| CS (2 - 1) | 2.2110-10 | 1.1510-7 | 15827.86 | 1.99 | 0.01 |
| SO (2,3 - 1,2) | 4.5710-13 | 8.9310-9 | 8013.68 | 0.83 | 0.11 |
| c-C3H2 (21,2 - 10,1) | 1.1610-9 | 2.3310-9 | 19408.84 | 0.56 | 0.12 |
| c-C3H2 (21,2 - 10,1) extended | 3.6910-12 | 1.5610-9 | 1161.61 | 0.79 | 0.11 |
| HCO+ (1 - 0) extended | 7.8510-13 | 2.0010-10 | 884.21 | 1.64 | 0.11 |
| H2CO (21,2 - 11,1) | 1.1310-10 | 7.1610-8 | 19876.12 | 1.43 | 0.10 |
| H2CO (21,2 - 11,1) extended | 3.4110-14 | 7.8710-9 | 5219.93 | 1.11 | 0.09 |
| Model | ain | aout | fv | ||
|---|---|---|---|---|---|
| (au) | (km s-1) | ||||
| CS + C34S (2 - 1) | 6.9210-9 | 1.4410-7 | 9550.04 | 1.43 | 0.00 |
| SO + 34SO (2,3 - 1,2) | 1.0710-11 | 1.7110-8 | 5567.16 | 0.82 | 0.11 |
| HCO+ + HC17O+ (1 - 0) extended | 8.4410-11 | 9.8910-9 | 9428.14 | 0.87 | 0.11 |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsStellar, planetary, and galactic studies · Astro and Planetary Science · Astrophysics and Star Formation Studies
High sensitivity molecular line observations towards the L1544 pre-stellar core challenge current models 111Based on observations carried out with the IRAM 30m Telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain)
J. Ferrer Asensio
Centre for Astrochemical Studies, Max-Planck-Institut für extraterrestrische Physik, Giessenbachstr. 1, 85748 Garching, Germany
RIKEN Pioneering Research Institute, Wako-shi, Saitama, 351-0106, Japan
S. S. Jensen
Centre for Astrochemical Studies, Max-Planck-Institut für extraterrestrische Physik, Giessenbachstr. 1, 85748 Garching, Germany
S. Spezzano
Centre for Astrochemical Studies, Max-Planck-Institut für extraterrestrische Physik, Giessenbachstr. 1, 85748 Garching, Germany
P. Caselli
Centre for Astrochemical Studies, Max-Planck-Institut für extraterrestrische Physik, Giessenbachstr. 1, 85748 Garching, Germany
F. O. Alves
Centre for Astrochemical Studies, Max-Planck-Institut für extraterrestrische Physik, Giessenbachstr. 1, 85748 Garching, Germany
O. Sipilä
Centre for Astrochemical Studies, Max-Planck-Institut für extraterrestrische Physik, Giessenbachstr. 1, 85748 Garching, Germany
E. Redaelli
Centre for Astrochemical Studies, Max-Planck-Institut für extraterrestrische Physik, Giessenbachstr. 1, 85748 Garching, Germany
European Southern Observatory, Karl-Schwarzschild-Straße 2, 85748 Garching, Germany
Abstract
The increased sensitivity and spectral resolution of observed spectra towards the pre-stellar core L1544 are challenging the current physical and chemical models. With the aim of further constraining the structure of L1544 as well as assessing the completeness of chemical networks, we turn to radiative transfer modelling of observed molecular lines towards this source. We obtained high-sensitivity and high-spectral resolution observations of \ceHCO+ (J = 1 - 0), CS (J = 2 - 1), C34S (J = 2 - 1), \ceH2CO (J = 21,2 - 11,1), c-C3H2 (J = 21,2 - 10,1), SO (N,J = 2,3 - 1,2) and 34SO (N,J = 2,3 - 1,2) with the IRAM 30m telescope towards the dust peak of L1544. A non-Local Thermodynamic Equilibrium radiative transfer code is coupled to the Markov Chain Monte Carlo method to model the observations. We find that with just one transition for each isotope, the modelling cannot find a global minimum that fits the observations. The derived fractional abundance profiles are compared to those computed with chemical models. The observed transitions trace gas components with distinct dynamics at different distances along the radius of the core. Moreover, the results evidence a poor reproduction of sulphur chemistry by chemical modelling and stresses the need to include a more consistent S-depletion process to accurately reproduce the S-chemistry towards dense cores.
ISM: molecules - ISM: clouds - radio lines: ISM - stars: formation - radiative transfer
1 Introduction
Pre-stellar cores are high-density regions formed from the fragmentation of filamentary structures that constitute molecular clouds (Palmeirim et al., 2013; Hacar et al., 2022). While pre-stellar cores are in the early stages of gravitational contraction, they have not yet formed a protostar. These objects are characterised by high densities ( cm*-3*) and low temperatures (T K) at their centre (Keto & Caselli, 2008). The study of the physical, chemical and dynamical properties of these cores and their surroundings helps us gain a comprehensive understanding of the initial conditions that lead to core collapse and, therefore, to understand the earliest stages of the star formation process (André et al., 2014).
Molecular spectra have proven to be valuable tools for characterizing dense cloud cores. Molecular line intensities, widths, and profiles are influenced by the physical , chemical and dynamical conditions along the path of the emitted radiation. The information embedded within observed spectral lines is extracted through a comprehensive radiative transfer modelling of the emission accounting for the physical structure of the observed source, as well as the molecular abundance distribution computed with an astrochemical code (e.g., Caselli et al. 2002a; Sohn et al. 2007; Keto et al. 2015; Redaelli et al. 2019; Ferrer Asensio et al. 2022). Upon the fitting of the observed molecular transitions with a specific chemical model we can additionally gain information on the source’s evolutionary status.
L1544 is a well-known pre-stellar core in the Taurus Molecular Cloud at a distance of 170 pc (Galli et al., 2019). Extensive research has focused on the study of its physical and chemical properties, revealing it is centrally concentrated, with low central temperatures and signs of contraction motions (Ward-Thompson et al., 1999; Crapsi et al., 2005, 2007).
With access to increasingly high-sensitivity and high-spectral resolution spectra of molecules, we are challenging the current models and knowledge on L1544, learning about its sub-structure and the influence of its surroundings. For example, the analysis of high-spectral resolution observed spectra of CO, H2O, N2H*+, HCO+* and isotopologues towards this source revealed important information on the structure and kinematics of L1544 including the CO depletion towards the centre of the core and the the ionization fraction across the core (Caselli et al., 1999, 2002b). Moreover, past studies showed the need for adjusting the velocity profile of the physical model (Keto et al., 2015) to fit the profile of specific lines (Bizzocchi et al., 2013; Ferrer Asensio et al., 2022). In addition, spatially extended abundant molecules such as \ceHCO+ have shown the need to take into account a diffuse envelope around L1544 to be able to reproduce their molecular profiles (Redaelli et al., 2019). Moreover, Lin et al. (2022) found a local gas density enhancement towards L1544 deviating from the commonly assumed spherically symmetric Bonnor Ebert (BE) sphere structure, which could be the result of an accretion flow impact.
In this work, we deepen our comprehension of the physical, chemical and kinematic structure of L1544 complementing past works (Tafalla et al., 1998; Caselli et al., 1999, 2002a; Bizzocchi et al., 2013; Redaelli et al., 2019; Ferrer Asensio et al., 2022; Lin et al., 2022), through the radiative transfer modelling of new high-sensitivity and high-spectral resolution observations of molecular species towards this pre-stellar core. To probe the entirety of the core, we observe rare isotopologue transitions, C34S (J = 2 - 1) and 34SO (N,J = 2,3 - 1,2), expected to trace the whole core, and main isotopologue transitions, \ceHCO+ (J = 1 - 0), CS (J = 2 - 1), \ceH2CO (J = 21,2 - 11,1), c-C3H2 (J = 21,2 - 10,1) and SO (N,J = 2,3 - 1,2), which are expected to trace the outer parts of the core as they are optically thick. It is important to note that this is the first time that such high-sensitivity and high-spectral resolution data of these molecular transitions are available for this source. To adjust the pre-stellar core model as a consequence of the new information provided by recent observations (Ferrer Asensio et al., 2022), we adopt a non-Local Thermodynamic Equilibrium (non-LTE) approach coupled with a Markov Chain Monte Carlo (MCMC) method. This allows us to explore the parameter space of selected physical and chemical model variables to provide the highest-probability model solutions to fit the observations. Here, we assess the molecular fractional abundance profiles computed with the pseudo-time dependent chemical model (pyRate) described in Sipilä et al. (2015) by comparing them to simple step-abundance profiles found to give the best fits of observed transitions towards L1544.
This manuscript is structured as follows: the observations are presented in Section 2, the radiative transfer methodology used to reproduce the observed molecular transitions is described in Section 3, the modelling results are introduced in Section 4, these results are interpreted and placed in the context of past works in Section 5 and we summarize our conclusions in Section 6. Additionally, we include further tests in the Appendix.
2 Observations
Single-pointing observations of the \ceHCO+ (J = 1 - 0), CS (J = 2 - 1), C34S (J = 2 - 1), \ceH2CO (J = 21,2 - 11,1), c-C3H2 (J = 21,2 - 10,1), SO (N,J = 2,3 - 1,2) and 34SO (N,J = 2,3 - 1,2) rotational transitions are obtained towards the dust peak of L1544 ( = 05h04m17s.21, = +25∘10*′42′′.8). These observations were carried out in August and October of 2021, using the IRAM 30m telescope at Pico Veleta, and took a total of 5.7 hours of integration time. The telescope pointing was frequently checked against the nearby bright quasars B0316+413, B0851+202, and B0439+360, depending on the elevation of the source. We used the E090 and E150 bands of the EMIR receiver. The observations were performed in frequency-switching mode. We used the VESPA backend to achieve a spectral resolution of 10 kHz. The spectral velocity resolution for each transition is reported in Table LABEL:obsdat. Both horizontal and vertical polarisations were observed simultaneously. The ranged between 65 and 173 K. The summary of the observed transitions in this project alongside HC17O+* presented in Ferrer Asensio et al. (2022) can be found in Table LABEL:obsdat.
The observational data was processed and then averaged with the Continuum and Line Analysis Single-dish Software (class), an application from the gildas222https://www.iram.fr/IRAMFR/GILDAS/ software (Pety, 2005). The resulting spectra, whose intensities have been converted to main beam temperature (see Table LABEL:obsdat), are presented in Figure 1.
Signal-to-noise ratio (S/N) peak values are 11 for 34SO, 13 for the least bright component of \ceH2CO, 26 for \ceHCO+, 28 for \ceC34S, 48 for the least bright component of \ceCS, 49 for the least bright component of c-C3H2, and 69 for the least bright component of \ceSO.
3 Radiative transfer
In pre-stellar cores, where the molecular hydrogen volume density () and gas temperature () can span between 106 cm*-3* and 6 K at the centre and 102 cm*-3* and 20 K at the edge, we cannot assume LTE conditions throughout the core. The LTE regime is a reasonable assumption when the critical density of the transition (= / for optically thin lines, where is the Einstein A coefficient and is the collisional coefficient at a specific temperature, and where ‘’ and ‘’ are the upper and lower levels, respectively), is substantially lower than the volume density of the emitting region. Thus, for high transitions ( 105 cm*-3*), LTE applies only in a reduced area of the pre-stellar core centre where the exceeds . In non-LTE the energy level populations deviate from the Boltzmann distribution and need to be calculated solving the statistical equilibrium equations, as we describe below.
3.1 The LOC radiative transfer code + MCMC
Radiative transfer modelling of molecular transitions towards L1544 is typically aided by a detailed 1D physical pre-stellar core model presented in Keto et al. (2015). This model has been successful at reproducing many molecular transitions observations towards the dust peak of L1544 (Keto & Rybicki, 2010; Keto & Caselli, 2010; Caselli et al., 2012, 2017, 2022; Bizzocchi et al., 2013; Redaelli et al., 2018, 2019, 2021, 2022; Ferrer Asensio et al., 2022). The radial profiles of the pre-stellar core physical model are presented in Figure 2. The molecular hydrogen volume density (solid black line) ranges from 8 cm*-3* at the centre to 1 cm*-3* at the edge of the core (0.32 pc). The gas temperature (blue dashed line) ranges from 6 K at the centre to 18 K at the edge. Finally, the gas velocity (orange dash and dotted line) ranges from -0.14 km s*-1* at the velocity peak ( 0.01 pc) to -0.01 km s*-1* at the edge of the core. The negative sign of the velocity indicates that the pre-stellar core model is contracting.
In previous works, the adjustment of physical and chemical model parameters, such as the fractional abundance or the velocity profile scaling, aided the fit and understanding of the nature of some molecular emission towards L1544 (Bizzocchi et al., 2013; Ferrer Asensio et al., 2022; Redaelli et al., 2022). In this work, we explore the parameter space of specific physical and chemical parameters when performing the radiative transfer modelling of the observations to ensure the best possible fit. With this aim in mind, we combine non-LTE radiative transfer modelling with the Markov Chain Monte Carlo method (MCMC) to calculate probability distributions of the physical and chemical parameters, which produce the closest modelled spectra with respect to the observed spectra.
The radiative transfer modelling of the observed molecular lines is carried out with the Line transfer with OpenCL (LOC) code (Juvela, 2020). As a first step, LOC solves the statistical equilibrium equations with the radiative and collisional coefficients of the molecule with an accelerated lambda iteration (ALI) (Rybicki & Hummer, 1991). Then, the radiative transfer is computed using a 1D model of the source, taking into account the molecular hydrogen volume density, the gas temperature, the velocity of the gas, , the turbulent velocity dispersion, , and the fractional abundance profiles of the desired molecule. The 1D radial profiles of the physical properties are taken from the quasi-statically contracting Bonnor Ebert sphere described in Keto et al. (2015).
Additionally to the pre-stellar core physical model, LOC requires information on the fractional abundance profile of the molecules. The fractional abundance profiles for the molecules can be simulated using chemical models. In the past, fractional abundance profiles computed with pyRate were used when performing radiative transfer modelling of L1544 observations. These abundance profiles are simulated using the pre-stellar core physical model described in Keto et al. (2015), which is separated into concentric shells where chemical simulations, based on a chemical network, are conducted. The physical structure of the core is fixed, and the chemistry evolves with time. As one of the objectives is to assess the accurateness of fractional abundance profiles computed with pyRate, we use simple step abundance profiles described in Section 3.2 for the purposes of comparison.
The Markov Chain Monte Carlo method, as the name suggests, combines the Markov Chain process with Monte Carlo simulations. The MCMC code EMCEE is used (Foreman-Mackey et al., 2013). The MCMC method involves randomly sampling high-dimensional probability distributions while considering the probabilistic dependence between samples using the Markov Chain. The density estimation for probability distributions is done taking into account both the observations and the model. The value of the parameters that have the lowest minimization with respect to the observations is the value with the highest likelihood. The is computed channel wise and then summed over all channels.
First of all, we define the parameters to be explored by the MCMC and a sensible range for their values. Then, initial values, from where the random walk starts, are defined. The MCMC, starting from the initial values, generates a sequence of models using random parameter values. Depending on the minimization of the resulting model with respect to the observations, the parameters used are accepted or rejected. If they are accepted, these parameter values are used for the next step on the modelling chain. If the parameters are rejected, the previous values in the chain are used. As a product, we obtain the probability distribution of the parameter values as constrained by the comparison of the produced spectral models with the observed spectra in the form of histograms. In the corner plots, the parameter histograms are arranged to give information on the correlation between the parameters in pairs.
The models presented in this manuscript were run with 40 walkers for total chains of 10000 steps with a burn-in of 500 steps. The results correspond to the median values. For a well-constrained model, the median values will converge to the model with the highest probability (lowest ). In such cases, the histogram in a corner plot will often resemble a Gaussian profile, depending on the underlying posterior distribution. In cases where the parameter cannot be properly constrained, the histograms are often flat profiles, indicating that the model cannot constrain the parameter, either due to the lack of information or poor assumptions in the model. In the cases where the variable is not constrained, as seen as a flat, instead of Gaussian, histogram in the corner plot, the stated values correspond to the median as the highest probability were not able to be determined. The resulting modelled spectra are convolved with the size of the IRAM 30m beam at each frequency to enable the comparison with the observed spectra. The spectrum simulated by LOC corresponds inherently to a pencil beam, and to account for the beam dilution, the spectrum of the 1D spherically symmetric model is first projected onto a 2D map of the sky. A 2D Gaussian profile, with a FWHM corresponding to the HPBW of the IRAM 30m telescope at the observed frequency, typically extending from –3FWHM to +3FWHM, is then convolved with the intensity distribution of this 2D model. A 1D spectrum is subsequently extracted toward the centre of the core in the convolved 2D map.
3.2 Models and Parameter Space
Past works required small adjustments of the velocity profile scaling to reproduce some high-sensitivity observations of molecular transitions towards L1544 (Bizzocchi et al., 2013; Ferrer Asensio et al., 2022). Thus, for the modelling in this manuscript, we set the velocity profile scaling, , as a parameter to be explored by the MCMC. Moreover, we include the turbulent velocity dispersion, , as another parameter, as we expect the turbulence to vary in different parts of the core and surroundings. Turbulence is known to decrease within dense cores (e.g. Fuller & Myers 1992; Goodman et al. 1998; Pineda 2010; Choudhury et al. 2021), likely due to the dumping of magneto-hydrodynamic waves (e.g. Pineda et al. 2021; Caselli et al. 2002a).
Some articles in the literature have also shown the need to scale the fractional abundance profiles computed from pseudo-time-dependent chemical models (Redaelli et al., 2019; Ferrer Asensio et al., 2022). In order to gain insight into the molecular radial distributions directly from observations, we use simple step abundance profiles for the modelling. This choice of abundance profile takes into account the almost total molecular depletion observed towards the centre of the core (Caselli et al., 2022), which has an impact on the observed molecular emission.
Depending on the nature of the observed molecular transition, we use either a ”non-extended” or an ”extended” model. The ”non-extended” model corresponds to a pre-stellar core physical and chemical model 1D profiles that extend to 0.32 pc as used for the majority of the radiative transfer modelling of molecules towards L1544 (Giers et al., 2022; Ferrer Asensio et al., 2022). This model is from Keto et al. (2015) but sets the and as parameters to be explored by the MCMC. The molecular abundance profile is constructed with two constant abundance values, one representing the abundance in the inner part of the core, , and another for the outer part of the core, . The separation of the inner abundance and outer abundance regions is defined by a radius, . Consequently, the abundance profile has two zones: between 0 and pc and between and 0.32 pc (red stripe in Figure 3).
The ground state spectra of some abundant molecules, such as HCO*+*, trace also less dense gas around the pre-stellar core (Redaelli et al., 2022). Thus, for molecules expected to be present in this surrounding extended structure, we adopt the “extended” model. The extended model consists of the non-extended model structure (0 - 0.32 pc) with an additional layer spanning from 0.32 to 1 pc. The physical parameter profiles in this extended layer are assumed to be constant with the values given at the edge (0.32 pc) of the Keto et al. (2015) model. The parameter scales the velocity in the entire model from 0 to 1 pc. The abundance is kept constant in this external layer and is set to the observed value towards diffuse clouds. Thus, for the extended model, the abundance profile is defined by (between 0 and pc), (between and 0.32 pc), and diffuse cloud observed abundance (between 0.32 and 1 pc) (blue stripe in Figure 3). From this point on, the spectra resulting from the non-extended model are plotted with a solid red line, and the spectra resulting from the extended model are plotted with a dash-dotted blue line.
We have 5 parameters explored by the MCMC: , , , and . The prior probability distributions are set as:
- •
10*-15* 10*-6.5*
- •
10*-15* 10*-6.5*
- •
500 20000 (au)
- •
0 2
- •
0 0.35 (km/s).
These limits are set from previous knowledge of the physical, kinematic and chemical nature of the source. For example, scaling the velocity profile more than the limit specified above would be incompatible with previous observations and modelling of CO and H2O towards L1544 (Keto et al., 2015). We also impose the condition that cannot be larger than as we expect the molecules studied here to be depleted towards the centre of the core and, therefore, be less abundant in that area. The prior distribution is a uniform distribution, meaning that all values in the ranges mentioned above have the same prior probability.
4 Results
The modelling results for the observed molecular transitions are separated into two Sections: rarer isotopologue transitions and main isotopologue transitions (Sections 4.1 and 4.2, respectively). The LOC + MCMC parameter values for the models described in the following sections are summarised in Table LABEL:modres.
4.1 Rarer Isotopologue Transitions
Due to the intrinsically lower fractional abundances of the rare isotopologues, C34S, 34SO and HC17O*+* are expected to emit only from the core (i.e. within 0.32 pc). We do not expect their transitions to trace the extended structure, as traced with the HCO*+* = 1 - 0 transition (Redaelli et al., 2022). Thus, for the modelling of rarer isotopologue transitions presented in the next subsections, we use the non-extended model described in Section 3.2.
4.1.1 C34S
The C34S (2 - 1) transition presents an asymmetric line profile indicative of infall (Figure 1). To the best of our knowledge, this is the first time that the C34S (2 - 1) line has been seen to show clear self absorption. In Tafalla et al. (1998) this self absorption was hinted at (see Figure 3) but the S/N ratio was relatively low. Although the C34S (2 - 1) line has a profile typical of an optically thick line tracing a contracting core, we are including this transition in this Section, as it is a rarer isotopologue and its optical depth ( 4.22) is however significantly lower than the CS (2 - 1) line ( 105.1). For further information refer to Appendix A. The LOC + MCMC approach reproduces well the observations (observations in black, model in red; top panel in Figure 4). The LOC + MCMC parameter values and their uncertainties are listed in Table LABEL:modres. The resulting corner plot is presented in Figure 17.
The comparison of the fractional abundance profile used for the line profile from LOC + MCMC and the fractional abundance profile computed with pyRate is presented in the bottom panel of Figure 4. As pyRate does not take into account sulphur isotope chemistry, we scale down the CS fractional abundance profile by the [32S/34S] = 22 ratio (Wilson & Rood, 1994).
The radius where the C34S abundance drops due to depletion in the centre of the core found with the LOC + MCMC modelling is shifted to larger radii with respect to the radius computed with chemical models. The C34S value found with the LOC + MCMC modelling is a factor 9 larger than the highest C34S abundance from pyRate. The high abundances derived with the LOC + MCMC are connected to the high optical depth derived for this species.
In order to see the effect of the abundance profile used on the output line profile, we compare the results of the LOC + MCMC modelling approach, which we will from now on address as Model 0 (red solid line, top panel Figure 4), with the spectrum computed with the abundance profile from pyRate and the resulting LOC + MCMC and values (Model 1; dotted grey line in the top panel of Figure 4). The Model 1 spectrum shows a single-peaked line profile and under-reproduces the observed line intensity by a factor of 6 with respect to the observations.
Moreover, to see the effect of the LOC + MCMC and values on the modelled line profile with respect to the line profile computed by using the original velocity profile , from the physical core model described in Keto et al. 2015 and commonly assumed = 0.075 km/s, we compute a third spectrum by using both the abundance profile from pyRate, the velocity profile from Keto et al. 2015, and by setting = 0.075 km/s (Model 2; dashed grey line in the top panel of Figure 4). This spectrum shows a line profile tending to a single-peak morphology and underestimates the line intensity with respect to the observations by a similar factor as Model 1.
4.1.2 34SO
The 34SO (2,3 - 1,2) transition has a single-peaked component with an asymmetric shoulder at 6.9 km/s (Figure 1). Note that this line profile is not related to infall asymmetry. The infall asymmetry profile is seen as two-peaked line profiles with a less-intense red component compared to the blue component. Here, the asymmetric shoulder appears towards lower VLSR (blue part). We looked for possible line candidates that fall close to the 34SO rest frequency without success. As of now the physical reason behind this asymmetric shoulder is not clear.The central velocity of this transition is shifted by 0.23 km/s with respect to the VLSR of the core when using the frequency value taken from the CDMS (97715.3170 0.0500 MHz). Thus, the 34SO (2,3 - 1,2) transition was remeasured with the CAS Absorption Cell (CASAC) in the laboratory of the Center for Astrochemical Studies (CAS) in Garching bei München, Germany. The CASAC setup provides higher frequency accuracy measurements compared to the previously reported measurements for this molecule in the literature. The rest frequency was measured to be 97715.395 0.023 MHz (Table LABEL:obsdat).
The resulting model fits the observed transition well, excluding the asymmetric shoulder (observations in black, Model 0 in red, top panel in Figure 5). Nevertheless, most of the variables are not constrained, making the corner plot deviate from Gaussian profiles to more flattened profiles (Figure 18). The resulting Model 0 parameters are presented in Table LABEL:modres. The parameter was not constrained, appearing as a flat profile in the corner plot. Moreover, the , and parameters do not show strong constraints. The , on the other hand, appears better constrained with a visible peak in the corner plot. These uncertainties can also be seen represented with a broad red band in the bottom panel of Figure 5. The asymmetric shoulder may be adding uncertainty to the modelled parameters, especially the ones directly linked with the line width.
The comparison of the fractional abundance profile used for the LOC + MCMC line profile modelling and the fractional abundance profile computed with pyRate is presented in the bottom panel of Figure 5.
As mentioned above, the radius where the 34SO abundance drops due to depletion in the centre of the core found with the LOC + MCMC modelling presents large errors (seen as a broad red vertical band in the bottom panel of Figure 5). The depletion radius is larger in the LOC + MCMC model than in the pyRate abundance profile, even when taking into account the error bars. The 34SO value found with the LOC + MCMC modelling is a factor of 3 higher than the highest 34SO abundance from pyRate. Nevertheless, the LOC + MCMC value and the pyRate abundance in this region agree within errors.
As done for C34S, we compare the spectra computed with the different abundance profiles and physical parameter combinations, Model 0, 1 and 2 (for more information refer to Section 4.1.1). Model 1 (dotted grey line in top panel, Figure 5) overestimates the intensity with respect to the spectrum computed with Model 0 by a factor of 1.4 (red solid line in top panel, Figure 5). Model 2 is double-peaked and overestimates the line intensity by a factor of 1.5 (dashed grey line in top panel, Figure 5).
As mentioned above, the spectra computed with the abundances predicted by pyRate, Models 1 and 2 (grey spectra in top panel, Figure 5), are more intense than the one computed with the step abundance profile constrained by the LOC + MCMC approach; Model 0 (red spectrum in top panel, Figure 5). This is due to the fact that the pyRate abundance profile is significantly below that derived by the LOC + MCMC (dash-dotted grey and white lines, respectively, in bottom panel, Figure 5). The line intensity results from a combination of factors including the of the emitting region. The spectra computed with the abundances predicted by pyRate, shown in grey in top panel, Figure 5, arise from a region closer to the centre of the core (see the peak of the abundance profile; dash-dotted grey line in bottom panel, Figure 5), with respect to the LOC + MCMC computed spectrum (red in top panel, Figure 5), which emits from larger radii (white line in bottom panel, Figure 5).
4.1.3 HC17O*+*
The HC17O*+* (1-0) transition, first presented in Ferrer Asensio et al. (2022), is split into three hyperfine components ( 5/2 5/2, 5/2 7/2 and 5/2 3/2, from lower to higher velocity), which are in turn displaying a double-peaked line profile (Figure 1). The hyperfine structure of the HC17O*+* (1-0) transition occurs because the nuclear spin of 17O (I = 5/2) and the molecular rotation are coupled, resulting in the splitting of the J = 1 rotational level. In Ferrer Asensio et al. (2022), the observations were successfully reproduced with LOC using the pre-stellar core physical model described in Keto et al. (2015) and the HC17O*+* fractional abundance profiled scaled from the HCO*+* fractional abundance profile computed with pyRate. In this work, we test the LOC + MCMC approach using a step abundance profile (solid white line over a red solid band, Figure 6).
The LOC + MCMC fits the line double-peaked profile slightly overestimating the intensity at the centre, not fully reproducing the depth of the dip, nevertheless, this discrepancy falls within the noise (Model 0, Figure 6). The variable values are presented in Table LABEL:modres and the corner plot in Figure 19.
The comparison of the LOC + MCMC fractional abundance profile and the fractional abundance profile computed with pyRate is presented in the bottom panel in Figure 6. As pyRate does not take into account oxygen isotope chemistry, we scaled down the HCO*+* fractional abundance profile by the [16O/17O] = 2044 ratio (Penzias, 1981; Wilson & Rood, 1994). The HCO*+* abundance profile used in this comparison is an updated abundance profile with respect to Ferrer Asensio et al. (2022), as this better fits the HCO*+* observations towards L1544 (Redaelli et al., 2022).
The radius where the HC17O*+* abundance drops due to depletion in the centre of the core found with the LOC + MCMC modelling agrees with the radius at which the fractional abundance profile computed with chemical models presents a drop. Nevertheless, the abundance drop towards the core centre, as computed with the chemical models, is less steep compared to other molecules, and may not be well reproduced with the simplified step-abundance profile. The HC17O*+* value found with the LOC + MCMC modelling is a factor 5 larger than the highest HC17O*+* pyRate abundance. Nevertheless, the LOC + MCMC value and the pyRate abundance in this region agree within errors.
A comparison between the simulated line profiles from the three models is shown in Figure 6. Model 1 (dotted grey line) underestimates the intensity by a factor of 3 and gives double-peaked hyperfine components. Model 2 is also less intense by a similar factor, and the hyperfine components appear single-peaked (dashed grey line in the top panel in Figure 6).
4.2 Main Isotopologue Transitions
The molecules included in the next subsections are abundant enough to be potentially present in the extended structure beyond 0.32 pc. Our approach is first to model the main isotopologue transitions arising from abundant molecules with the non-extended model and then test the extended model if the previous fails to reproduce the observations.
4.2.1 CS
The CS (2 - 1) transition is double-peaked with a dip that almost reaches the baseline (Figure 1). This transition was originally presented in Tafalla et al. (1998). The blue and red peaks have similar intensities. This transition was fit with the non-extended model described in Section 3.2, and the results are presented in the red spectrum in the top panel in Figure 7.
The model fits well the observed line profile (black spectrum in the top panel in Figure 7). Moreover, the variables are highly constrained compared to the ones found for the C34S line (see corner plot in Figure 20 and abundance profile on the bottom panel in Figure 7). The resulting LOC + MCMC parameters are presented in Table LABEL:modres. The , and parameters in particular present low errors, which translate to sharp columns in the corner plot (see Table LABEL:modres and Figure 20). Moreover, the and values were found to be 1.99 and 0.01 , respectively which are close to the boundaries (2 and 0) of the used prior distribution.
The comparison of the fractional abundance profile used for the non-extended model and the fractional abundance profile computed with pyRate is presented in the bottom panel in Figure 7.
The radius where the CS abundance drops due to depletion in the centre of the core found with the non-extended LOC + MCMC modelling appears shifted towards larger radii with respect to the chemical models. The CS value found with the LOC + MCMC modelling is a factor of 20 higher than the largest abundance from the chemical model.
We compare the spectra computed with Model 0 with Models 1 and 2 in Figure 7. Model 1 (dotted grey line, top panel in Figure 7) underestimates the line strength by a factor of 1.4 but reproduces well the self-absorption dip. Model 2 (dashed grey line, top panel in Figure 7) also underproduces the line intensity by a factor of 1.4, but less than for the model using the resulting LOC + MCMC and values.
As the non-extended model seems to be sufficient to reproduce the self-absorption dip, which is expected to be caused by the outermost layers of the core, we do not model CS (2 - 1) with the extended model. With the idea of better constraining the CS model parameters, we perform a combined modelling of the CS and C34S (2 - 1) transitions. This test does not bring additional constraints with respect to the CS and C34S separated models. For more information, see section C.1 in the Appendix.
4.2.2 SO
The SO (2,3 - 1,2) transition presents the characteristic infall asymmetry line profile (Figure 1). Unlike other optically thick lines introduced in this work (CS, HCO*+* and H2CO), SO presents a shallow self-absorption feature. The SO (2,3 - 1,2) transition is fitted with the non-extended model. As we can see in the top panel in Figure 8, Model 0 (in red) fits well the observations (in black). The model parameters are presented in Table LABEL:modres and in the corner plot in Figure 21.
The comparison of the fractional abundance profile used for the LOC + MCMC and the fractional abundance profile computed with pyRate are presented in the bottom panel in Figure 8.
The radius where the SO abundance drops due to depletion in the centre of the core found with the LOC + MCMC modelling is shifted to larger radii with respect to the radius of the fractional abundance profile computed with chemical models. The SO value found with the LOC + MCMC modelling is a factor of 2 lower compared to the abundance of SO in the chemical model.
The comparison of the spectra computed with Model 0 with Models 1 and 2 is presented in the top panel in Figure 8. Model 1 results in a single-peaked line and overestimates the line intensity by a factor of 2. Model 2 (dashed grey line in top panel in Figure 8) also overestimates the line intensity by a factor of 1.7.
Furthermore, as for the case of CS, we perform a combined fit of SO and 34SO to try to better constrain the parameters found with the individual models of the molecules. This test does not add additional constraints (see Section C.2 in the Appendix).
4.2.3 HCO*+*
The HCO*+* (1 - 0) line shows a pronounced self-absorption dip that reaches the spectral baseline (third-column, first-row panel in Figure 1). Moreover, there is a blue excess feature centred at 6.9 km/s. With the LOC RT model utilized here, we reproduce the HCO*+* (1-0) fit performed by Redaelli et al. (2022) using the RT code MOLLIE (Keto & Field, 2005). This is shown in section C.2.1 of the Appendix.
In this Section, we test the LOC + MCMC methodology to fit the HCO*+* observations. As previous works pointed out the need for an external layer (Redaelli et al., 2019, 2022), we modelled HCO*+* with the extended model described in Section 3.2 with some adjustments. The velocity and were fixed to -0.05 km/s and 27 cm*-3*, respectively (Redaelli et al., 2022). The abundance of HCO*+* in the external layer is kept constant with the value of the abundance profile predicted by pyRate at 0.32 pc, as done in Redaelli et al. (2022).
This approach reproduces the HCO*+* (1 - 0) fairly well, except for the separation between the two peaks, which is overestimated by the model (top panel, Figure 9). The model parameters are presented in Table LABEL:modres as well as in the corner plot in Figure 22. The parameters are fairly well constrained with the exception of , as expected from an optically thick transition, which does not trace the core inner regions.
The comparison of the fractional abundance profile used for the LOC + MCMC and the fractional abundance profile computed with pyRate are presented in the bottom panel in Figure 9. The value agrees within errors with the pyRate abundance profile drop. The parameter found is of the same order of magnitude as the pyRate one.
In the top panel of Figure 9 the spectra computed with Models 0, 1 and 2 are presented. Both Models 1 and 2 overestimate the line intensity by a factor of 3 and do not reproduce the depth of the absorption feature. Moreover, Model 1 results in a broader line profile than the one observed.
In addition, we also tested fitting HCO*+* and HC17O*+* together. This test does not add additional information than the one found for the single HCO*+* and H17CO*+* model fits (see Section C.2.2 in the Appendix).
4.2.4 \ceH2CO
The H2CO (21,2 - 11,1) transition presents a self-absorption dip that reaches the baseline (fourth-column, first-row panel in Figure 1). Moreover, the blue component is brighter than the red component. The asymmetry of the blue and red peaks indicates that H2CO traces contracting gas. The results for the non-extended modelling are presented in red in Figure 10. The LOC + MCMC modelling reproduces well the observations. The corner plot presents for all of the parameters, except for , two peaks indicating two possible solutions (Figure 23). The model parameter values are listed in Table LABEL:modres.
The comparison of the fractional abundance profile used for the LOC + MCMC and the fractional abundance profile computed with pyRate are presented in the bottom panel in Figure 10. The value is shifted to larger radii with respect to the drop in abundance in the pyRate abundance profile. The parameter value found agrees within errors with the highest value of the pyRate abundance profile.
As done with the other molecules, we compare the spectra computed with Model 0 with Models 1 and 2 (dotted grey line in top panel in Figure 10). Both Models 1 and 2 overestimate the line intensity by a factor of 2.5, do not reproduce the depth of the absorption feature and have broader line profiles than the observed one.
In view of these results, we test the extended model (see Figure 11). The constant H2CO abundance profile in the external layer is fixed to the value 3.710*-9* as found for diffuse clouds (Snow & McCall, 2006) The extended model over-predicts the line intensity and the separation of the two peaks (blue dash-dotted line, top panel, Figure 11). Nevertheless, this model presents an improvement with respect to the non-extended model when it comes to reproducing the self-absorption dip. This result implies that H2CO needs to be present in the diffuse surrounding cloud to reproduce the depth of the self-absorption feature, as seen for HCO*+* (1 - 0) (Redaelli et al., 2022) The corner plot for the extended model (Figure 24) shows a larger scatter
The comparison of the fractional abundance profile used for the LOC + MCMC and the fractional abundance profile computed with pyRate are presented in the bottom panel in Figure 11. The pyRate abundance profile extends to 0.32 pc only. Moreover, as the is a fixed value only, a sky blue line is plotted without uncertainties. The value agrees within error with the radii of the pyRate abundance profile. The parameter value is a factor of 4 lower than the highest value of the pyRate abundance profile.
The resulting spectra computed with Models 1 and 2 are presented in Figure 11). Both Models 1 and 2 overestimate the line intensity by a factor of 1.8 and do not reproduce the depth of the absorption feature and have broader line profiles than observed
4.2.5 c-C3H2
The c-C3H2 (21,2 - 10,1) transition also presents a self-absorption feature but not as pronounced as for CS, HCO*+* and H2CO (fourth-column, second-row in Figure 1). The blue and red peaks show similar intensities, hinting at this transition arising from a static layer, which was previously deduced from the same transition in Tafalla et al. (1998). This is also seen for CS (2 - 1) (Section 4.2.1).
The results of the non-extended model are shown with a solid red line in the top panel in Figure 12. The model reproduces well the intensity of the blue peak but underestimates the intensity of the red peak by a factor of 1.5. The parameters are presented in Table LABEL:modres as well as in the corner plot in Figure 25.
The comparison of the fractional abundance profile used for the LOC + MCMC and the fractional abundance profile computed with pyRate are presented in the bottom panel in Figure 12. The value is shifted to larger radii with respect to the main drop in abundance in the pyRate abundance profile. The parameter value found is 3 times higher than the maximum value of the pyRate abundance profile.
The comparison of the spectra computed with Model 0 with Models 1 and 2 can be found in the top panel in Figure 12. Both Models 1 and 2 underestimate the line intensity by a factor of 2. Moreover, in Models 1 and 2 the double-peaked line profile is only hinted.
We also tested the extended model to see if this would improve the fit with respect to the non-extended model (see Figure 13). We set the c-C3H2 constant abundance value in the external layer to 6.410*-10*, as seen in diffuse clouds (Snow & McCall, 2006). The fits of the extended and non-extended models are similar (see Figure 12). The parameters are presented alongside the corner plot in Figure 26. With the exception of the and , the parameters found for the non-extended and extended models are similar.
The comparison of the fractional abundance profile used for the LOC + MCMC and the fractional abundance profile computed with pyRate are presented in the bottom panel in Figure 13. The value partially agrees with the drop in the pyRate abundance profile. The parameter value found is 2 higher than the maximum value of the pyRate abundance profile.
Both Models 1 and 2 underestimate the line intensity by a factor of 2. In this case the double-peaked line profile of Model 2 becomes more apparent.
The c-C3H2 (21,2 - 10,1) transition is not constrained with the LOC + MCMC modelling approach. Additional transitions are needed to improve the fit and decrease the uncertainties in the model parameters.
5 Discussion
In general the LOC + MCMC modelling is successful in reproducing the observed spectra. Nevertheless, not all of the parameters are constrained. Two types of physical models are tested in this manuscript; the extended and non-extended models. The non-extended model was used for all of the transitions except for the HCO*+* (1 - 0) one, which was directly fitted with the extended model as previous works showed the need to take into account an extended structure to fit the observed line profile (Redaelli et al., 2022). All of the observed spectra are reasonably well fitted with the non-extended model. As a test, H2CO (21,2 - 11,1) and c-C3H2 (21,2 - 10,1) were also modelled with the non-extended model. These lines were selected because their line profiles show self-absorption, suggesting they may trace an extended structure, and because their fits with the non-extended model showed room for improvement. The extended model fits improve the results from the non-extended model for the H2CO (21,2 - 11,1) case. Thus, both H2CO (21,2 - 11,1) and HCO*+* (1 - 0) seem to be tracing material outside 0.32 pc, commonly considered in radiative transfer modelling of L1544 as the edge of the cloud. Additionally, with the aim of constraining the parameters explored by the MCMC by giving more initial constraints, another approach was tested. This additional test consisted in fitting the main and rare isotopologues together. This does not improve the constraints on the modelled fractional abundance profiles (see Appendix).
The fact that the models fit well the observations visually but the parameters explored with the MCMC are not constrained indicates that there may be other possible parameter combinations that result in a good fit of the spectra. Thus, high-resolution spectra of single lines do not provide enough information to converge to one unique solution. A way to improve the constrain on the parameters would be to fit multiple lines for each species and isotopologue as done in Redaelli et al. (2019) and Spezzano et al. (2025). Future work should aim to include multiple lines of each species, to provide more constraints for the fitting (see, e.g., Lin et al. 2022; Jensen et al. 2024).
5.1 Evaluation of Rarer Isotopologue Fits
From the rarer isotopologue lines presented in Section 4.1 the one with the most constrained is HC17O*+* (1-0). This, alongside the higher value of found, which is interpreted as an enhanced gas contraction velocity towards the core centre (Ferrer Asensio et al. 2022), could indicate that HC17O*+* (1-0) is the line tracing the innermost part of the core out of the three rarer isotopologues, and is possibly the most optically thin. Another interesting result from the HC17O*+* (1-0) modelling is the highly unconstrained parameter. The inability for the model to constrain may come from the abundance profile decreasing in a more gradual fashion as predicted by pyRate (dash-dotted profile on the lower panel of Figure 6). This may also be the reason for the unconstrained , which depends on .
The unconstrained of C34S (2 - 1) seems to indicate this line may be slightly optically thick, as traced by the blue asymmetric line profile observed, and not tracing the innermost part of the core as seen for HC17O*+* (1-0). An additional test to determine the optical depth of this transition (see Appendi A) returned a value of 4.22, which is unexpectedly high. This is most likely due to the found for C34S which is almost an order of magnitude higher than the peak abundance computed by pyRate (see dash-dotted grey line in the bottom panel of Figure 4).
The most striking result from the 34SO (2,3 - 1,2) LOC + MCMC modelling is the significant displacement of with respect to the abundance profile computed with pyRate (see dash-dotted grey line in the bottom panel of Figure 5). Specifically, the radius where 34SO depletes computed with the LOC + MCMC lays at larger radii than the one shown in the fractional molecular abundance profile computed with chemical modelling. One has to keep in mind that was not well-constrained, but these results could also be pointing at a problem with the chemical network used to compute the SO, and consequently 34SO, fractional abundance profile.
From these results, C34S and 34SO seem to trace a region at similar radii in the inner part of L1544, while HC17O*+* traces deeper regions towards the centre of the core. These transitions provide constraints on the molecular fractional abundance profiles and on the kinematic structure of the inner parts of L1544.
5.2 Evaluation of Main Isotopologue Fits
From the main isotopologue transitions, CS (2 - 1) is the one with the most constrained parameters. Nevertheless, upon inspection of the parameter values found by the LOC + MCMC modelling, the physical correctness of this solution seems doubtful. First of all, the values for and , 1.99 and 0.00 respectively, match the bounds set in the prior distribution of 2 and 0 respectively. This seems to indicate the assumed model is too simplistic for a convergence to be found. Moreover, as seen for C34S, the value of is 2 orders of magnitude higher than the maximum abundance computed with pyRate (see dash-dotted grey line on the bottom panel of Figure 7). This seems to be the reason why the calculated CS optical depth is unexpectedly high ( 105.1, see Appendix A).
As seen for 34SO, the calculated by the LOC + MCMC approach of SO is also significantly shifted towards larger radii (see dash-dotted grey line in the bottom panel of Figure 8).
Observationally, SO has an extended morphology in some pre-stellar cores (e.g. in L183, Swade 1989; Dickens et al. 2000 and Pagani et al. 2005) and it could also be the case for L1544 (see Figure A.3 in Spezzano et al. 2017). The drop in the outer parts of the core in the SO fractional abundance profile computed with PyRate, shown in a grey dashed line in the bottom panel in Figure 8, can result from a combination of factors: photodissociation, ”depleted” initial abundance assumption and uncertainties in sulphur chemistry. The ”depleted” initial abundance assumption involves reducing the initial S abundance to avoid an excess of sulphur-bearing species in the inner part of the core to reproduce the observed sulphur depletion known as ”the sulphur depletion problem” (e.g. Ruffle et al. 1999). This results in an underestimation of the sulphur-bearing abundances on the outer parts of the core. At the edge of the core, where the visual extinction is estimated to be of the order of 1 mag, 61% of SO destruction is mostly due to reactions such as SO + C*+, with 38% of the SO destruction resulting from photodissociation. In regions with Av = 2 mag most of the destruction comes from SO + C+* and only 4% of the destruction is due to photodissociation. This seems to point out at a non-negligible effect of the initial S depletion abundance as well as possibly the incompleteness of the sulphur chemical network on the computed low SO abundance at the edge of the core.
The HCO*+* (1 - 0) transition is not that well reproduced compared to the other main isotopologue transitions. The complicated line profile seems to be affecting the line width of Model 0.
The H2CO (21,2 - 11,1) transition is fairly well fitted by Model 0. Nevertheless, the parameter exploration by the MCMC resulted into a double peaked distribution indicating that there are two equally possible solutions to fit this line profile. In view of these results we tested the extended model, which resulted in a better reproduction of the self-absorption feature. Nevertheless, the constraint of the parameters was not improved.
The c-C3H2 (21,2 - 10,1) transition is fairly well fitted by Model 0 except for the red peak whose intensity is under-reproduced. As done for H2CO we have tested the extended model to see if the fit would be improved, but in this case this resulted in a similar fit than the non-extended model.
5.3 Spatial Distribution of Molecular Emission in L1544
All of the transitions presented in this manuscript trace different parts of the L1544 core, which have different physical, chemical and kinematic properties. To illustrate this point, we sketch in Figure 14 the relative location of the studied molecular transitions within L1544. The structure in Figure 14 results from the interpretation of the observed line profiles. As not all of the LOC + MCMCs’ radii were well constrained, exact spatial distributions of the molecules cannot be determined. That is also the reason why precise modified gas velocity and turbulence profiles to fit all of the transitions together cannot be constructed. The main factors which decide the position of the layers are the presence or absence of self absorption and its magnitude and the double peaked line profile tracing contraction. The presence of a self-absorbed feature requires a foreground layer that absorbs the radiation emitted on a background layer. Thus, molecular transitions that show a self absorbed line profile trace spatially more extended regions than molecular transitions that do not show self absorption, being then present in external layers. The more material in the foreground layer, the more emission can be absorbed resulting in a deeper self-absorption feature.
In the outermost layer we are placing HCO*+* (1 - 0), as it is the only transition whose model requires the molecule to be present beyond the commonly assumed edge of the cloud (0.32 pc). The HCO*+* (1 - 0) line profile is characterised by a very deep self-absorption feature that reaches the baseline level, as well as a pronounced blue line asymmetry. This indicates that the molecule is abundant in a foreground layer capable of absorbing background emission and that this layer is undergoing contraction. The second transition with a comparably deep self-absorption feature is H2CO (21,2 - 11,1). The self-absorption in this transition also reaches the baseline, and like HCO*+*, the self-absorption dip is just fitted by including material beyond the core’s boundary. Thus, it is interpreted as tracing material beyond the dense core boundary. The blue asymmetry in its line profile similarly indicates contraction in the layer it traces. Following the self-absorption dip trend, the next transition is CS (2 - 1), which shows a deep but shallower self-absorption feature that does not reach the baseline. The comparable brightness of the blue and red peaks suggests that CS is tracing a static layer in the outer parts of the core.
The CS layer is followed by that of c-C3H2 (21,2 - 10,1), which also shows a self-absorption dip and similar blue and red peak brightness, indicating that this molecule is tracing a static middle-to-outer layer of the core. Next is SO (2,3 - 1,2), which displays the prototypical blue asymmetry line profile, suggesting contraction. However, the relatively weak self-absorption dip implies that this line originates from an intermediate region within the L1544 core. In addition, SO shows a significant shift in its depletion radius () compared to the predictions from chemical modelling, possibly pointing to incomplete sulphur chemistry or underestimated abundances in the outer layers.
The inner layers of the core are traced by the rarer isotopologues. Continuing inward, we find 34SO (2,3 - 1,2), which is the only single-peaked line in this study, suggesting that it is not abundant enough in the outer layers to produce self-absorption, nor is it tracing the very innermost, high-velocity region. It is therefore assigned to an intermediate position in the diagram. C34S (2 - 1) and HC17O*+* (1 - 0) probe the deepest layers of L1544. Both lines show double-peaked profiles that suggest they trace the contracting central regions, where freeze out is almost complete (Caselli et al., 2022). Among these, HC17O*+* shows the highest velocity profile scaling (), indicating it likely originates from the very centre of the core. It is thus placed in the innermost layer of the diagram, followed by C34S at slightly larger radii.
The comparison between the spectra computed with the LOC + MCMC approach, with the original fractional abundance profiles computed with pyRate and the resulting LOC + MCMC and values, and the original fractional abundance and velocity profiles and = 0.075 km/s, showed the need of modifying the abundance profiles and physical parameters used for the modelling of the high-sensitivity and high-spectral resolution observed lines towards L1544.
6 Conclusions
This work has presented new high-sensitivity and high-spectral resolution observations of the \ceHCO+ (J = 1 - 0), CS (J = 2 - 1), C34S (J = 2 - 1), \ceH2CO (J = 21,2 - 11,1), c-C3H2 (J = 21,2 - 10,1), SO (N,J = 2,3 - 1,2) and 34SO (N,J = 2,3 - 1,2) as well as HC17O*+* (J = 1 - 0) from Ferrer Asensio et al. (2022) rotational transitions towards the dust peak of the L1544 pre-stellar core.
The modelling successfully reproduces most line profiles, but not all parameter distributions were well constrained by the MCMC, highlighting degeneracies in the modelling approach. The only two transitions requiring an extended envelope beyond 0.32 pc for a good fit are HCO*+* (1–0), consistent with previous works, and H2CO (21,2 - 11,1). c-CH2 showed self-absorption but did not benefit from an extended model. The combined modelling of the main and rare isotopologues did not improve parameter constraints, suggesting that single-line fits lack the depth needed to yield unique solutions. Although the LOC + MCMC approach allows for flexible modelling of molecular abundances and kinematics, it remains sensitive to prior assumptions and observational limitations. To overcome these limitations, future efforts should focus on multi-line fitting for each isotopologue.
A layered picture of the L1544 core was developed based on the morphology of the observed line profiles and modelling outcomes. Transitions showing deep self-absorption and/or extended profiles (e.g., HCO*+, H2CO) trace the outer layers, while rarer isotopologues such as C34S and HC17O+* reveal the kinematics of the inner regions through their double-peaked profiles. Among these, HC17O*+* stands out as a probe of the innermost contracting layers, evidenced by its high velocity scaling and relatively constrained abundance drop-off radius.
A comparison of the fractional abundance profiles of the molecules obtained via LOC + MCMC models or by chemical modelling points to the problems of the chemical networks due to the uncertainties in elemental abundances input in chemical models, in particular sulphur. The poor reproduction of the sulphur chemistry hints at the need of the inclusion of a more consistent sulphur-depletion process, starting from the elemental cosmic abundance of S, to be able to model sulphur-bearing molecules towards pre-stellar cores.
The results presented in this manuscript show that all of the molecular transitions trace different parts of the L1544 core with different physical, chemical and kinematic properties, revealing the sub-structure of the L1544 pre-stellar core. Future 3D models will be considered to assess whether the deviation from spherical symmetry is needed. Moreover, these results point to the need to use 3D radiative transfer modelling to be able to interpret increasingly sensitive and spectrally resolved spectra towards L1544.
Acknowledgments. The authors gratefully acknowledge the support of the Max Planck Society. J. Ferrer Asensio thanks RIKEN Special Postdoctoral Researcher Program (Fellowships) for financial support. Moreover, the authors thank the anonymous referees for their constructive comments and helpful suggestions, which improved the clarity and quality of the manuscript.
Appendix A Optical Depth Estimations
Estimations of the optical depth of selected transitions performed with the LOC software are presented. The optical depths are computed at different values of VLSR resulting into an optical depth profile. Figure 15 shows the C34S (2 - 1) observed and LOC + MCMC modelled spectra (in black and red, respectively) alongside the optical depth profile (in dashed grey). The maximum optical depth of the line has a value of 4.22. The optical depth profile appears double peaked. Despite being from a rarer isotopologue, this line is sufficiently optically thick to present self absorption which, alongside the contraction of the core, results into the characteristic infall asymmetry profile. On the other hand, the optical depth profile of the CS (2 - 1) transition (dashed grey line in Figure 16) appears single peaked, with a maximum optical depth value of 105. The maximum optical depth values computed from the LOC + MCMC for these two molecules are unexpectedly high. This is probably due to the high found by the LOC + MCMC approach. This result also stresses the fact that, with the lack of enough constraints, the result found by LOC + MCMC may not be the ”correct” or a physically accurate result.
Appendix B Corner Plots
This Section introduces the corner plots for the LOC + MCMC modelling results presented in Section 4. In all of the corner plots, the top left subplot corresponds to the inner abundance values explored in the modelling. The second column is the outer abundance. The third corresponds to the radius separating the abundance regions. The fourth shows the velocity profile scaling, and the fifth the turbulent velocity dispersion. The rest of the panels show the correlation between the different variables. On top of each column the resulting model parameter value and its errors are displayed. The model parameter values are additionally listed in the main text in Table LABEL:modres.
Appendix C Additional LOC + MCMCs
Complementary models to those introduced in Section 4 are presented in the following Sections. The LOC + MCMC parameter values are summarised in Table LABEL:modresa.
C.1 CS + C34S modelling
In an attempt to improve the fit and parameter constraints presented in Section 4.2.1, we construct a new model which fits the CS and C34S (2 - 1) observed lines simultaneously. The idea behind this combined modelling is to have the same amount of variables explored with the MCMC with more observational constraints coming from two observed spectral lines instead of one. This model assumes that the main and rarer isotopologue share the same abundance profile, and the rarer isotopologue profile is scaled down by the 32S/34S = 22 ratio (Wilson & Rood, 1994). Thus, we have the same variables as for the CS and C34S separate models: ain, aout, , and . We test this combined modelling for both the non-extended and the extended models.
The combined non-extended model results for CS and C34S can be found in Figure 27. The model fit parameters are presented in Table LABEL:modresa as well as in the corner plot in the bottom panel in Figure 27. The CS (2 - 1) modelled transition fits the observations worse than the one shown in Section 4.2.1. The model slightly over-predicts the intensity of the blue peak with respect to the observations (top panel, Figure 27). The C34S (2 - 1) modelled transition fit is worse (Middle Panel, Figure 27) when compared to the individual C34S model in Section 4.1.1 by underproducing the red part of the line and overproducing the self-absorption dip. The resulting model parameters, and specifically , and show values close to the ones found for the CS individual model. The combined non-extended model mainly fits CS, and thus, the modelled parameters are closer to the CS individual fit model ones.
Notice that the corner plot of this combined fit also resembles the corner plot for the CS individual fit with really narrow distributions. Moreover, the value found in the combined fit of 0.00 also equals the prior distribution limit, being a sign that there may have been a problem with the exploration.
C.2 SO + 34SO modelling
With the same aim as for CS + C34S, we construct a new model to fit SO and 34SO (2,3 - 1,2) simultaneously. As for CS + C34S (Section C.1) the model for SO and 34SO assumes that the main and rarer isotopologue share the same abundance profile but the rarer isotopologue profile is scaled down by the 32S/34S = 22 ratio (Wilson & Rood, 1994). As the SO line fits the observations with the non-extended model, we try the SO + 34SO combined fit with the non-extended model.
The fit results can be found in Figure 28. The model parameters are presented in Table LABEL:modresa as well as in the corner plot in the bottom panel in Figure 28. The combined model fits SO and 34SO (Top and Middle Panels, respectively, Figure 28) similarly as the respective single models do (Figures 8 and 5, respectively). The and parameters are higher than for both the SO and 34SO individual fits. Both the and parameters are lower than the SO and 34SO individual models. Finally, the value is the same than for the individual fits.
C.2.1 HCO*+* Redaelli et al. 2022 reproduction.
Using the same modelling approach as Redaelli et al. (2022), which involves using the HCO*+* abundance profile from pyRate, we can reproduce their HCO*+* 1 - 0 transition fit (right panel in Figure 7, Redaelli et al., 2022), as shown in Figure 29. The model intensity increases by 14% in our model.
C.2.2 HCO*+* + HC17O*+* modelling
The combined HCO*+* + HC17O*+* (1 - 0) fit has been constructed similarly as the CS + C34S and the SO + 34SO pair fits. The abundance profile used is the same as the one described in Section 4.2.3 for HCO*+* and it is scaled down by [16O/17O] = 2044 (Penzias, 1981; Wilson & Rood, 1994) for HC17O*+. As for the individual HCO+* fit, we only test the extended fit.
The HCO*+* + HC17O*+* extended model fits are shown on the Top and Middle Panels, respectively, in Figure 30. The model parameters are presented in Table LABEL:modresa as well as in the corner plot in the bottom panel in Figure 30. The HCO*+* modelled spectrum does not reproduce the observed relative intensity of the blue and red components but reproduces the blue and red peak separation. The HC17O*+* modelled spectrum does not fit the observed line intensity nor its profile. The HCO*+* + HC17O*+* combined fit does not improve the constraints on the fractional abundance profile compared to the separate HCO*+* and HC17O*+* fits.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1André et al. (2014) André, P., Di Francesco, J., Ward-Thompson, D., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 27, doi: 10.2458/azu_uapress_9780816531240-ch 002 · doi ↗
- 2Bizzocchi et al. (2013) Bizzocchi, L., Caselli, P., Leonardo, E., & Dore, L. 2013, A&A, 555, A 109, doi: 10.1051/0004-6361/201321276 · doi ↗
- 3Bogey et al. (1981) Bogey, M., Demuynck, C., & Destombes, J. L. 1981, Chemical Physics Letters, 81, 256, doi: 10.1016/0009-2614(81)80247-3 · doi ↗
- 4Bogey et al. (1986) —. 1986, Chemical Physics Letters, 125, 383, doi: 10.1016/0009-2614(86)85177-6 · doi ↗
- 5Brünken et al. (2003) Brünken, S., Müller, H. S. P., Lewen, F., & Winnewisser, G. 2003, Physical Chemistry Chemical Physics (Incorporating Faraday Transactions), 5, 1515, doi: 10.1039/B 301657 A · doi ↗
- 6Caselli et al. (1999) Caselli, P., Walmsley, C., Tafalla, M., Dore, L., & Myers, P. 1999, Ap J, 523, L 165, doi: 10.1086/312280 · doi ↗
- 7Caselli et al. (2002 a) Caselli, P., Walmsley, C., Zucconi, A., et al. 2002 a, Ap J, 565, 331, doi: 10.1086/324301 · doi ↗
- 8Caselli et al. (2002 b) Caselli, P., Walmsley, C. M., Zucconi, A., et al. 2002 b, Ap J, 565, 344, doi: 10.1086/324302 · doi ↗